Uncertainty propagation of fatigue crack growth life commonly aims to provide the probability distribution of the lifespan needed for probabilistic damage tolerance analysis and for structural integrity assessment. This paper presents a novel methodology for efficiently estimating the parameters of the probability distribution of fatigue lifespan considering the Pearson distribution family. First, the full secondorder approach for expected value and variance prediction of probabilistic fatigue crack growth life is extended to predict higher order statistical moments of the underlying distribution. That is, the expected value (first raw moment) and the variance (second central moment) equations are complemented with the probabilistic formulations for the skewness and for the kurtosis (third and fourth central standardized moments, respectively). Then, from these moments, the Pearson distribution type is automatically determined. Finally, the parameters of the particular Pearson distribution type are estimated making the statistical moments of the constructed lifespan distribution match the first four prescribed moments predicted by the probabilistic equations. The validity of the proposed method is verified by a numerical example regarding the fatigue crack growth in a railway axle under random bending loading. It is proven that the probability density function of the lifespan is properly derived by the methodology, without knowing or assuming the output probability distribution beforehand. The methodology presented enables an efficient and an accurate quantification of the lifespan uncertainties via its probabilistic distribution. This probabilistic description of fatigue crack growth life can be subsequently used in reliability studies or in damage tolerance assessment.
Keywords: Probabilistic fatigue crack growth, uncertainty propagation, lifespan probability distribution, pearson distribution family, NASGRO
: crack depth / normal semiaxis
: final crack depth
: initial crack depth
: real roots of the quadratic equation in the denominator of the integral of the Pearson solution
: superscript above referring to the bending loading case
: superscript above referring to the bending plus interference loading case
: parameter of the crack growth equation in the Paris region
: random variables number
: crack growth rate
: scale parameter of the lognormal distribution
: Newman’s crack opening function
: probability density function for the standardized beta distribution
: probability density function for the standardized lognormal distribution
: probability density function of the random variable fatigue life
: first partial derivative of with respect to
: second partial derivative of with respect to and
: evaluation of at the mean value vector of
: step increment
: superscript above referring to the interference loading case
: index from 1 to random variables
: stress intensity factor (general)
: critical stress intensity factor for static unstable crack growth
, : maximum and minimum
: bending moment
: exponent of the crack growth equation in the Paris region
: steps number
: number of applied cycles
: parameter describing the sigmoidal shape of the crack growth equation in the threshold region
: Pearson probability density function
: parameter describing the sigmoidal shape of the crack growth equation in the toughness region
: stress ratio
: radial coordinate direction at the axle surface
: set of random variables
: general multivariate function
: shape parameters of the beta distribution
: square of skewness of
: kurtosis of
: skewness of
: gamma function
: stress intensity factor range
: threshold stress intensity factor range
: location parameter of the lognormal and the beta distributions.
: expected value of or first raw moment simply denoted
: third central moment of
: fourth central moment of
: stress (general)
: shape parameter of the lognormal distribution
: standard deviation of
: variance of also simply denoted
: probability density function for the standard normal distribution
FOSM: first order second moment
FEM: finite element method
SIF: stress intensity factor
MC: Monte Carlo
Pr. Eq.: probabilistic equation
PDF: probability density function
The fatigue crack growth mechanism is affected by many sources of variability such as the statistical variations of the mechanical properties of solid materials [1–4], random loading conditions [5–7] and uncertainties inherent to geometrical parameters [8]. Therefore, the stochastic nature of the fatigue crack growth phenomenon requires a proper uncertainty processing for reliability assessment and durability prediction. Uncertainty propagation methods in the fatigue life prediction commonly characterize the lifespan after a certain crack growth through a probability distribution. In consequence, probabilistic fatigue crack growth analyses aim to determine the probability distribution of the random response, lifespan, as a result of the randomness of the input variables. In areas such as damage tolerance analysis in the railway sector, predicting the probability distribution of fatigue crack growth life is becoming an important aspect. However, there is still a lack of harmonization of probabilistic procedures to define inservice inspection intervals in train axles for crack detection [9,10]. This particular application still deserves some developments via the extension of the current numerical methodologies.
Probabilistic fatigue crack growth is an extensive area of research. Numerous studies develop probabilistic models based on deterministic crack growth equations such as the Paris’ law [11]. For instance, some probabilistic approaches founded on the Paris law are presented in [12–21], being the model most frequently used. Over recent years, the first advances considering the complete crack growth curve have been made because of the increasing importance of taking into account the early stages of crack growth for an adequate estimation of the lifespan. The state of the art equation in the field of the fatigue crack propagation problem is a modified version of the Paris’ law referred to as NASGRO equation [22]. The NASGRO equation is commonly applied in the evaluation of components as in the case of railway axles [18,23–25]. Additionally, the use of the firstorder secondmoment (FOSM) probabilistic method [26–29] has enabled some new perspectives in the fatigue crack growth analysis of components. Some examples are the fatigue life estimation according to the Paris' law [13,14], the crack nucleation stage analysis based on the CoffinManson model [30], the multiaxial fatigue assessment based on the virtual strain energy damage model of Liu [31,32] and especially the fatigue crack propagation based on the NASGRO equation in [25] where an extended version of the FOSM to full secondorder approach is developed for predicting statistical moments of the fatigue lifetime in a railway axle under random loading.
The probabilistic fatigue crack growth strategies above provide certain statistics that describe the response distribution, for instance, mathematical moments like the expected value and the variance. Then, the problem under investigation is further analyzed according to these statistics. One such interesting analysis is to construct the whole probability density function (PDF) of fatigue crack growth life, namely, the uncertainty propagation in fatigue crack growth life. The PDF can be used to describe the relative likelihood that the value of the lifespan would be equal to a particular value. Generally, the random output variable is assumed to follow a common probability distribution, such as the normal or the lognormal distribution, and subsequently the parameters of the distribution are estimated considering the statistics of the response as constraints. An illustrative example of the previous strategy is presented in [25] where the full secondorder approach provides the expected value and the variance of the random output variable lifespan and then the normal and the lognormal PDFs are constructed according to these two statistics of the lifespan. The discussion in the example mentioned above highlights some drawbacks regarding the need to assume in advance a probability distribution of two parameters, since only the expected value and the variance of the underlying lifespan probability distribution are calculated, and also that the accuracy or similarity between the constructed PDFs and the results of the Monte Carlo method (MC) can be improved. A general methodology that considers the statistics of the response and automatically provides an appropriate PDF of the response distribution will be more useful than a predefined one. It would thus be of interest to explore some alternatives to overcome the limitations aforementioned to contribute to a better knowledge of the distribution of fatigue life.
As shown in the previous literature, the scientific community has studied the uncertainty propagation in fatigue crack growth life for many years. These investigations are usually intended to provide the probability distribution of the lifespan which, as is the case of train axles, is a fundamental aspect for the maintenance planning definition to keep the probability of failure as small as possible, as far as it is economically viable. Addressing the problem in an efficient and precise manner is, therefore, a challenge. The construction of a PDF under moment constraints has historically been of great interest. In some cases, the underlying probability distribution manifests asymmetry or presents a characteristic heaviness of the tails relative to the rest of the distribution or a combination of both. Therefore the need to manage the shape of the constructed PDF arises. The shape of the constructed distribution can be handled by means of controlling its skewness and kurtosis moments. Early contributions considered this topic, devising a system of distributions, the Pearson distribution family [33], which has the appeal of encompassing several wellknown distributions as members. This distribution has a rich flexibility in shape, covering a wide skewnesskurtosis region. Furthermore, the location, scale and shape parameters of a particular Pearson distribution can be estimated in terms of the first four prescribed moments, that is, the expected value, first raw moment; the variance, second central moment; the skewness, third central standardized moment, and the kurtosis, fourth central standardized moment. As a consequence, the Pearson family of distributions has been used for modelling purposes by those engaged in the probabilistic analysis of structures. However, an overall treatment of the uncertainty propagation in fatigue crack growth life based on the NASGRO equation combined with the versatile Pearson distribution family is missing. To the best of the authors’ knowledge, the work presented in this paper is the first attempt to develop an effective, efficient and practical uncertainty propagation approach that is capable of predicting the first four moments of the lifespan distribution according to the NASGRO equation and using them to construct a probability density function based on the Pearson distribution family.
The objective of this work is to providea novel uncertainty propagation methodology for efficiently estimating the parameters of the probability distribution of fatigue lifespan considering the Pearson distribution family. For this purpose, the pertinent Pearson distribution type is automatically determined depending on the stochastic moments of the underlying lifespan distribution according to the NASGRO equation. Subsequently, the parameters of the distribution are calculated making the statistical moments of the constructed lifespan distribution match the first four prescribed moments. To estimate these first four required moments in advance, the full secondorder approach for expected value and variance prediction presented in [25] is extended to predict the skewness and the kurtosis of the underlying lifetime distribution. As a result, the methodology developed enables an efficient estimation of the parameters of the particular Pearson distribution type and, in the meantime, provides the first four moments of the underlying lifetime distribution. The combination of both, contributes to overcome the drawbacks previously highlighted, avoiding the need to assume a probability distribution of two parameters and also increasing the accuracy of the PDF constructed. The proposed uncertainty propagation methodology is illustrated and discussed in detail through an application example specifically oriented on the fatigue crack growth in a metallic railway axle.
The paper is organized as follows. Section 2 presents the methods used for the uncertainty propagation methodology. Firstly, the full secondorder approach for the skewness and kurtosis moments of functions of random variables is described in detail. It also includes the probabilistic expressions regarding its integration with the NASGRO model. Secondly, an overview of the Pearson distribution family is given and the estimation of the distribution parameters from prescribed moments is summarized. Section 3 illustrates and discusses the proposed uncertainty propagation procedure by means of an application example in a metallic railway axle. The results of the full secondorder approach are compared with those obtained by the Monte Carlo method in terms of skewness and kurtosis. Furthermore, the quality of PDF constructed with the Pearson distribution in comparison with the histogram from the Monte Carlo method is discussed. Finally, Section 6 concludes the main findings in this paper.
The objective of the uncertainty propagation methodology is to efficiently estimate the parameters of the probability distribution of a random variable of interest. To establish a general methodology, a key point is to consider the Pearson distribution family because of the following two reasons combined. Firstly, it is a rich family that is able to adjust a broad range of distribution shapes and therefore it can closely represent a portion of the reality. Secondly, the Pearson distribution family is easy to use in practice because it has mathematical characteristics that make the different distribution parameters be expressed as function of the first four moments.
To obtain the first four required moments of the underlying distribution, the full secondorder approach for expected value and variance prediction of probabilistic fatigue crack growth life is extended to predict higher order statistical moments of the underlying distribution. That is, the expected value (first raw moment) and the variance (second central moment) probabilistic equations (Pr. Eq.) are complemented with the Pr. Eq. for the skewness and for the kurtosis (third and fourth central standardized moments, respectively). The reasoning followed for the expected value and the variance estimation is used to derive the Pr. Eq. for the skewness and the kurtosis prediction.
The full secondorder approach for the moments of functions of random variables presented in [25] enables the prediction of the expected value and the variance of the lifespan of interest. Here, the approach is extended to predict the skewness and the kurtosis of the probabilistic fatigue crack growth life. This section recalls the definition of skewness and kurtosis, then derives the full secondorder approximation for the third and fourth central moments of a general function and finally provides the probabilistic formulations for propagating the third and fourth central moments through the fatigue crack growth NASGRO model.
The third central moment is a measure of the distribution symmetry or lack of symmetry. Any symmetric distribution will have a third central moment of zero. The third central moment divided by the standard deviation to the power of three is referred to as normalized or standardized third central moment and commonly known as skewness. The skewness is denoted by as it is shown in Eq. (1)

(1) 
As the third central moment has units of the random variable to the power of three, after normalization, skewness is dimensionless. Usually, the square of the skewness is denoted by as it is expressed in Eq. (2)

(2) 
The fourth central moment is a measure of the combined weight or heaviness of the tails relative to the rest of the distribution. As the variance, it measures the spread of the data but is more dependent on the behavior of outliers in the tails. The fourth central moment divided by the standard deviation to the power of four is known as normalized or standardized fourth central moment and generally referred to as kurtosis. The kurtosis is frequently denoted by as it is presented in Eq. (3), being dimensionless for the same reason mentioned above

(3) 
The skewness and the kurtosis are calculated from the third and fourth central moments as indicated in the previous definitions. Following the reasoning presented in [25] to predict the expected value and the variance, the third central moment approximation is derived and as a result the skewness can be calculated. The full secondorder approximation for the third central moment of a general function is given by Eq. (4)

(4) 
Note that indicates second central moment , indicates the third central moment , indicates the fourth central moment and so on. In other words the central multivariate moments of continuous random variables are denoted with consecutive indexes, two for the second central moment, three indexes for the third central moment, four indexes for the fourth central moment etc.
Likewise, following the reasoning presented in [25] the fourth central moment approximation is derived and as a result the kurtosis can be calculated. The full secondorder approximation for fourth central moment of a general function is given by Eq. (5)

(5) 
At this point, to bridge the gap between the full secondorder approach presented and the fatigue crack growth life prediction, the probabilistic formulations for propagating the third and fourth central moments through the fatigue crack growth NASGRO model are deduced. To begin with, the NASGRO model described from a probabilistic point of view in [25] is adopted. This probabilistic form is adapted to obtain the distribution of the number of load cycles to reach a given crack depth . It basically starts from the original NASGRO crack propagation rate , isolates and using the discretised version for every crack growth increment leads to Eq. (6)

(6) 
where, is the stress ratio, is the stress intensity factor (SIF) range from the maximum and minimum , i.e. and , is the crack opening function, is the threshold stress intensity factor range, is the critical stress intensity factor and , , , and are material empirically derived constants. For a more detailed description of the previous equation refer to [22,25].
The third central moment of the fatigue life is obtained by applying the third central moment definition to the total number of cycles which is the summation of the life increments up to a final crack depth . Using the formula for the third central moment of the sum of random variables leads to Eq. (7)

(7) 
In the same manner as above, the fourth central moment of the fatigue life is obtained by using the formula for the fourth central moment of the sum of random variables giving Eq. (8)

(8) 
In Eqs. (7) and (8) the moments and of are obtained by applying the full secondorder approximation developed to obtain the stochastic central moments of third and fourth order Eqs. (4) and (5) respectively, on the function Eq. (6). Eqs. (7) and (8) are referred to as the probabilistic NASGRO equations for the third and fourth central moments respectively. In summary, for each crack growth increment the third and fourth central moments of the fatigue lifetime increment are calculated by applying the full secondorder approximation method on the discretised version of the NASGRO equation, requiring the first to eighth order moments of the random input variables and the first and second partial derivatives of the NASGRO equation with respect to the random input variables.
Many statistical analyses consider the normal distribution when there is not more information available about the underlying probability distribution, however, in some cases, this assumption does not reflect the reality. Among the different statistical distributions that can be chosen, the Pearson distribution family is considered in the methodology as it is a versatile family that covers a broad range of distribution shapes and therefore it is able to model a wide range of data accurately. Moreover, it enables the expression of the distribution parameters as a function of the first four moments of the distribution without a priori hypotheses and thus it is an efficient and objective procedure.
The Pearson distribution is a family of continuous probability distributions [33]. A Pearson probability density function is defined to be any valid solution to the first order linear differential equation presented in Eq. (9)

(9) 
with

where, is a location parameter, is the square of the skewness, is the kurtosis, and is the second central moment of the distribution. The solution to the above differential equation is shown in Eq. (10)

(10) 
The variety of solutions differ in the values of the coefficients , , and . These solutions are classified into different Pearson distribution types depending on the quantities and . In other words, the type of distribution to which the data belong is entirely determined by the square of the skewness presented in Eq. (2) and the kurtosis expressed in Eq. (3). The Pearson distribution types correspond to common probability distributions. The following types and their common distribution associated arise: type I (beta distribution), type II (symmetrical beta distribution), type III (gamma distribution), type IV (Cauchy distribution), type V (inversegamma distribution), type VI (beta prime distribution), type VII (Student's tdistribution) and the limit of type I, III, IV, V (normal distribution).
Once the Pearson distribution type is determined, the two, three or four parameters of the particular distribution type can be calculated as a function of the expected value, variance, skewness and kurtosis, i.e. from the first four moments. The formulas to calculate the parameters for each type of Pearson distribution are given in [34]. As an example, the formulas for the Pearson distribution type I that corresponds to the beta distribution are included here. A random variable following a beta distribution with the shape parameters and , can be further parametrized as with the scale parameter and with a location parameter . The shape parameters and are calculated according to the Pearson distribution family as expressed in Eq. (11)

(11) 
where are the real roots of the quadratic function in the denominator of Eq. (10). Hence, they are calculated as shown in Eq. (12), and additionally they set the scale parameter

(12) 
Finally, the location parameter calculated according to the Pearson distribution family is shown in Eq. (13)

(13) 
where is the expected value of the distribution of the original variable to model, i.e. its first raw moment.
The aim of this example is to illustrate the uncertainty propagation methodology described in the present paper which extends the full secondorder approach for expected value and variance prediction to predict the skewness and the kurtosis of the fatigue lifetime based on NASGRO model and estimates the parameters of a particular Pearson distribution type based on the first four prescribed moments. The methodology was applied to the example 1 in [25] to verify that the proposed improvements contribute to a better knowledge of the distribution of fatigue life. Briefly, the numerical example dealt with the fatigue crack growth in the railway axle shown in Figure 1 under random bending moment. The axle was 173 mm in diameter and it was made of EA1N steel defined in the EN 13261 standard [35].
Figure 1. General axle view of a nonpowered wheelset with a postulated crack in the Ttransition 
A semicircular initial crack of 2 mm was postulated at the Ttransition, as indicated in the crosssection of the Figure 1. The crack grows describing the crack propagation direction defined according to the radial coordinate system at the axle surface in Figure 1. The crack kept a semielliptical shape while growing up to a final crack depth of 50 mm. The fatigue crack growth material parameters used in the NASGRO model were those collected in [25].
For the present study, the following two types of loads were considered: the bending moment loading mainly due to the vehicle weight and cargo and the pressﬁt loading produced by the wheel mounting with interference. The bending moment in the railway axle was assumed as a random input variable normally distributed with a standard deviation equal to the 5% of the mean value. The parameters of the bending moment distribution were: mean value and variance . The bending moment level selected corresponded to the highest load amplitude in the spectrum of a railway axle which is caused by the weight of the wagon and cargo, 22.5 tonnes per axle, plus additional forces, generated when the train goes through curved track, over crossovers, switches, rail joints, braking efforts, etc. Additionally, the wheel was pressfitted with 0.286 mm interference in diameter. The reference bending stress amplitude for the mean value of bending moment and the interference stress normal to the crack surface needed for the stress intensity factor and evaluation were calculated via the finite element method (FEM). The stresses at the interference , at the reference bending moment and at the combination of both are shown in Figure 2. In this adapted version of the original figure in [25], a schematic representation of the axial stress distributions originated at the Ttransition in a solid axle under a pure pressfit with interference and under a pure bending load case is also included for clarity.
Figure 2. Schematic representation and calculated stress distributions at the Ttransition in axial direction at the pressfit with interference, at the bending and at the superposition of both load cases 
In the example, there was a single source of variability from the bending moment that induced variability on the bending stress and therefore there was variability in the combination of the bending with the interference stress. The stress intensity factors and depend on the maximum and minimum stresses and in consequence they were also stochastic. The combination of stresses was obtained considering the fact that the railway axle rotates as the train moves. Consequently, every revolution the crack goes through all possible positions with respect to the applied moment since the bending load always acts in the vertical plane and hence the resulting load type is rotary bending. The maximum and minimum stresses were related to each other as a result of the stresses symmetry. This fact together with the cyclic rotary character of the load implies that the random variables and were related to each other too. These two random variables used to describe the input variability were correlated variables so they were not independent. It is important to emphasize that and are explicitly involved in the NASGRO equation.
The full secondorder approach presented was applied to calculate the third moment and the fourth central moments of at every crack depth, Eq. (6), with the two random input variables and . Then, the skewness and the kurtosis of the fatigue life were obtained from the moments, providing a continuous result along the crack depth . The results provided by the proposed methodology were compared with the results of 10 000 MC simulations. The skewness and kurtosis results provided by the MC were obtained by further processing the MC simulations performed in [25].
To check the accuracy of the method in terms of skewness, the history values provided by the MC and the Pr. Eq. are compared in Figure 3.
Figure 3. Skewness of history values provided by the Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr.Eq.) 
The values of the skewness curves for 550 mm crack depths are listed in Table 1.
MC  Pr. Eq.  Pr. Eq.MC  

[mm]  [dimensionless]  [dimensionless]  [%] 
5  0.844  0.781  7.52% 
10  0.777  0.722  7.14% 
20  0.733  0.682  6.91% 
30  0.712  0.664  6.70% 
40  0.698  0.653  6.47% 
50  0.689  0.646  6.24% 
The skewness computed from the MC simulations distribution is considered as the framework of reference for comparison. It is observed that the skewness values are always positive, indicating that the right tail, in a sense, is longer and or heavier than the left one. Additionally, it can be seen that as the crack depth increases the skewness decreases reaching a value greater than 0.6 in both calculations. The error in the Pr. Eq. ranges from 7.5% to 6.2%, reducing as the crack depth increases. The Pr. Eq. provides a slightly lower magnitude than the MC for every crack depth.
To check the accuracy of the method in terms of kurtosis, the history values provided by the MC and the Pr. Eq. are compared in Figure 4.
Figure 4. Kurtosis of history values provided by the Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr.Eq.) 
The values of the kurtosis curves for 550 mm crack depths are listed in Table 2.
MC  Pr. Eq.  Pr. Eq.MC  

[mm]  [dimensionless]  [dimensionless]  [%] 
5  4.33  3.84  11.33% 
10  4.13  3.73  9.53% 
20  4.00  3.66  8.38% 
30  3.94  3.63  7.83% 
40  3.90  3.61  7.47% 
50  3.88  3.60  7.22% 
Once again, the kurtosis provided by the MC simulations is considered as the framework of reference. It is observed that the larger the crack the lower kurtosis and, in both calculations, a stabilization in magnitude is reached. The error in the probabilistic equation ranges from 11% to 7%, decreasing as the crack depth increases. The maximum difference of the Pr. Eq. calculates a kurtosis about 0.5 smaller.
At this point, the first part of the methodology presented in Section 2.1 is completed. As a result, the first to fourth moments of the underlying lifespan probability distribution are available after applying the full secondorder approach. The expected value, first raw moment, and the variance, second central moment are taken from [25], whereas the skewness and the kurtosis, third and fourth central standardized moments, respectively, have been calculated above. Recapitulating, the moments of fatigue life provided by the probabilistic NASGRO equations for a crack depth equal to 50 mm are enclosed in Table 3.
Pr. Eq.  Pr. Eq.  Pr. Eq.  Pr. Eq.  

[mm]  [Cycles]  [Cycles]  [dimensionless]  [dimensionless] 
50  398 853  56 103  0.646  3.60 
Note that the square root of the variance, i.e. the standard deviation , was reported instead of the variance . This is because the standard deviation is more meaningful to interpret since it is expressed in units of the original variable.
The next part of the methodology addresses the problem of constructing a PDF from prescribed moments as described in Section 2.2. In this case, it is possible to construct probability distributions with more than two parameters given that the first four moments are available. To check the benefits of the uncertainty propagation methodology presented in terms of life distributions, three scenarios were considered: (i) the lifespan was assumed to be normally distributed; (ii) the lifespan was assumed to be lognormally distributed; (iii) the Pearson distribution family was used to model the lifespan, thus avoiding the need of assuming a distribution in advance. The first two scenarios were discussed in [25], however, the latter one was not investigated as it requires the first four moments and only the expected value and the variance were calculated in the previous procedure. In the present work the third scenario considering the Pearson distribution family is also explored. The Pearson family is compared with the other two scenarios to verify that the proposed improvements contribute to a better knowledge of the distribution of fatigue life. For that purpose, the three probability density functions were constructed as described below from the moments in Table 3.
In case (i), the probability density function for the lifespan is assumed to be normally distributed, and therefore it is derived from the wellknown probability density function for the standard normal distribution . Specifically, the probability density function is parameterized in terms of a location parameter which is directly the expected value of the fatigue life , and a scale parameter which is directly the square root of the variance, i.e. the standard deviation of the random output variable fatigue life . As a result, the lifespan distributed normally is characterised by the probability density function , where and are fixed parameters.
In case (ii), the probability density function for the lifespan is assumed to be lognormally distributed, and therefore it is derived from the standardized form for the lognormal probability density function with as a shape parameter. Commonly, the probability density function is further parameterized in terms of a location parameter and a scale parameter . The location parameter was set equal to zero, hence only two parameters remain to be estimated. The physical meaning of the location parameter is realistic since it indicates that crack growth occurs after any given cycle and so the support of the distribution is nonnegative. The shape and the scale parameters of the lognormal are functions of the expected value and the variance of the fatigue life . As a consequence, the lifespan distributed lognormally is described by the probability density function where , and are fixed parameters.
In case (iii), the Pearson distribution type was automatically determined based on the quantities and as it is described in Section 2.2. In this example, the procedure leads to the Pearson type I that corresponds to the beta distribution. Afterwards, the parameters of the beta distribution were estimated making the statistical moments of the constructed lifespan distribution match the first four prescribed moments predicted by the Pr. Eq. In this case, the probability density function for the lifespan betadistributed is derived from the standardized form for the beta probability density function in Eq. (14) with , as shape parameters and where is the gamma function . Once again, the probability density function is further parameterized by introducing two parameters representing the location and the scale. The shape parameters and , the scale and the location parameter were calculated by using the closedform expressions derived in Section 2.2 that express the parameters as function of the first four moments of the fatigue life distribution. Recall that in these expressions the expected value or first raw moment is denoted as . As a result, the lifespan following a beta distribution is characterized by the probability density function given in Eq. (15)

(14)  

(15) 
Note that in the probability density function the , , and are fixed parameters.
The resulting parameters for the normal, the lognormal and the beta PDFs are collected in Table 4.
Probability Distr.  Shape  Location  Scale 

Normal    398 853  56 103 
Lognormal  0.140  0  394 964 
Pearson type I (Beta)  8.49 ,  232 076  4 825 029 
237.33 
The normal distribution, the lognormal distribution, the beta distribution and the MC histogram of the fatigue life for a crack depth equal to 50 mm are compared in Figure 5.
Figure 5. Histogram of fatigue life provided by the Monte Carlo (MC) and PDFs of the normal, the lognormal and the beta distributions constructed from the moments provided by the probabilistic NASGRO equations (Pr.Eq.) for 50 mm crack depth 
The histogram of fatigue life provided by the MC simulations is taken as a benchmark for the comparison of the three PDFs constructed. The differences between the normal distribution and the MC histogram highlight the degree of nonnormality and nonsymmetry of the underlying lifespan distribution. The lognormal distribution includes a certain degree of asymmetry and therefore it is a more convenient choice over the symmetrical normal distribution. It can be observed that the lognormal is advantageous to the normal distribution in describing the lower and higher tails of the lifespan, although the degree of asymmetry in the lognormal is not directly defined. On the other hand the beta distribution was determined from the Pearson distribution family according to the quantities related to shape and . It can be observed that the beta distribution agrees well with the MC histogram for all of the lifespan range, including the tails and also the peak. The superiority of the beta distribution over the normal and the lognormal distributions to represent the MC results is clear, especially when describing the lower tail of the distribution of lives, which is of great importance in reliability and in damage tolerance assessment. The behaviour of the normal and the lognormal distributions in the lower tail region slightly underestimate the lifespan. This aspect is clearly observed in the inset of Figure 5 which shows a zoom of the lower tail. The better tail performance of the beta distribution can be attributed to the fact that the first four moments of the lifespan are matched when estimating the parameters of the Pearson distribution. That is to say that the PDF using the Pearson distribution family is sensitive to variations in skewness and kurtosis. Besides, the Pearson distribution type is selected depending on the moments that are related to the shape of the underlying distribution, and afterwards the corresponding probability distribution parameters are calculated. On the contrary, the construction of the normal and lognormal distributions do not take into account the skewness or kurtosis changes given that their third and fourth moments are not matched to the calculated lifespan moments. In fact, the normal distribution has skewness and a kurtosis equal to inherently, and the skewness and the kurtosis of the lognormal distribution are determined indirectly once the distribution is defined in terms of the expected value and variance of the lifespan . Consequently, the quality of the normal and lognormal PDFs is influenced by the similarity between them and the actual lifespan distribution, which is not always known. In other words, the quality can only be good if the lifespan is close to the selected distribution. It must be emphasized that in some applications the normal or the lognormal assumption may be enough, however, the Pearson distribution family still offers a perceptible improvement as the goodness and quality of the PDF constructed is not compromised by an a priori selection or assumption of a probability distribution.
Additionally, a KolmogorovSmirnov test [36] was used to measure quantitatively the goodness of fit between the constructed probability distributions and the MC results. The test statistic quantifies the distance between the cumulative distribution functions of the constructed distribution and the cumulative MC histogram of reference, and therefore the lowest value indicates the best representation of the underlying distribution. The test statistics obtained were 0.052, 0.027 and 0.017 for the normal, lognormal and beta distributions, respectively. These values confirm the advantage of the beta distribution over the normal and the lognormal distributions to represent the MC results.
The following outcomes are drawn from the analysis of the results obtained by applying the uncertainty propagation methodology:
This paper presents a novel uncertainty propagation methodology for efficiently estimating the probability distribution of the fatigue crack growth lifetime. It implements the Pearson distribution family and takes advantage of the statistical moments of the lifetime predicted via existing and innovative approaches based on NASGRO model. In the present work the full secondorder approach for expected value and variance prediction of probabilistic fatigue crack growth life presented in [25] is revisited and extended to predict higher order statistical moments of the underlying distribution. Based on that approach, the probabilistic equations for the skewness and kurtosis of the fatigue crack growth lifetime are provided. These two moments along with the expected value and variance are used to estimate the parameters of the Pearson distribution to represent the underlying lifespan distribution. The probability distribution constructed from the first four prescribed moments is helpful to describe the fatigue crack growth phenomenon under stochastic conditions such as under a random bending moment loading.
The present uncertainty propagation methodology is found to be advantageous compared to the previously applied in [25]. The main advantages that are worth mentioning include: (i) the ability to automatically determine the pertinent distribution shape using a simple and uniform criterion for choosing the corresponding Pearson distribution type and thus avoiding the need to assume a probability distribution in advance; (ii) the versatility of constructing probability distributions with more than two parameters if required given that in the current methodology the skewness and the kurtosis are calculated in addition to the expected value and variance predicted in the initial approach; (iii) the significant improvement in the accuracy of the probability density function constructed in the analysis of the fatigue crack growth in a railway axle under random bending loading, especially in the lower tail of the distribution of life which is of notable interest in reliability assessment. It demonstrates that in practical engineering it is relevant to determine the third and fourth central moments for a precise description of the probabilistic fatigue crack growth life.
In general, the uncertainty propagation methodology presented has a positive and comprehensive effect in engineering applications such as in the probabilistic life prediction. The lifespan probability distribution provided has the potential to be used in probabilistic damage tolerance design and in the integrity assessment of components and structures. Therefore, the novel uncertainty propagation methodology in this paper is expected to be an asset in a broad range of engineering problems dealing with random variables.
The authors acknowledge the Spanish Ministry of Economy, Industry and Competitive through the National Programme for Research Aimed at the Challenges of Society that financially supported the project RTC201648134.
[1] Cervello S. Fatigue properties of railway axles: New results of fullscale specimens from Euraxles project. Int J Fatigue, 86:2–12, 2016. https://doi.org/10.1016/j.ijfatigue.2015.11.028.
[2] Novosad M., Fajkoš R., Řeha B., Řezníček R. Fatigue tests of railway axles. Procedia Eng, 2:2259–68, 2010. https://doi.org/10.1016/j.proeng.2010.03.242.
[3] Beretta S., Carboni M. Variable amplitude fatigue crack growth in a mild steel for railway axles: Experiments and predictive models. Eng Fract Mech, 78:848–62, 2011. https://doi.org/10.1016/j.engfracmech.2010.11.019.
[4] Mädler K., Geburtig T., Ullrich D. An experimental approach to determining the residual lifetimes of wheelset axles on a fullscale wheelrail roller test rig. Int J Fatigue, 86:58–63, 2016. https://doi.org/10.1016/j.ijfatigue.2015.06.016.
[5] UIC B 169 RP 36:201312 Defect Tolerance Concept (DTC) Permitted defects and surface properties for axles under maintenance. UIC Standard, 2013.
[6] Pokorný P., Náhlík L., Hutař P. Comparison of different load spectra on residual fatigue lifetime of railway axle. Procedia Eng, 74:313–6, 2014. https://doi.org/10.1016/j.proeng.2014.06.269.
[7] Pokorný P., Hutař P., Náhlík L. Residual fatigue lifetime estimation of railway axles for various loading spectra. Theor Appl Fract Mech, 82:25–32, 2016. https://doi.org/10.1016/j.tafmec.2015.06.007.
[8] Traupe M., Landaberea A. EURAXLES  A global approach for design, production and maintenance of railway axles: WP2  development of numerical models for the analysis of railway axles. Mater Werkst, 48:687–98, 2017. https://doi.org/10.1002/mawe.201600570.
[9] Zerbst U., Beretta S., Köhler G., Lawton A., Vormwald M., Beier H.Th., et al. Safe life and damage tolerance aspects of railway axles – A review. Eng Fract Mech, 98:214–71, 2013. https://doi.org/10.1016/j.engfracmech.2012.09.029.
[10] Zerbst U., Klinger C., Klingbeil D. Structural assessment of railway axles – A critical review. Eng Fail Anal, 35:54–65, 2013. https://doi.org/10.1016/j.engfailanal.2012.11.007.
[11] Paris P.C., Erdogan F. A critical analysis of crack propagation laws. J Basic Eng, 85:528–33, 1963. https://doi.org/10.1115/1.3656900.
[12] Akama M., Ishizuka H. Reliability analysis of Shinkhansen vehicle axle using probabilistic fracture mechanics. JSME Int J Ser Mech Mater Eng, 38(3):378–83, 1995. https://doi.org/10.1299/jsmea1993.38.3_378.
[13] Bea J.A. Simulación del crecimiento de grietas por fatiga aleatoria mediante elementos probabilistas. PhD Thesis, Universidad de Zaragoza, 1997.
[14] Bea J.A., Doblaré M., Gracia L. Evaluation of the probability distribution of crack propagation life in metal fatigue by means of probabilistic finite element method and Bmodels. Eng Fract Mech, 63:675–711, 1999. https://doi.org/10.1016/S00137944(99)000533.
[15] Bea J.A., Doblaré M., Gracia L. Fiabilidad de elementos metálicos en crecimiento de grieta por fatiga aleatoria mediante elementos finitos probabilistas y modelos B. Rev Int Métod Numér para Cálculo Diseño en Ing, 15(1):85–112, 1999. https://www.scipedia.com/public/Bea_et_al_1999a.
[16] Hillmansen S., Smith R.A. The management of fatigue crack growth in railway axles. Proc Inst Mech Eng Part F J Rail Rapid Transit, 218:327–336, 2004. https://doi.org/10.1243/0954409043125879.
[17] Hong Y.J., Xing J., Wang J.B. A secondorder thirdmoment method for calculating the reliability of fatigue. Int J Press Vessels Pip, 76:567–70, 1999. https://doi.org/10.1016/S03080161(99)000137.
[18] Náhlík L., Pokorný P., Ševčík M., Fajkoš R,, Matušek P,, Hutař P. Fatigue lifetime estimation of railway axles. Eng Fail Anal, 73:139–57, 2017. https://doi.org/10.1016/j.engfailanal.2016.12.014.
[19] Wang L., Liang J., Yang Y., Zheng Y. Timedependent reliability assessment of fatigue crack growth modeling based on perturbation series expansions and interval mathematics. Theor Appl Fract Mech, 95:104–18, 2018. https://doi.org/10.1016/j.tafmec.2018.02.010.
[20] Yasniy O., Lapusta Y., Pyndus Y., Sorochak A., Yasniy V. Assessment of lifetime of railway axle. Int J Fatigue, 50:40–6, 2013. https://doi.org/10.1016/j.ijfatigue.2012.04.008.
[21] Zhu S.P., Liu Q., Huang H.Z. Probabilistic modeling of damage accumulation for fatigue reliability analysis. Procedia Struct Integr, 4:3–10, 2017. https://doi.org/10.1016/j.prostr.2017.07.012.
[22] Forman R.G., Mettu S.R. Behavior of surface and corner cracks subjected to tensile and bending loads in Ti6Al4V alloy. ASTM STP 1131 Am Soc Test Mater Phila PA, 519–46, 1990.
[23] Beretta S., Carboni M. Experiments and stochastic model for propagation lifetime of railway axles. Eng Fract Mech, 73:2627–2641, 2006. https://doi.org/10.1016/j.engfracmech.2006.04.024.
[24] Beretta S., Villa A. A RV approach for the analysis of fatigue crack growth with NASGRO equation. Proceedings of the 4th International ASRANet Colloquium, 1–7, 2008.
[25] Mallor C., Calvo S., Núñez J.L., RodríguezBarrachina R., Landaberea A. Full secondorder approach for expected value and variance prediction of probabilistic fatigue crack growth life. Int J Fatigue, 133:105454, 2020. https://doi.org/10.1016/j.ijfatigue.2019.105454.
[26] Mayer M. Die Sicherheit der Bauwerke und ihre Berechnung nach Grenzkräften anstatt nach zulässigen Spannungen. Berlin, Springer Verlag, 1926.
[27] Mayer M. La seguridad en las construcciones y su cálculo aplicando los esfuerzos límite en lugar de las tensiones admisibles = Safety in constructional works and its design according to limit states instead of permissible stresses. Madrid, INTEMAC, 1975.
[28] Cornell C.A. A Probabilitybased structural code. ACI J Proc, 66:974–85, 1969. https://doi.org/10.14359/7446.
[29] Cornell C.A. Structural safety specifications based on secondmoment reliability analysis. IABSE Symposium on Concepts of Safety of Structures and Methods of Design, 1969. https://doi.org/10.5169/seals5948.
[30] Núñez J.L. Análisis del fenómeno de la fatiga en metales en etapa de nucleación mediante la utilización de modelos estadísticos de daño acumulado y elementos finitos probabilistas. PhD Thesis, Universidad de Zaragoza, 2003.
[31] Calvo S. Determinación de la probabilidad de fallo en componentes métalicos sometidos a estados multiaxiales de tensión mediante la utilización de elementos finitos probabilistas y modelos estadísticos de daño acumulado. PhD Thesis, Universidad de Zaragoza, 2008.
[32] Calvo S., Canales M., Gómez C., Valdés J.R., Núñez J.L. Probabilistic formulation of the multiaxial fatigue damage of Liu. Int J Fatigue, 33:460–5, 2011. https://doi.org/10.1016/j.ijfatigue.2010.10.003.
[33] Pearson K. X. Contributions to the mathematical theory of evolution.—II. Skew variation in homogeneous material. Philos Trans R Soc Lond A, 186:343–414, 1895. https://doi.org/10.1098/rsta.1895.0010.
[34] Johnson N.L., Kotz S., Balakrishnan N. Continuous univariate distributions. WileyInterscience, Vol. 1, 2nd edition, 784 pages, 1994.
[35] EN 13261:2009+A1:2010: Railway applications – Wheelsets and bogies – Axles – Product requirements. European Committee for Standardization, 2010.
[36] Haldar A., Mahadevan S. Probability, reliability, and statistical methods in engineering design. John Wiley, 320 pages, 1999.
Published on 16/07/20
Accepted on 06/07/20
Submitted on 12/03/20
Volume 36, Issue 3, 2020
DOI: 10.23967/j.rimni.2020.07.004
Licence: CC BYNCSA license
Are you one of the authors of this document?