Full second-order approach for expected value and variance prediction of probabilistic fatigue crack growth life

International Journal of Fatigue 133 (2020) 105454

Contents lists available at ScienceDirect

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

Full second-order approach for expected value and variance prediction of probabilistic fatigue crack growth life

T

C. Mallora, , S. Calvoa, J.L. Núñeza, R. Rodríguez-Barrachinaa, A. Landabereab ⁎

a b

Area of Research and Development, Instituto Tecnológico de Aragón, María de Luna, 8-50018 Zaragoza, Spain CAF S.A., J.M. Iturrioz, 26-20200 Beasain, Spain

ARTICLE INFO

ABSTRACT

Keywords: Probabilistic Fatigue crack growth Fatigue life prediction NASGRO Statistical moments

The NASGRO crack growth equation for fatigue life estimation is ordinarily used in its deterministic standard form. This paper presents a new probabilistic formulation for fatigue crack propagation based on the NASGRO equation providing a stochastic approach for predicting statistical moments of fatigue lifetime. The methodology approximates the expected value and variance, of the output random variable fatigue lifetime (). These moments are obtained from the approximation via the Taylor series up to the quadratic terms, full second-order, of the NASGRO equation with respect to the random input variables taken into account. The validity of the proposed method is verified by two numerical examples regarding the fatigue crack growth in a railway axle under two different random loading conditions. The first one applies the probabilistic formulations under a random bending loading. The second one takes into account a service loading spectrum acting on the axle subjected to a random load. Then, the statistical moments calculated are checked by comparison with Monte Carlo simulations. The probabilistic model developed enables an efficient estimation of statistical moments, providing accurate probabilistic results that can be used in design stages, reliability studies or in damage tolerance assessment.

1. Introduction It is currently recognized that the statistical nature of the mechanical properties of solid materials, the scattering of loads and uncertainties inherent to geometrical parameters, affect their structural behaviour and, therefore, probabilistic numerical approaches for reliability assessment and durability prediction are needed. Probabilistic analyses are a useful complement to deterministic analyses that can be used in the design of components or structures and to support the decision-making process of defining maintenance inspections. Those applications still deserve some developments via the extension of the prevailing numerical methodologies. One such area is damage tolerance in the railway sector, where no systematic and detailed probabilistic analysis is provided so far, especially to define in-service inspection intervals in train axles for crack detection [1,2]. Damage tolerance analyses use fracture mechanics to describe the fatigue crack growth process and combined with the definition of a periodic inspection plan, provide a certain level of safety. That means that fatigue cracks are allowed to appear as long as they can be detected with a certain probability before they reach their critical size. Fatigue crack growth is stochastic in nature. The randomness in components under fatigue conditions is inherent to at least three ⁎

aspects. Firstly, it is clearly justified with the experimental tests, where the measurements quite naturally exhibit statistical variations [3–6]. Secondly, the structural components used in industrial sectors as the automotive, aeronautics, transport or railway, are often under variable loading and a large number of cycles during their operation which means scattering of loads [7–9]. The third reason is the fact that geometrical dimensions may have some variations that could significantly affect the functioning of a component [10]. These aspects determine the complexity and the way to address any calculation, estimation or experimental analysis related to the fatigue process. However, as mentioned above, fatigue life prediction is a fundamental aspect for the design and the maintenance planning of components, as is the case of train axles. Consequently, a probabilistic approach that considers input statistical distributions and provides an output response distribution will be more useful than a deterministic one. Addressing the problem from a probabilistic point of view that is also efficient and precise is, therefore, a challenge. In some applications it is enough to obtain certain statistics of the response distribution, for instance, mathematical moments as the expected value and the variance. Over recent decades, probabilistic fatigue crack growth methodologies are becoming more popular in the literature. Some approaches consist in acquiring S–N diagrams for specimens taking into account the

Corresponding author. E-mail address: [email protected] (C. Mallor).

https://doi.org/10.1016/j.ijfatigue.2019.105454 Received 25 July 2019; Received in revised form 17 December 2019; Accepted 21 December 2019 Available online 23 December 2019 0142-1123/ © 2019 Elsevier Ltd. All rights reserved.

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

Nomenclature

p

c C Cthm Cthp d da/ dN f

crack depth/normal semiaxis final crack depth initial crack depth El-Haddad’s parameter crack deepest point crack width at the surface crack surface point superscript above referring to the bending loading case superscript above referring to the bending plus interference loading case crack tangential semiaxis parameter of the crack growth equation in the Paris region parameter for the Kth0 – R relationship parameter for the Kth0 – R relationship random variables number crack growth rate Newman’s crack opening function

g, j

first partial derivative of g with respect to Xj =

g, jk

second partial derivative of g with respect to Xj and

a afin aini a0 A b B B B+I

P q R Smax / x X Y a b K K eff Kth Kth0

µjk µjkl µjklm µY

( )

(

Xk =

2g Xj Xk

g Xj

)

Y 2 Y

evaluation of g (X ) at P vector step increment superscript above referring to the interference loading case j, k, l, m index from 1 to d random variables stress intensity factor (general) K Kc critical stress intensity factor for static unstable crack growth Kmax ,Kmin maximum and minimum K m (x , a) weight functions M bending moment n exponent of the crack growth equation in the Paris region ns steps number N number of applied cycles

gµ i I

0

parameter describing the sigmoidal shape of the crack growth equation in the threshold region mean value vector of X (=(µ X1 , µ X2 , , µ Xd )) parameter describing the sigmoidal shape of the crack growth equation in the toughness region stress ratio (=Kmin/ Kmax ) ratio of the maximum applied stress to the flow stress radial coordinate direction at the axle surface set of random variables (={X1 , X2 , , Xd }) general multivariate function (=g (X )) plane stress/strain constraint factor crack depth increment at point a crack depth increment at point b stress intensity factor range effective stress intensity factor range threshold stress intensity factor range threshold SIF range at R = 0 second central moment (=µ 2 (Xj , Xk )) third central moment (=µ3 (Xj , Xk , Xl )) fourth central moment (=µ4 (Xj , Xk , Xl , Xm )) expected value of Y (=E [Y ]) stress (general) standard deviation of Y (=SD (Y )) variance of Y (=Var (Y ))

Abbreviations FOSM SOSM FEM SIF MC Pr. Eq. Det. NDT POD PDF

variability of fatigue life providing probabilistic fatigue S–N curves for the subsequent component integrity assessment [11,12]. Many other probabilistic models are based on deterministic crack growth equations such as the Paris’ law [13]. For example, several probabilistic studies founded on the Paris law are presented in [14–22], being the model most frequently used. In recent years, due to the importance of the early stages of the crack growth to achieve appropriate life estimation, the first advances considering the complete crack growth curve have been made, especially for its application in components. A modified version of the Paris’ law referred as NASGRO equation [23] represents the state of the art in the field of the fatigue crack propagation problem, being commonly used in railway axles evaluation. Some probabilistic approaches use the Monte Carlo (MC) method on the NASGRO equation to quantify the material uncertainties and loading conditions in the fatigue behaviour of railway axles [24,25]. Additionally, a comparison of the fatigue life calculation in railway axles according to the Paris’ law and to the NASGRO equation is performed in [19]. In this article, for both models, several material parameters levels are considered to quantify how they affect the dispersion of the fatigue life estimation. Early contributions developed general basic probabilistic methods, revealing probabilistic fundamental ideas of the structural reliability problem [26,27]. These probabilistic fundamental ideas were further developed establishing the theory of the first–order second-moment (FOSM) method for reliability analysis [28,29]. This method has been strongly used for many years for those engaged in the probabilistic analysis of structures. Recent research has contributed with new

first order second moment second order second moment finite element method stress intensity factor Monte Carlo probabilistic equation deterministic non-destructive testing probability of detection (of a crack) by NDT probability density function

perspectives on the analysis of components in practical applications, taking advantage of the FOSM method. It is worth highlighting the works made to the field of fatigue crack growth life estimation according to the Paris' law [15,16], in the fatigue crack nucleation stage based on the Coffin-Manson damage model [30] and to evaluate components probabilistically subjected to multiaxial stress states by using the multiaxial fatigue virtual strain energy damage model of Liu [31,32]. However, few attempts to apply the FOSM method considering the NASGRO model have been made. A statistical description of the NASGRO model applying the FOSM method is derived in [33]. This study indicates that the first-order expansions provided are not sufficient to correctly describe the variability of the NASGRO law. The authors discard the use of the second or higher-order expansions alleging that it could produce very complicated formulation becoming the method not viable for practical applications. For that reason, in that article, some simplifications regarding the relationships among the random variables in terms of probability are made to overcome the difficulty. However, it remains unclear whether the FOSM or a higherorder extended version, has the potential to be applied in a successful way. It would thus be of interest to check the validity of the second or higher-order expansions to describe the probabilistic fatigue crack growth based on the NASGRO equation. As shown in the previous literature, the research community has studied the fatigue crack growth process and the stochastic lifetime prediction for many years. Moreover, in engineering practical applications, in-service loading is complex and the experimental information is 2

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

quite limited, so that to solve probabilistic crack growth problems, hypotheses are generally assumed that may imply important limitations and restrictions. Several papers take into account the Paris’ Law or more advanced models as the NASGRO equation. Most of the works use statistical methods that sometimes, as in the MC method, its application implies a high computational cost and the methodology becomes not suitable for practical applications where the computational efficiency is a crucial issue. However, an overall treatment of the fatigue crack growth based on the NASGRO equation combined with the efficient FOSM or higher-order method is missing. Given these premises, the objective of this work is to provide a new probabilistic formulation of the fatigue crack growth phenomenon based on the NASGRO equation. For this purpose, a generalized extension of the FOSM to a full second-order method is combined with the NASGRO model. The methodology derived enables an efficient estimation of the first two statistical moments, the expected value and variance, of the fatigue crack growth lifetime. Two comprehensive application examples of the proposed probabilistic methodology specifically oriented on the fatigue crack growth in a metallic railway axle are presented and discussed. The paper has been organized as follows. Section 2 presents the methods used for the probabilistic analysis. Firstly, the Monte Carlo method is briefly summarized; then the full second-order approach for the moments of functions of random variables is described in detail. In Section 3, an overview of the deterministic NASGRO model is provided. In addition, the probabilistic expressions regarding its integration with the full second-order approach are included. Section 4 illustrates and discusses the stochastic approach proposed by means of two application examples in a metallic railway axle. The results of the full second-order approach are compared with those obtained by the Monte Carlo method with the purpose of checking the performance of the developed methodology in terms of expected values and variances. Section 5 gives an outlook on extensions of the full second-order approach, suggesting future developments. Finally, Section 6 concludes the main novelties introduced by this paper. An Appendix is also provided at the end of the paper to enclose the equations used to obtain the expected value and the variance of the fatigue life according to the NASGRO equation.

Fig. 1. FOSM on a function with one random input variable.

deterministic model obtaining many deterministic results that interpreted as a whole, form a probability distribution. Afterwards the response can be analyzed in the desired statistical terms, for example, expected values and variances. 2.2. Full second-order approach for the moments of functions of random variables The FOSM method is a probabilistic method to determine stochastic moments of a function with random input variables by using Taylor expansions, provided that the function is sufficiently differentiable and that the moments of the input variables are known. The objective is to determine the effect of the input variables randomness on the randomness of the function based on them. As scientific articles use different notations for vectors, random variables or derivatives, the convention and the associated definitions adopted in this paper have been listed in the nomenclature section. Fig. 1 illustrates the simplest version of the FOSM, a function Y = g (X ) with only one random input variable X assumed normally distributed and the first-order Taylor series approximation of g (X ), that is, a linear equation, used to map the input randomness onto the y-axis. The real output shape mapped on y-axis would be in some extent distorted and the distribution would be asymmetric, certainly not normal. When using a first-order Taylor series, the linear mapping provides a normal distribution for Y as represented. To enable the estimation of a non-symmetrical or symmetric non-normal shape accounting the effect of the nonlinear function, a second-order Taylor expansion is proposed. Consider an arbitrary general multivariable function Y = g (X ) of d random variables X = {X1 , X2 , , Xd } . Y is randomly distributed as X is random. To relate the moments of the output variable to the moments of the random input variables, the function Y is rewritten using the multivariate second-order Taylor series expansion about the means vector of the random input variables P = (µ X1 , µ X2 , , µ Xd ) . Integration of the Taylor series, needed for moments estimation, can be performed term by term hence it is straightforward. The generalized second-order Taylor approximation of Y becomes Eq. (1), where gµ is the evaluation of g (X ) at the P vector, g, j is the first partial derivative of g (X ) with respect to Xj and g, jk is the second partial derivative of g (X ) with respect to Xj and Xk .

2. Probabilistic analysis methods In general, it may not be feasible to calculate the expected value and the variance of the response by means of the direct use of the expectation operator on the function which relates the random input and output variables if that function is not relatively simple. In such cases, approximation techniques as the Monte Carlo (MC) and the first-order second-moment (FOSM) are resorted to address the stochastic problem. Once the latter method has been applied, it is common practice to check its results against the MC. Both methods are approximations in their essence, which one is more suitable, will depend on the requirements of problem considered. Probabilistic structural analysis is conditioned mainly by two characteristics, the efficiency and the accuracy. The FOSM method is a low-cost technique in terms of computation time while the MC is an expensive method. If the differences between the results provided are small enough, both methods are considered to be equally accurate and thus the use of the FOSM is recommended for applications that require fast estimations such as technical tools used in design or as an aid to the decision making process of defining periodic inspections of maintenance.

d

Y

gµ +

g, j (x j j=1

2.1. Monte Carlo method

µ Xj ) +

1 2

d

d

g, jk (x j j =1 k=1

µ Xj )(xk

µ Xk )

(1)

The first raw moment of Y , µ Y , is defined as the expected value E [Y ]. The full second-order approximation for the expected value of a general function is shown in Eq. (2), where µjk indicates second central moment µ 2 (Xj , Xk ) .

The MC method is well-known and used in many fields of engineering for solving problems of random variables providing a probabilistic interpretation [34]. The essential idea consists of generating input variability via repeated random sampling and to evaluate a 3

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

E [Y ] = µ Y

gµ +

d

1 2

intensity factor (SIF) range from the maximum and minimum K , Kmax and Kmin , f is the crack opening function, Kth is the threshold stress intensity factor range, K c is the critical stress intensity factor, K eff is the effective stress intensity factor range and C , n , p , and q are material empirically derived constants. For a more detailed description refer to [23]. For a sufficiently small step increment, the crack propagation rate in Eq. (4) can be approximated for the ratio between finite increments Eq. (13), enabling the use of an iterative scheme suitable for computing Eq. (14), where i is the step increment up to the ns steps number.

d

g, jk µjk

(2)

j=1 k=1

The second central moment of Y is the variance, and is usually denoted as Var (Y ) or Y2 . The full second-order approximation for the variance of a general function is given by Eq. (3), where µjkl and µjklm indicate the third and the fourth central moments, that is, µ3 (Xj , Xk , Xl ) and µ4 (Xj , Xk , Xl , Xm ) respectively.

Var (Y ) d

=

2 Y

gµ 2 +

d

d

j=1 k=1

d

j =1 k=1

g, j g, k µjk +

d

g, jk µjk +

d

1 4

d

d

d

d

j=1 k=1 l=1 m=1

g, jk g, lm µjklm + gµ

da dN

d

j =1 k=1 l=1

g, j g, kl µjkl

µY

2

(3)

The first-order second-moment (FOSM) method derives its name from the fact that it uses a first-order Taylor approximation of the function and uses only the first and second moments of the random input variables to determinate the expected value and the variance of the response. When the second-order Taylor is used, it is named as second-order second-moment (SOSM) method, requiring up to secondorder moments of the input variables for the expected value estimation, but for the variance up to fourth-order moments. The moments of the random input variables required to approximate the first to the fourth order moments are given in the Table 1. In this paper the method is referred as full second-order approach, not naming the order of the moments of the input random variables needed for computation as it depends on the order of the moment to approximate.

a N

(13)

ai

ai + 1

ai

Ni

N i+1

Ni

(14)

The crack propagation is calculated at two different points, the deepest point named A and the crack surface point named B as depicted in Fig. 2. To use the NASGRO equation efficiently, fast SIFs calculations, which are dependent on the instant crack shape, is therefore fundamental. To do that, the weight functions method is used. According to the method, for a one-dimensional variation of stresses acting across the crack plane, the relation between the stress intensity factor and the stress distribution is given by the Eq. (15), where (x ) is the stress distribution in the uncracked component and m (x , a) are the weight functions, which are dependent on the instant crack shape and on the position along the crack front in which the SIFs are calculated.

K=

a

(x )·m (x , a) dx

0

(15)

The NASGRO crack growth model is presented from two viewpoints: the deterministic and the probabilistic one.

To determine the crack evolution, the SIFs are evaluated in two points, which determine the crack shape in the next step of the calculation. Those points are the deepest crack point A and the surface point B of Fig. 2. The weight functions used at those points come from [35].

3.1. Deterministic NASGRO equation

3.2. Probabilistic NASGRO equations

3. Fatigue crack growth model

In stochastic fatigue crack growth analyses, the following probability distributions are achievable: the distribution of the crack size at any given number of load cycles and the distribution of the number of load cycles to reach any given crack size. In damage tolerance analyses of railway axles, the number of load cycles distribution is preferred as it enables the definition of interval inspections verifying crack sizes. Introducing Eq. (5) into Eq. (4), isolating dN and using the discretised version gives Eq. (16).

NASGRO expressions are enclosed in Eqs. (4)–(12).

(

)

Kth p K q K

1 da = C ( K eff ) n dN 1

(

K eff =

1 1

f=

Kc

)

(4)

f ( K) R

Kth = Kth0

(5) a a + a0

(

(1

Cthp Cthm

Cth =

max

)

1 + Cth·R

1 f A0 )(1

R)

if R 0 if R < 0

A0 = (0.825

0.34 + 0.05 2) cos

A1 = (0.415

0.071 )

A2 = 1

A0

Smax 0

(8)

Smax /

0

A3

1

fi

1

Ri

Ki

)

n

1

q

i Kth

p

(16)

Ki

Table 1 Moments of the random input variables required using first order (FO) and second order (SO) approximations.

1/

2

(

i Kmax Kc

The expected value of the fatigue life E [N ], is obtained by applying the expectation operator to the total number of cycles N which is the summation of the ns life increments dN i up to a final crack depth afin . Applying the linearity property of the expectation operator, leads to Eq.

(9) (10) (11)

1

A1

C

(7)

max (R, A 0 + A1 R + A2 R2 + A3 R3) if R > 0 A0 + A1 R if R 0

A3 = 2A 0 + A1

dN i =

(6)

1

dai

(12)

where da/ dN is the crack propagation rate, N is the number of applied cycles, a is the crack depth, R is the stress ratio, K is the stress 4

Moment to approximate

Moment required

Ordinal

Significance

FO

SO

1 2 3 4

Expected value (Raw) Variance (Central) Skewness (Central Normalised) Kurtosis (Central Normalised)

1 1, 2 1, 3 1, 4

1, 1, 1, 1,

2 2, 3, 4 2, 3, 4, 5, 6 2, 3, 4, 5, 6, 7, 8

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

relatively long cracks larger than 2 mm in approximately 95% cases. That is, a crack of 2 mm length is not detected in 5% cases. Therefore, there is a risk that an existing crack is not detected, and so the existence of a crack should be considered. In the cases of study, a semicircular initial crack aini of 2 mm was postulated at the T-transition, as indicated with a red line in Fig. 3, in which is also defined the radial coordinate system x at the axle surface. The fatigue crack growth material parameters required for the NASGRO model evaluation are obtained from the bibliography and given in Table 2. 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-fit loading produced by the wheel mounting with interference. The bending moment level selected M , equal to 70.32 MN mm, 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. This assumption implied a worst case scenario since such stress level corresponded to the maximum one for axle bodies according to the EN 13103 standard [43]. The wheel was press-fitted with 0.286 mm interference. The stress distributions normal to the crack surface needed for the SIF calculation were obtained via the finite element method (FEM). The stresses at the interference I , at the reference bending moment B and at the combination of both B + I are shown in Fig. 4. The railway axle rotates as the train moves, and the bending load always acts in the vertical plane, as a result, every revolution the crack goes through all possible positions with respect to the applied moment and therefore the resulting load type is rotary bending. In consequence, in the next application examples, the total stress distribution on the crack depended on its angle location with respect to the axle rotation as it is shown in Fig. 5. As mentioned before, the validity of the proposed method is verified by the following two numerical examples considering two different random loading conditions.

Fig. 2. Axle cross-section with a postulated crack.

(17). ns

ns

dN i =

E[N ] = E i=1

E[dN i]

(17)

i=1

The fatigue life variance Var (N ) , is obtained by applying the variance definition to N . Using the formula for the variance of the sum of random variables gives Eq. (18). ns

ns

dN i =

Var (N ) = Var i=1

ns 1

ns

Var (dN i1) + 2 i1= 1

Cov (dN i1, dN i2) i1= 1 i2 = i1+ 1

(18) In Eq. (17), is obtained by applying the full second-order approximation developed to obtain the stochastic moments of first order Eq. (2), on the dN i function Eq. (16). Similarly, the moments of dN i in Eq. (18), are obtained by applying the full second-order approximation developed to obtain the stochastic second central moment Eq. (3). Eqs. (17) and (18) are referred to as the probabilistic NASGRO equations for the expected value and variance calculation respectively. Summarizing, for every ith crack growth increment the first raw and the second central moments of the fatigue lifetime increment dN i are calculated by applying the full second-order approximation method on the discretised version of the NASGRO equation, requiring the first and second partial derivatives of the NASGRO equation with respect to the random input variables, and the first to fourth order moments of the random input variables.

E[dN i]

4.1. Example 1: Railway axle under random bending moment The purpose of this case is to elucidate the strategy proposed and make clear that it can well predict the expected value and variance of the fatigue lifetime based on NASGRO model. In this illustrative application example, the bending moment in a railway axle was assumed as a random input variable normally distributed with a standard deviation equal to the 5% of the mean value. The bending stress distribution shown in Fig. 4 was considered as the reference bending stress amplitude for the mean value of bending moment B . The parameters of the distribution are shown in Table 3. As normally distributed, the third central moment was zero and the fourth central moment was equal to 3 M4 . It is worth noting that the bending moment distribution used in this example to examine and demonstrate the validity of the proposed method corresponded to an extreme situation that was not

4. Numerical results and discussion In this section, two examples are given to verify the validity of the proposed full second-order method. The first example is intended to clearly illustrate the stochastic approach proposed in a railway axle under random bending moment. However, under realistic service conditions, the axles of rail vehicles do not experience long periods of uniform amplitude loading, but they are instead subjected to a stochastic stress spectrum. In order to further verify the viability of the proposed method, a second example takes into account the loading spectrum acting on a railway axle over its service, as well as the randomness of the cargo into the vehicle. Prior to those examples, several characteristics that are common to both cases are described, regarding the railway axle geometry, the initial crack considered, the material properties and the reference bending and interference stresses. The probabilistic analysis methods presented were applied to the fatigue crack growth of the railway axle shown in Fig. 3, under different random loading conditions that are presented hereafter in two examples. The axle body was 173 mm in diameter and it was made of EA1N steel defined in the EN 13261 standard [36]. The probability of detection (POD) of a crack depends not only on the non-destructive inspection (NDI) technique but also on the actual crack size [37]. For instance, magnetic particle inspection can detect

Fig. 3. General axle view of a non-powered wheelset. 5

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

Table 2 NASGRO fatigue crack growth material parameters. Parameter C

a

Value

Reference

Cthp

3.3197E−10 MPa 2.09 1.3 0.001 2434.9 MPa mm 233.7 MPa mm 1.442

Smax / a0

2.5 0.3 0.0381 mm

n p q Kc Kth0

mm1

[38]

n/2

[38] [38] [38] [39] [38] [38]

−0.02

m Cth

a

n

0

[38]

[40,41] [41] [41,42]

Exponents in the C parameter units make reference to n .

Fig. 6. Crack depth evolution versus the number of cycles and histograms, frequency normalized, for different crack depths.

scheme described in Section 3.1, increasing the crack depth and keeping a semielliptical shape while growing up to a final crack depth afin of 50 mm. The final length considered was shorter than the critical length determined from fracture toughness. The definition of the railway axle geometry, Fig. 3; the initial crack of 2 mm; the material properties, Table 2 and the reference bending and interference stresses, Fig. 5; are included at the beginning of the Section 4. With the aim of checking the statistical moments computed via the full second-order approach, 10 000 MC simulations were performed providing the results shown in Fig. 6. The lower plot shows the different evolutions of the crack depth versus the number of cycles for each realization and it also represents by a red line the computed mean value. The upper graph shows six histograms, frequency normalized, of the random output variable fatigue life N for 5, 10, 20, 30, 40 and 50 mm crack depth levels. The probability distribution of the fatigue life N can be analyzed for every crack depth a into the 2–50 mm interval. It is observed from the histograms that as the crack depth increases, the width of the histogram, that is, the variability of N , increases too. The underlying probability distribution of N differs from a normal distribution according to a normality test based on [44]. To give insight into the non-normally distributed fatigue life N , a normal probability plot is shown in Fig. 7. The deviations from the straight line indicate departures from normality and they evidence the typical inverted C shape of a right-skewed distribution. The full second-order approach presented was applied to calculate the first moment and the second central moments of dN i at every crack i i depth, Eq. (16), with the two input random variables Kmax and Kmin . Then, the expected value µN and the variance N2 of the fatigue life N were obtained from the ith moments, providing a continuous result along the crack depth a . See the Appendix A for the complete mathematical formulation of the expected value and the variance of the discretised NASGRO equation derived for the application example. To check the accuracy of the method in terms of central tendency, the µN histories provided by the Monte Carlo (MC) and by the probabilistic equation (Pr. Eq.) are compared in Fig. 8. It is also included the deterministic (Det.) number of cycles for the mean value of the bending moment µM . For the crack depth levels used in the MC results, the values of the previous curves are listed in Table 4. The expected value computed from the MC simulations distribution is considered as the framework of reference for comparison. The error in the Pr. Eq. is about 0.3%, staying constant as the crack depth increases. The maximum difference of the Pr. Eq. leads to an expected value about 1000 cycles higher. The percentage difference of the Det. number of cycles is about the 2% and the maximum difference is less than 6000 cycles, being the Det. value lower. In the latter case, the

Fig. 4. Stress distributions at the T-transition in axial direction at the press-fit with interference, at the bending and at the superposition of both load cases.

Fig. 5. Illustration of the bending and interference stresses at the axle cross section in two different angle rotations. Table 3 Normal distribution parameters of the random variable bending moment. Probability Distr.

µM [MN mm]

Normal

70.32

2 M

[MN mm]2

12.37

representative of the actual load spectrum acting on a railway axle over its service life. From Fig. 5, it may be inferred that the maximum and minimum stresses were obtained at = 0 and = . The maximum and minimum stresses were related to each other as a result of the stresses symmetry. 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. Accordingly, the SIF was also stochastic, and given the axle cyclic rotary character, it was described as two random variables Kmax and Kmin that were related to each other. The two random variables were correlated variables so they were not independent. It is noteworthy that Kmax and Kmin are explicitly involved in the NASGRO equation. The combination of the loading conditions was applied repeatedly according to the iterative 6

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

variable, so it is more meaningful to interpret. For the crack depth levels used in the MC results, the values of the previous N curves are listed in Table 5. Once again, the standard deviation provided by the MC simulations is considered as the framework of reference for comparison. It is trivial that the deterministic calculation does not account the variability of the response, and then it does not provide standard deviation results. It is observed that the larger the crack, the more the variability in cycles. This behavior has previously been inferred from the histograms shown in Fig. 6. The Pr. Eq. also reproduces this tendency. The error in the probabilistic equation ranges from 1% to 3%, increasing as the crack depth increases. The maximum difference of the Pr. Eq. calculates a standard deviation about 1400 cycles higher. At this point, the problem of constructing a probability density function (PDF) from prescribed moments arises. It is worth to recall that the current Pr. Eq. formulations provide only the expected value and the variance of the underlying lifespan probability distribution. Therefore, it is possible to infer probability distributions of two parameters. To check the viability of the strategy presented in terms of life distributions, two scenarios were considered: (i) the lifespan was assumed to be normally distributed; (ii) the lifespan was assumed to be log-normally distributed. Despite the fact that the underlying probability distribution of N is not normally distributed as shown in Fig. 7, the normal distribution is considered since it is the most commonly assumed distribution when there is not more information available. On the other hand, the log-normal distribution is more versatile than the normal because it has a range of shapes to represent asymmetric data and therefore is generally accepted for cycles-to-failure in fatigue reliability analysis. Both probability density functions were inferred from the expected value and the variance provided by the Pr. Eq. formulations. The parameters of the normal distribution are directly the expected value and the variance, whereas the scale and shape parameters of the log-normal are functions of them. The location parameter was set equal to zero, resulting in the two-parameter log-normal distribution. The physical meaning of the location parameter is realistic, since it indicates that crack growth occurs after any given cycle. The normal distribution, the log-normal distribution and the MC histogram of the fatigue life N for a crack depth equal to 50 mm are compared in Fig. 10. The superiority of the log-normal over the normal distribution to represent the MC data is clear from the behavior of the tails and the peak. It is better in describing the lower tail of the lifespan distribution where the skewed data causes the symmetrical normal distribution to slightly underestimate the lifespan. The overall similitude between the log-normal and the MC histogram confirms the good selection of the distribution.

Fig. 7. Normal probability plot of the number of cycles N provided by the Monte Carlo (MC) for crack depth of 50 mm and the least squares fit line of the data.

Fig. 8. Expected value of N history values provided by the Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr.Eq.) and the deterministic (Det.) number of cycles. Table 4 Expected value of N provided by Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr. Eq.) and fatigue life according to the deterministic calculation (Det.). a

MC µN

Pr. Eq. µN

Det. N

Pr. Eq.-MC Error

[mm]

[Cycles]

[Cycles]

[Cycles]

[%]

Det.-MC Diff . [%]

5 10 20 30 40 50

158 230 301 344 375 397

158 231 301 345 376 398

155 226 296 339 370 392

0.32% 0.27% 0.25% 0.24% 0.24% 0.24%

−1.96% −1.67% −1.49% −1.41% −1.36% −1.32%

249 755 204 724 677 901

749 384 945 544 566 853

154 906 724 879 587 650

deterministic N is lower than the MC expected value for every crack depth. It can be explained due to the existing right-skewed distribution of N observed previously. To check the accuracy in terms of dispersion, the standard deviation history values N provided by the MC and by the Pr. Eq. for the variance estimation are compared in Fig. 9. The standard deviation is reported instead of the variance even though the variance is more convenient when developing the probabilistic formulations, since the standard deviation has the convenience of being expressed in units of the original

Fig. 9. Standard deviation of N history values provided by the Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr.Eq.) 7

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

Table 5 Standard deviation of N provided by Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr. Eq.) for the variance calculation. a [mm]

MC

[Cycles]

[Cycles]

Pr. Eq.-MC Error [%]

5 10 20 30 40 50

26 35 44 48 52 54

27 36 44 49 53 56

1.17% 1.43% 1.64% 1.88% 2.18% 2.54%

N

Pr. Eq. N

750 828 031 947 360 715

062 342 752 869 504 103

Fig. 11. Stress spectrum (mileage 15 000 km) considered in the probabilistic analysis.

acting on a railway axle over its service, as well as the randomness of the load into the vehicle. The load spectrum was derived from the one available in the UIC B 169/RP 36 report [7]. The spectrum was obtained from measurements taken on wheelsets in representative traffics with an axle load of 22.5 tonnes and an axle diameter of 174 mm. The stress reversals collective was discretised into 29 classes or loading blocks, accumulating the number of cycles by the common cycle counting rain-flow method. In this example, the stress spectrum was scaled to reproduce the stresses at the axle T-transition of interest, Fig. 3. To achieve this, the loading block with the highest amplitude was adjusted to the bending stress distribution represented in Fig. 4 and the other blocks were scaled accordingly. The spectrum was also transformed into a 5 million load alternations sequence which corresponded to a mileage of 15 000 km in relation to a wheel diameter of 920 mm. Finally, the loading blocks were sorted in descending sequence, that is, from highest to lowest load amplitudes. The resulting service loading spectrum is shown in Fig. 11, preserving its original outline despite the adaptations aforementioned. The iterative calculation scheme for crack growth takes into account the spectrum via varying the bending stress component and applying the pertinent number of cycles at each block, both described by the load spectrum. This procedure evaluates each loading block in the preestablished sequence, computing the stress intensity factor range K and the threshold stress intensity factor range Kth at the current crack size. If K is above Kth , the crack grows and the number of cycles is accumulated. If K is below Kth , the current spectrum block does not contribute to crack growth, and therefore the crack shape remains unchanged for the next iteration, however, the number of cycles is accumulated. In other words, depending on the stress intensity factors ranges which depend on the current crack size and on the stress level at the current block of the spectrum, the crack grows or not but the number of cycles applied in the block is always accumulated. The uncertainties of the cargo in the vehicle at each block of the spectrum were represented by a random bending stress distribution centered at the different stress amplitude levels of the load spectrum. It means that the stress spectrum shown in Fig. 11 was considered as the reference bending stress amplitude for the mean value of the load at each block. These random variables at every class are related to a common origin, the randomness of the load into the vehicle. In consequence, all the classes were considered as normally distributed perfectly correlated random variables with the same coefficient of variation equal to the 5%. Thus, all the classes have the same ratio of their

Fig. 10. Histogram of fatigue life N provided by the Monte Carlo (MC) and PDFs of the normal and the log-normal distributions inferred from the probabilistic NASGRO equations (Pr.Eq.) statistics for 50 mm crack depth.

Taken the observations from the results together, the next outcomes are obtained:

• The expected values and the standard deviation or variance pro• • • •

vided by the MC and by the probabilistic NASGRO equations are really close, giving errors between them small enough to consider both methods to be equally accurate. The expected value and variance provided by the Pr. Eq. are useful to infer the PDF of the crack growth distribution under stochastic conditions. The output distribution of N shown in the example is not normally distributed, even with normal inputs. It is closer to the log-normal distribution. The key advantage of the proposed full second-order method is the lower computational time, similar to that employed for a deterministic one. So, this method is considered to be advantageous compared to the conventional MC method. The constant amplitude bending moment considered in this example, and so the stresses do not reproduce the actual load state of a railway axle. Therefore, it should be noted that the results obtained in terms of number of cycles are not representative of the lifespan of the component in service.

4.2. Example 2: Railway axle under random load into the vehicle and stress spectrum During train operation the railway axles are subjected to cyclic loading with variable amplitude. The variability is caused by several regimes of train operation related to: the weight level (ranging from unloaded to fully laden train), the velocity of the vehicle, braking efforts or to the track layout (e.g. the train goes through a straight or curved rail, over crossovers, switches, rail joints, etc.). For a proper life assessment of components subjected to variable amplitude fatigue loads, the shape of the stress spectrum and the stochasticity of the inputs are key factors. The present example considers the load spectrum 8

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

standard deviation to its mean. The assumptions taken go one step further by considering a complex scenario where the concepts of stress spectrum and a random load are combined and encompassed in the cycle-by-cycle integration iterative scheme. The stress spectrum represents the amplitude variation of the bending moment due to dynamic effects such as the velocity of the vehicle or by the track layout as the train goes through straight, curved tracks, crossovers, etc. On the other hand, the random load represents the uncertainties related to the amount of cargo in the vehicle, which originates the bending moment. The purpose of the assumptions taken is to generalize the consideration of the load at each stress level, providing the load with a probabilistic character and, at the same time, representing it in a more realistic way. The assumptions essentially corresponded to the actual loading conditions acting on a railway axle over its service life. The random load and stress spectrum combination was applied repeatedly in a loop, increasing the crack depth progressively up to the final crack depth afin equal to 10 mm. The final crack depth was shorter than the critical length determined from fracture toughness. Nevertheless, due to the high crack propagation rate in the last stage of fatigue crack propagation, a longer final crack depth does not lead to significantly higher number of cycles. The stress spectrum applied every loop, was every time sorted in the same fashion as presented in Fig. 11, i.e. in descending sequence. The definition of the railway axle geometry, Fig. 3; the initial crack of 2 mm; the material properties, Table 2 and the reference bending and interference stresses, Fig. 5; are included at the beginning of the section 4. The full second-order approach presented was applied to calculate the first moment and the second central moments of dN i at every crack i i depth, Eq. (16), with the two input random variables Kmax and Kmin , accumulating the number of cycles according to the stress spectrum. Then, the expected value µN of the fatigue life N was obtained from the ith moments, providing a continuous result along the crack depth a. See the Appendix A for the complete mathematical formulation of the expected value and the variance of the discretised NASGRO equation derived for the example. The results provided by the proposed methodology were compared with the results of a series of 10 000 MC simulations. To check the accuracy of the method in terms of central tendency, the µN histories provided by the Monte Carlo (MC) and by the probabilistic equation (Pr. Eq.) are compared in Fig. 12. It is also included the deterministic (Det.) number of cycles for the case of considering the mean value of the load into the vehicle at each block of the spectrum. First of all, the staircase shape in the Pr. Eq and in the Det. curves was caused by blocks of cycles in the spectrum with load amplitudes that did not contribute to the crack growth as the K was below Kth , at the pertinent crack size. Hence, the crack depth a remained unchanged until the next damaging loading block, and the number of cycles in the block was accumulated. Nevertheless, it is important to remark that not damaging load levels in the stage of a short crack could become damaging with a growing crack in its stage of a long crack. The staircase shape phenomenon was not observed in the MC curve given that it was the computed mean value of the MC simulation series. Moreover, the fatigue crack growth curves are usually presented in a logarithmic x-axis, where the staircase effect is not noticed. The above curves are shown using natural scale to recognize better the differences between the results of the three curves. The values of the previous curves for 3–10 mm crack depths are listed in Table 6. The expected value computed from the MC simulations distribution is considered as reference for comparison. The error in the Pr. Eq. is less than 4% at every crack depth, fluctuating as a result of the staircase shape. The percentage difference of the Det. number of cycles ranges from 24 to 28%. In the MC and Pr. Eq. calculations, the residual lifetime from the initial crack depth of 2 mm to 10 mm is more than 100 million cycles corresponding to 314 and 310 thousand kilometers respectively, in relation to a wheel diameter of 920 mm. On the other hand, the Det. calculation estimates around 82 million cycles which leads to 240 thousand kilometers. The expected value of the lifetime provided by the

Pr. Eq. is about 4000 km lower than the MC, considered negligible, while the lifetime provided by the Det. is 74 000 km lower. The estimation of the lifespan dispersion taking into account the stress spectrum combined with the randomness of the load is challenging. This is because each block of the spectrum contributes to crack growth or not depending on the level of the random load amplitude and on the pertinent crack shape. Additionally, not damaging load levels in the stage of a short crack could become damaging in its stage of a long crack. Both circumstances induce additional variability that requires further developments. Overcoming this issue will make the construction of the probability density function of the underlying lifespan distribution possible, based on the expected value and the variance provided by the Pr. Eq. A damage-tolerance-based inspection interval definition consists in two major elements: facture-mechanics-based residual lifetime, i.e. crack depth a vs. cycles N , and the probability of detecting (POD) a crack of a certain size by non-destructive testing (NDT). The calculation of inspection periodicity is beyond the scope of the present research, however, in a simple manner, the longer the expected lifetime the longer the feasible inspection intervals. For a detailed description of the definition of inspection intervals in railway axles see [1]. Usually, the deterministic lifetime calculation is considered for computing the inspection intervals. As previously observed, the Det. evaluation could lead to an overly conservative estimation, and therefore it would lead to short inspection intervals with excessive maintenance costs and a reduced performance life due to premature axle replacement. Alternatively, the expected value of the lifetime provided by the Pr. Eq. can be considered for the calculation of the inspection intervals. As a result, it could enable an optimization, in other words, an extension of the inspection intervals of railway axles that could reduce the total life-cycle cost of wheelsets. The choice between the previous two residual lifetime curves, the Det. and the expected value provided by Pr. Eq., implies a different degree of conservativeness. The probability of failure (crack depth a reaches the final crack depth afin ) according to the probability density function of the underlying lifespan distribution can be defined as a metric designed to quantify the degree of conservativeness of any residual lifespan curve. Specifically, the probability of failure associated to any residual lifetime would be provided. Gathering the observations from the results, the next outcomes are achieved:

• The expected values of the lifespan provided by the MC and by the •

probabilistic NASGRO equations are really close, even considering the load spectrum acting on a train axle over its service, as well as the randomness of the load into the vehicle. From an engineering perspective, the accuracy of the proposed approach is fully acceptable. Damage tolerance assessment based on a probabilistic fatigue crack

Fig. 12. Expected value of N history values under the load spectrum provided by the Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr.Eq.) and the deterministic (Det.) number of cycles applying the load spectrum. 9

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

Table 6 Expected value of N under the load spectrum provided by Monte Carlo (MC) and by the probabilistic NASGRO equation (Pr. Eq.) and fatigue life according to the deterministic calculation (Det.) applying the load spectrum. a

MC µN

Pr. Eq. µN

Det. N

Pr. Eq.-MC Error

[mm]

[Cycles]

[Cycles]

[Cycles]

[%]

Det.-MC Diff . [%]

3 4 5 6 7 8 9 10

99 172 651 104 906 936 106 753 685 107 561 536 108 018 632 108 358 866 108 618 244 108 839 185

95 598 667 101 339 372 107 086 970 107 157 160 107 242 613 107 337 457 107 433 821 107 533 702

71 77 81 81 82 82 82 82

−3.60% −3.40% 0.31% −0.38% −0.72% −0.94% −1.09% −1.20%

−27.62% −26.60% −23.26% −23.77% −24.02% −24.17% −24.27% −24.33%

• •

776 006 920 989 073 166 261 359

784 829 891 815 732 911 648 920

•

growth estimation, can optimize the definition of inspection intervals of railway axles that could reduce the total life cycle cost of wheelsets. The probabilistic formulations provided can lead to a better description of the full complexity of fatigue crack propagation processes than the deterministic one. The probabilistic approach presented aims to reflect the railway axle lifetimes under service conditions. Further developments of the accumulation of lifespan variance taking into account a stress spectrum combined with the randomness of the load are required. More precisely, this is due to the fact that depending on the level of the random load amplitude and on the pertinent crack shape, the block contributes to crack elongation or not. Also, not damaging load levels in the stage of a short crack could become damaging in its stage of a long crack. Both circumstances induce additional output variability that requires further treatment. As a result, an overall accumulation of the variance considering these effects in the crack growth predictions remains as an open point for future developments.

•

5. Outlook This section summarizes some proposals for future work that would allow the method presented to be extended, contributing to a better knowledge of the distribution of fatigue life.

• Consideration of more input variables of interest as random due to

the correlations between random parameters could lead to incorrect results. In consequence, taking account of the correlations between parameters is of major importance for a proper probabilistic analysis. The strategy to include the variability of the aforementioned parameters into the architecture of the method is the same as the i one applied in the examples for two input random variables Kmax i and Kmin . First, the full second-order approach is applied to calculate the first and the second central moments of dN i , Eq. (16), at every crack depth, with the four input random variables Kthi0, C i, ni and K Ci . These moments are obtained by applying the equations developed to obtain the stochastic moments of first order, Eq. (2), and the stochastic second central moment, Eq. (3). Then, the expected value µN and the variance N2 of the fatigue life N are obtained from the previous ith moments, applying the probabilistic NASGRO equations for the expected value and variance calculation, Eqs. (17) and (18) respectively, providing a continuous result along the crack depth a . Prediction of more statistics of the underlying lifespan probability distribution. A key point to construct a probability density function (PDF) from prescribed moments is the selection of the probability distribution to adjust. Alternative statistical distributions can be adopted and their suitability evaluated. The choice does influence the goodness of the strategy presented. It can also consistently determine a conservative or non-conservative estimation. In some circumstances, depending on the uncertainties involved, a twoparameter distribution may not be enough. Predicting higher order moments of the underlying distribution would enable the adjustment of probability distributions with more than two parameters if needed. For instance, statistical moments such as the third and the fourth central moments that can be related to the shape of the distribution could be predicted for that purpose. Further development of the accumulation of lifespan variance considering the load spectrum variability contribution. As a result of the spectrum combined with the randomness of the load, different amount of blocks are eventually damaging, i.e., contribute to the crack growth, and it generates additional variability. Although the full second-order approach is able to properly calculate the variance contribution related to the randomness of the load into vehicle as verified in the example 1, the aforementioned additional variability requires further treatment. This improvement would give information in terms of dispersion about the predicted lifespan in railway axles under service conditions.

6. Conclusions

the inherent randomness of the variables involved in the fatigue crack growth process. The approach presented provides the uncertainties of the lifespan in terms of two statistical moments according to the NASGRO fatigue crack growth law, via input variables modelled as random variables. In other words, a stochastic fatigue crack growth model is provided, with respect to whatever parameters selected. Therefore, the applicability of the full secondorder approach goes beyond the illustrated examples that only consider the variability in the loading conditions. As an example, the scatter of the fatigue crack growth curve can be modelled by means of considering the material parameters involved in the NASGRO equation as random variables. The consideration of them as random is straightforward in the architecture of the methodology presented. The material properties such as the threshold stress intensity factor range at R = 0 , Kth0 , the intercept C and slope n in the Paris region [41], and the critical stress intensity factor K c seem to have experimental scatter [24,40,45,46]. The statistical variability of the parameters among specimens can be characterized via their statistical moments as the empirical mean value, the empirical variance, etc. In addition, these parameters can be tested for possible correlations. It is important to remark that the method proposed takes into account the relationships among the random parameters in terms of probability in addition to the statistical moments of the parameters individually. Notice that disregarding

This work presents a full second-order probabilistic formulation to predict the statistical moments, expected value and variance, of the fatigue crack growth lifetime based on NASGRO model using information about the input random variables distribution. Those moments are helpful to describe the crack growth phenomenon and infer its probability density function, under stochastic conditions such as under a random bending moment loading or under a real service load spectrum acting on a component combined with a random load. The accuracy of the fatigue lifetime moments obtained via the full second-order approach presented, and its efficiency when compared to an equivalent Monte Carlo method analysis, prove the good performance of the proposed approach. Therefore, it is remarkable that the method can be used for predicting the stochastic characteristics of crack length history in applications that require a low computation time as in reliability studies implemented in design stages. Damage tolerance assessment of components such as railway axles, can benefit from the results provided to define interval inspections with certain level of safety, optimizing maintenance costs. Declaration of Competing Interest The authors declare that they have no known competing financial 10

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

interests or personal relationships that could have appeared to influence the work reported in this paper.

Industry and Competitive through the National Programme for Research Aimed at the Challenges of Society that financially supported the project RTC-2016-4813-4.

Acknowledgement The authors acknowledge the Spanish Ministry of Economy,

Appendix A. Full second-order approximation for the expected value and for the variance of the discretised NASGRO equation used in the examples In this Appendix, the equations that enable the calculation of the expected value and the variance of the discretised NASGRO equation used in the i i examples are presented. Starting from the dN i function Eq. (16), with two random input variables Kmax and Kmin and rewriting it by means of the i i ) leads to Eq. (A.1). , µ Kmin multivariate Taylor series up to second-order about the mean value vector at every ith increment P i = (µ Kmax

dN i C

(

1

dai 1

fi

1

Ri

i (Kmax

Ki

)

n

i Kmax Kc

q

i Kth

p

1

i µ Kmax )2 +

dN Kmax

+ µ i K

Ki

µ i K

max µK i min

max µK i min

1 2dN 2 Kmin2

i (Kmin

µ i K

max µK i min

i (Kmax

i µ Kmin )2 +

dN Kmin

i ) + µ Kmax

2dN 1 2 2 Kmax Kmin

µ i K

max µK i min

i (Kmin

max

1 2dN 2 Kmax 2

µ i K

max

µK i

min

i (Kmax

µ i K

i ) + µ Kmin

i i µ Kmax )(Kmin

i ) µ Kmin

(A.1)

µK i

min

The expected value of , is obtained by applying the equation developed to obtain the stochastic moments of first order Eq. (2). In the case study, considering the approximated dN i function Eq. (A.1), gives Eq. (A.2).

dN i

E [dN i] = µdN i

dai C

1 1

fi Ri

E[dN i],

1

i Kmax Kc

q

p

1

i Kth Ki

n

Ki

+ µ i Kmax µ i Kmin

2dN Kmax 2

1 2

µ i Kmax µ i Kmin

i Var (Kmax )+

2dN Kmin 2

µ i Kmax µ i Kmin

i Var (Kmin )+2

2dN Kmax Kmin

µ i Kmax µ i Kmin

i i Cov (Kmax , Kmin )

(A.2)

Likewise, the variance of dN i , Var (dN i ) , is calculated by applying the equation developed to obtain the stochastic second central moment Eq. (3). The full second-order approximation of the variance in the application example is shown in Eq. (A.3). 2

Var (dN i) =

2 dN i

dai C

1 1

fi Ki Ri

1

i Kmax Kc

q

i Kth

p

n

1

2

+

2

dN µ i Kmax Kmax

i Var (Kmax )+2

dN µ i Cov(Kimax , Kimin) + Kmin Kmax µ i Kmin

µ i Kmin

µ i Kmin

µ i Kmax µ i Kmin

Ki

dN µ i Kmax Kmax

dN µ i Kmin Kmax

Var

µ i Kmin

2

1 (Kimin) + 4

2dN

Kmax 2

µ i Kmax µ i Kmin

i , Ki , Ki , Ki ) + 4 (Kmax max min min

i i i i , Kmax , Kmax , Kmax )+4 µ4 (Kmax

2dN µ i Kmax 2 Kmax

2dN i i i i ) +2 µ i , Kmax , Kmax , Kmin µ4 (Kmax Kmax Kmin Kmax

µ i Kmin

2dN µ i Kmax Kmin Kmax µ i Kmin

µ i Kmin

2dN i , Ki , Ki , Ki ) + 4 µ i µ4 (Kmax min max min Kmax Kmin Kmax

2 i i , Ki , Ki ) + (Kmax , Kmin min min

2

2dN

Kmin2

µ i Kmax µ i Kmin

2dN i , Ki ) + µ i Cov (Kmax min Kmax Kmin Kmax µ i Kmin

i , Ki , Ki ) + (Kmax max min

dN µ i Kmax Kmax µ i Kmin

i , Ki , Ki ) + 2 (Kmin max max

dN µ i Kmin Kmax µ i Kmin

2dN µ i Kmax Kmin Kmax

µ i Kmin

dai

i , Ki , Ki , Ki ) + µ4 (Kmin min min min

C

2dN i ) + µ i Var (Kmin Kmin2 Kmax µ i Kmin

1

fi Ki Ri

p

1

i Kth Ki

n

dN µ i Kmax Kmax

2dN i , Ki , Ki ) + µ i µ3 (Kmax min min Kmin 2 Kmax µ i Kmin

1

i Kmax Kc

q

1

µ i Kmin

dN µ i Kmin Kmax µ i Kmin

2dN i , Ki , Ki ) + µ i µ3 (Kmin min max Kmax Kmin Kmax µ i Kmin

11

µ i Kmin

µ i Kmax µ i Kmin

µ i Kmin

2dN µ i µ4 Kmin2 Kmax µ i Kmin

2dN µ i µ4 Kmin2 Kmax µ i Kmin

2dN i µ i Var (Kmax )+ Kmax 2 Kmax µ i Kmin

2dN i , Ki , Ki ) + 2 µ i µ3 (Kmax max max Kmax 2 Kmax µ i Kmin

2dN µ i Kmax 2 Kmax

dN µ i Kmax Kmax µ i Kmin

2dN µ i µ3 Kmax Kmin Kmax µ i Kmin

2dN µ i µ3 Kmax 2 Kmax µ i Kmin

dN µ i Kmin Kmax µ i Kmin

2dN i , Ki , Ki ) µ i µ3 (Kmin min min Kmin2 Kmax µ i Kmin

(E [dN i])2

(A.3)

International Journal of Fatigue 133 (2020) 105454

C. Mallor, et al.

References

[23] Forman RG, Mettu SR. Behavior of surface and corner cracks subjected to tensile and bending loads in Ti-6Al-4V alloy. ASTM STP 1131 Am Soc Test Mater Phila PA 1990:519–46. [24] Beretta S, Carboni M. Experiments and stochastic model for propagation lifetime of railway axles. Eng Fract Mech 2006;73:2627–41. https://doi.org/10.1016/j. engfracmech.2006.04.024. [25] Beretta S, Villa A. A RV approach for the analysis of fatigue crack growth with NASGRO equation. 4th Int ASRANAT Colloq 2008:1–7. [26] Mayer M. Die Sicherheit der Bauwerke und ihre Berechnung nach Grenzkräften anstatt nach zulässigen Spannungen. Berlin: Springer Verlag; 1926. [27] La Mayer M. 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 CA. A probability-based structural code. ACI J Proc 1969;66:974–85https:// doi.org/10.14359/7446. [29] Cornell CA. Structural safety specifications based on second-moment reliability analysis. IABSE symposium on concepts of safety of structures and methods of design; 1969. doi: 10.5169/seals-5948. [30] Núñez JL. 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 JR, Núñez JL. Probabilistic formulation of the multiaxial fatigue damage of Liu. Int J Fatigue 2011;33:460–5. https://doi.org/10. 1016/j.ijfatigue.2010.10.003. [33] Corbetta M, Sbarufatti C, Manes A, Giglio M. On dynamic state-space models for fatigue-induced structural degradation. Int J Fatigue 2014;61:202–19. https://doi. org/10.1016/j.ijfatigue.2013.11.008. [34] Sundararajan CR. Probabilistic structural mechanics handbook: theory and industrial applications. Springer Science & Business Media; 2012. [35] Wang X, Lambert SB. Stress intensity factors for low aspect ratio semi-elliptical surface cracks in finite-thickness plates subjected to nonuniform stresses. Eng Fract Mech 1995;51:517–32. https://doi.org/10.1016/0013-7944(94)00311-5. [36] EN 13261:2009+A1:2010: Railway applications – wheelsets and bogies – axles – product requirements. European Committee for Standardization; 2010. [37] Zerbst U, Mädler K, Hintze H. Fracture mechanics in railway applications––an overview. Eng Fract Mech 2005;72:163–94. https://doi.org/10.1016/j. engfracmech.2003.11.010. [38] Madia M, Beretta S, Zerbst U. An investigation on the influence of rotary bending and press fitting on stress intensity factors and fatigue crack growth in railway axles. Eng Fract Mech 2008;75:1906–20. https://doi.org/10.1016/j.engfracmech. 2007.08.015. [39] Zerbst U, Vormwald M, Andersch C, Mädler K, Pfuff M. The development of a damage tolerance concept for railway components and its demonstration for a railway axle. Eng Fract Mech 2005;72:209–39. https://doi.org/10.1016/j.engfracmech. 2003.11.011. [40] Zerbst U, Schödel M, Beier HTh. Parameters affecting the damage tolerance behaviour of railway axles. Eng Fract Mech 2011;78:793–809. https://doi.org/10.1016/ j.engfracmech.2010.03.013. [41] NASGRO. Fracture mechanics and fatigue crack growth analysis software. Reference manual; 2002. [42] El Haddad MH, Topper TH, Smith KN. Prediction of non propagating cracks. Eng Fract Mech 1979;11:573–84. https://doi.org/10.1016/0013-7944(79)90081-X. [43] EN 13103:2009+A2:2012. Railway applications – wheelsets and bogies – non powered axles – design method. European Committee for Standardization; 2012. [44] D’Agostino R, Pearson ES. Tests for departure from normality. Empirical results for the distributions of b 2 and √b. Biometrika 1973;60:613–22. https://doi.org/10. 1093/biomet/60.3.613. [45] Paolino DS, Cavatorta MP. Sigmoidal crack growth rate curve: statistical modelling and applications. Fatigue Fract Eng Mater Struct 2013;36:316–26. https://doi.org/ 10.1111/ffe.12001. [46] Annis C. Probabilistic life prediction isn’t as easy as it looks. J ASTM Int 2004;1:1–12. https://doi.org/10.1520/JAI11557.

[1] Zerbst U, Beretta S, Köhler G, Lawton A, Vormwald M, Beier HTh, et al. Safe life and damage tolerance aspects of railway axles – a review. Eng Fract Mech 2013;98:214–71. https://doi.org/10.1016/j.engfracmech.2012.09.029. [2] Zerbst U, Klinger C, Klingbeil D. Structural assessment of railway axles – a critical review. Eng Fail Anal 2013;35:54–65. https://doi.org/10.1016/j.engfailanal.2012. 11.007. [3] Cervello S. Fatigue properties of railway axles: new results of full-scale specimens from Euraxles project. Int J Fatigue 2016;86:2–12. https://doi.org/10.1016/j. ijfatigue.2015.11.028. [4] Novosad M, Fajkoš R, Řeha B, Řezníček R. Fatigue tests of railway axles. Proc Eng 2010;2:2259–68. https://doi.org/10.1016/j.proeng.2010.03.242. [5] Beretta S, Carboni M. Variable amplitude fatigue crack growth in a mild steel for railway axles: experiments and predictive models. Eng Fract Mech 2011;78:848–62. https://doi.org/10.1016/j.engfracmech.2010.11.019. [6] Mädler K, Geburtig T, Ullrich D. An experimental approach to determining the residual lifetimes of wheelset axles on a full-scale wheel-rail roller test rig. Int J Fatigue 2016;86:58–63. https://doi.org/10.1016/j.ijfatigue.2015.06.016. [7] UIC B 169 RP 36:2013-12 Defect Tolerance Concept (DTC) Permitted defects and surface properties for axles under maintenance. UIC Standard; 2013. [8] Pokorný P, Náhlík L, Hutař P. Comparison of different load spectra on residual fatigue lifetime of railway axle. Proc Eng 2014;74:313–6. https://doi.org/10.1016/j. proeng.2014.06.269. [9] Pokorný P, Hutař P, Náhlík L. Residual fatigue lifetime estimation of railway axles for various loading spectra. Theor Appl Fract Mech 2016;82:25–32. https://doi.org/ 10.1016/j.tafmec.2015.06.007. [10] 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 2017;48:687–98. https://doi.org/10.1002/ mawe.201600570. [11] Beretta S, Regazzi D. Probabilistic fatigue assessment for railway axles and derivation of a simple format for damage calculations. Int J Fatigue 2016;86:13–23. https://doi.org/10.1016/j.ijfatigue.2015.08.010. [12] Wu SC, Xu ZW, Kang GZ, He WF. Probabilistic fatigue assessment for high-speed railway axles due to foreign object damages. Int J Fatigue 2018;117:90–100. https://doi.org/10.1016/j.ijfatigue.2018.08.011. [13] Paris PC, Erdogan F. A critical analysis of crack propagation laws. J Basic Eng 1963;85:528–33. https://doi.org/10.1115/1.3656900. [14] Akama M, Ishizuka H. Reliability analysis of Shinkhansen vehicle axle using probabilistic fracture mechanics. JSME Int J Ser Mech Mater Eng 1995;38(3):378–83. https://doi.org/10.1299/jsmea1993.38.3_378. [15] Bea JA. Simulación del crecimiento de grietas por fatiga aleatoria mediante elementos probabilistas. PhD thesis. Universidad de Zaragoza; 1997. [16] Bea JA, Doblaré M, Gracia L. Evaluation of the probability distribution of crack propagation life in metal fatigue by means of probabilistic finite element method and B-models. Eng Fract Mech 1999;63:675–711. https://doi.org/10.1016/S00137944(99)00053-3. [17] Hillmansen S, Smith RA. The management of fatigue crack growth in railway axles. Proc Inst Mech Eng Part F J Rail Rapid Transit 2004;218:327–36. https://doi.org/ 10.1243/0954409043125879. [18] Hong YJ, Xing J, Wang JB. A second-order third-moment method for calculating the reliability of fatigue. Int J Press Vessels Pip 1999;76:567–70. https://doi.org/10. 1016/S0308-0161(99)00013-7. [19] 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 2017;73:139–57. https://doi.org/10.1016/ j.engfailanal.2016.12.014. [20] Wang L, Liang J, Yang Y, Zheng Y. Time-dependent reliability assessment of fatigue crack growth modeling based on perturbation series expansions and interval mathematics. Theor Appl Fract Mech 2018;95:104–18. https://doi.org/10.1016/j. tafmec.2018.02.010. [21] Yasniy O, Lapusta Y, Pyndus Y, Sorochak A, Yasniy V. Assessment of lifetime of railway axle. Int J Fatigue 2013;50:40–6. https://doi.org/10.1016/j.ijfatigue.2012. 04.008. [22] Zhu S-P, Liu Q, Huang H-Z. Probabilistic modeling of damage accumulation for fatigue reliability analysis. Proc Struct Integr 2017;4:3–10. https://doi.org/10. 1016/j.prostr.2017.07.012.

12