Investigation of Fatigue Crack Growth in Full-Scale Railway Axles Subjected to Service Load Spectra: Experiments and Predictive Models

In this paper, a series of experimental investigations was performed on full-scale railway axles to analyze the fatigue crack growth behavior of EA4T steel under load spectrum derived from real operating conditions. The experimental results were compared to life predictions carried out adopting two models: (i) the conventional Nasgro equation and (ii) the cyclic R-curve concept implemented in the Modified Nasgro equation for describing the crack growth behavior of an arbitrary crack length. The results show that the life predictions performed by means of the Modified Nasgro equation coincide well with the experimental results with an underestimation of the residual lifetime less than 32%, while the traditional Nasgro equation leads to significant overestimation (≈120%) of the residual lifetime for load spectra close to the in service scenario.


Introduction
Railway axles are critical components, and their design assessment procedure, as well as the maintenance program, require special attention. The design route must consider the deterioration of the axles during service that can be promoted by (i) surface damage from ballast impacts [1][2][3], (ii) pitting from corrosion [4], and (iii) surface scars in press-fits [5][6][7]. All these types of damages can lead to fatigue crack growth under a sufficient level of local cyclic stresses. Accordingly, the traditional design approach based on the infinite life has to be complemented by a damage tolerant approach. Damage tolerance was originally developed for aeronautic applications in the seventies [8,9] and more recently introduced in the railway industry [10,11]. The key concept of damage tolerance is the assumption of pre-existing defects which can be generally treated based on suitable probabilistic approaches [2,3]. By means of appropriate crack growth algorithms, the residual lifetime and, consequently, an adequate inspection plan aimed to improve the reliability of the component can be defined. Establishing a quantitative procedure for determining residual axle lifetimes is of fundamental importance, and it is still the subject of scientific studies as demonstrated by several projects devoted to this issue, e.g., Widem [12], Euraxles [13], and Maraxil [14].
Numerous studies have addressed the experimental investigation of the fatigue crack growth behaviour of full-scale railway axles [15][16][17][18][19][20][21]. As a result, the need for a proper fatigue crack propagation model applicable for the in-service condition of railway axles is essential to plan the inspection intervals. Several studies have presented damage tolerance analyses of railway axles exploiting the relationship between the fatigue crack growth rate, da/dN, and the stress intensity factor (SIF) range ∆K acquired under constant and variable amplitude load histories [10,[22][23][24][25][26]. The Nasgro equation [27] has been widely used in residual lifetime prediction of the railway axles [28] because it allows a good description of the near-threshold regime, which is the most important region for an axle assessment [10]. Moreover, the Nasgro equation also allows describing the crack growth rates for positive and negative load ratios (R). This effect was shown to be very important for modeling the detrimental effect of positive mean stresses induced in notches by the near press-fits [29], the beneficial effect of cold-rolling [30,31], and induction-hardening [21,32].
However, there are two major limitations related to the adoption of the Nasgro equation for residual lifetime predictions of railway axles. First, the Nasgro equation does not adequately describe short crack growth below the long crack SIF threshold range ∆K th,lc . To overcome this limitation, a modified version of the Nasgro equation has been introduced by J. Maierhofer et al. [33]. This equation can describe the growth rate of arbitrary crack lengths under the condition of small-scale yielding adopting the concept of the cyclic R-curve [34][35][36]. The R-curve highlights the ∆K th dependency with crack extension ∆a, and this represents an attractive improvement for the conventional Nasgro equation [37] because it allows to obtain a physically-based model for short crack growth behaviour. Second, a recent study by P. Pokorny [38] revealed that the opening function f implemented in Nasgro equation is not able to sufficiently describe the crack closure in the threshold area. The results indicated that the roughness-induced and oxide-induced crack closure mechanisms play an important role in the threshold area of the EA4T steel and that R-curve concept leads to an adequate description of the phenomena [39]. In terms of application to railway axles, the differences between the adoption of Nasgro equation and R-curve concept for describing crack growth curves have been recently discussed [40], but an experimental full-scale validation adopting real load spectrum, and a comparison with the widely used Nasgro equation, is still missing.
This manuscript presents an extensive experimental investigation of the crack growth behaviour of full-scale railway axles subjected to variable amplitude loadings under load spectra representative of real service conditions. Then, the application of the conventional Nasgro model vs R-curve concept for evaluating the residual lifetime of the axle is discussed, and the life predictions are compared with the full-scale experiments. For the first time in the literature, results show non-conservative predictions of the traditional Nasgro model, respect to simulations based on the cyclic R-curve, at stress spectra close to real service scenario of railway axles.

Specimens and Test Details of Full-Scale Axle Tests
The adopted full-scale axle specimens are shown in Figure 1. The axles were designed to be tested under the three-point bending bench with 250 kN load capacity (Figure 1b,c), available at both Politecnico di Milano and Lucchini Rs R&D laboratories. The axles were made of 25CrMo4 (EA4T) steel according to the European standard BS EN 13261 [41], and they were quenched and tempered with the tensile properties shown in Table 1. Note that the tensile specimens were machined according to the EN 13261 standard from different locations of the axle. Two hubs representing the wheel and the brake were press-fitted to the axle with contact pressures similar to the actual full-scale components. Semi-elliptical artificial defects with initial depth a = 1.5 mm and aspect ratio a/c = 0.8 were introduced by electrical discharge machining (EDM), being c the superficial crack length. The semi-elliptical shape is dictated by the type of damage that is typically promoted by the impact of foreign objects, as shown and discussed also in the literature [1,2,42]. The dimensions, instead, are determined by the minimum detectable defect size of the employed non-destructive methodology [43] together with the statistics obtained in the RSSB T728 project [2].
The identification of the region where a prospective crack would determine the shortest fatigue life is critical for a reliable assessment of railway axles. Analysis of the in-service fatigue failure of the axles [44][45][46] has shown that the notch transitions and the axle body are the most susceptible areas for fatigue failure. Then, for the present axle geometry, it can be expected that the critical position would occur in one of the notches close to the press-fitted wheel due to the combination of (i) the maximum bending moment, (ii) the change of the axle diameter which induces a notch effect, and (iii) the effect of the press-fit. For these reasons, the defects were positioned in the section near the press-fitted wheel and brake disc (see Figure 1a).  The tests on the full-scale axles were interrupted to measure the crack extension by using the combination of three techniques: the plastic replica, the magnetic particle inspection method, and the optical microscope, examples are shown in Figure 2.
Test interruptions were always made at the occurrence minimum load during test sequences. The plastic replica enables capturing short crack extensions from the notch root. However, when the crack extension is higher, the crack length measurements were performed with the optical microscope. Finally, when the total superficial crack length reaches approximately 40 mm, the magnetic particle inspection method was adopted. The axles were subjected to variable amplitude loading sequences derived from a typical in-service load spectrum train axle applications and representative of 57,000 km service. The load spectrum and its discretization as load blocks are shown in Figure 3. For each test, a variant of this spectrum is employed by increasing the load classes with a given amplification factor. The crack initiation phase of the test, in which the crack extension ∆a was less than 200 µm, was disregarded and only the part of the test where the crack exceeded this length was considered.

Experimental Setup and Test Details of Fatigue Crack Propagation Tests
The single-edge notched bending (SEB) specimen geometry was adopted to determine the crack propagation properties. The present SEB geometry was characterized by a length of 110 mm, thickness of 12 mm, a width of 24 mm, and a notch depth of 6.5 mm. The specimens were all cut from the tested full-scale axles. To nucleate a crack with limited load history effects, the compression pre-cracking procedure was adopted [47,48]. This technique enabled nucleating pre-cracks with lengths in the order of 100 to 200 µm.
Following pre-cracking, a combination of different testing procedures were adopted to evaluate the fatigue crack growth rates, the long crack thresholds, and cyclic R-curve: 1. Compression Pre-cracking Constant Amplitude (CPCA) [49] 2. Compression Pre-cracking Load Reduction (CPLR) [48] 3. Compression Pre-cracking Constant-∆F and Constant-∆K methods used for the determination of the cyclic R-curve [50].
The tests were performed in bending configuration under a Rumul Craktronic resonant bending load frame characterized by a working frequency of approximately 130 Hz. A summary of the experimental campaign performed on the fatigue crack propagation specimens and the experimental techniques adopted for crack growth measurements is presented in Table 2. Note that an additional 9 specimens were used to determine the R-curve at R = −1 (further details related to these experiments are given in [50]).

Results of the Full-Scale Tests
The results of the full-scale tests were analyzed in terms of crack length vs. number of cycles. Table 3 summarizes all the full-scale crack propagation tests that were performed in the framework of this research. The table contains the calculated maximum ∆K by adopting the Wang-Lambert [51] weight function SIF solution at surface and deepest points normalized by the average ∆K th,lc obtained at R = −1 (see Table 2). Table 3 also provides an indication of the load classes which would yield to crack propagation according to a simple crack propagation model that adopts ∆K th,lc as the threshold condition. It can be seen that the highest level of the spectrum was applied for axle S064 with 2.3% and 35% of the spectrum above the ∆K th,lc for the deepest and surface points of the initial crack, respectively. In comparison, the lowest level of the spectrum was applied for axle S067, where 99.98% and 98.1% of the spectrum at the points a and c were below the ∆K th,lc at the beginning of the crack propagation phase. Following the full-scale tests, the axles were cut and broken to analyze the fracture surfaces. As an example, Figure 4 shows the fracture surface of the S064 axle. During crack propagation, the crack depth to surface length ratio a/c varies, but the basic shape remains semi-elliptical, as also observed in other experiments on railway axles [10,24]. As the crack growth during the tests was measured on the surface, the beach marks in the fracture surface were reconstructed using ellipses to assess the crack shape and subsequent crack depth a.

Fatigue Crack Propagation Tests Results
The trend of the ∆K th,lc against the stress ratio is reported in Figure 5a. Figure 5b shows the obtained experimental crack growth data for the load ratios R = −2, −1, 0.05, 0.7. The experimental cyclic R-curve at the load ratio R = −1 is illustrated in Figure 6a. The R-curve shows a steep increase in the threshold up to crack extension of 1 mm, then it gradually increases up to 3 mm and reaches asymptotic ∆K th,lc values. The fitting of the cyclic R-curve will be presented in Section 4.2.
(a) (b) Figure 5. Fatigue crack growth characterization of the present EA4T steel: (a) Evaluated Rdependence of long crack thresholds ∆K th,lc (values acquired from [42]) and comparison with the available data in literature [40] and (b) crack growth curves and interpolation with the Nasgro equation.

Crack Growth Models and Parameter Determination
The conventional Nasgro equation [52] and the Modified Nasgro equation for physically short cracks [33] have been chosen because they take into account different effects acting on the crack propagation rate in components made of metallic materials (i.e., R and load interaction effects). In the following sections, the calibration of the two models is presented and discussed.

Calibration of the Nasgro Equation
The Nasgro equation fully describes the three regions of the da/dN vs. ∆K diagram: where C, n, p, and q are empirically derived constants; ∆K th is the SIF threshold range; and K c is the critical SIF. The parameter f accounts for the plasticity-induced crack closure [53], and it is calculated with the Newman function provided in the Nasgro equation. ∆K th is determined from the following equations: ∆K 1 is the threshold SIF range for R → 1, A 0 is the parameter used in the crack closure equation of Newman, and C th+ and C th− are two empirical constants which control the ∆K th dependency on R and need to be considered separately for positive and negative R values. a 0 is the El-Haddad parameter for small cracks [54], its value was taken from previous analyses on EA4T [5]. Figure 5a,b shows the experimental ∆K th,lc and fatigue crack growth curves obtained by the LR and CA methods (see Section 3.2), respectively. The Nasgro Equations (1) and (2) were fitted applying the least-squares method considering S max /S o = 0.3 and α = 2.5. The present fit is compared with Nasgro fit obtained by EBFW3 project [40] and it can be seen that the results are similar.

Calibration of the Modified Nasgro Equation
The Nasgro equation (Equation (1)) was modified by Mairhofer and co-authors to improve the accuracy of the predicted crack growth rates in the short crack regime [33]. In particular, the authors proposed to adopt the concept of the cyclic R-curve for the description of the ∆K th -dependency with crack length. Accordingly, the initial ∆K th at a nominal crack extension of ∆a = 0 is considered to be equal to the intrinsic SIF range ∆K th,e f f . This condition is supposed to correspond to the initial condition of a defect in the material. The stabilized ∆K th for long cracks is ∆K th,lc , while the transition region is described with the following equation: with the constraint: The l i values can be interpreted as fictitious length scales for the formation of crack closure effects which must be determined in conjunction with the ν i by fitting the experiment R-curve data points [50].
For calculating the crack growth rate for an arbitrary crack length, the Nasgro equation (Equation (1)) was simplified as suggested in [33] by setting p = n and q = 1 independent of R, and it was transformed to The crack growth constant C and the Paris exponent n are determined by least-squares fitting of Equation (5), and the crack velocity factor F is calculated using the predetermined l i and ν i parameters: with in which f is the corresponding to the Newman function in Nasgro equation. Figure 6a illustrates the experimentally obtained R-curve at R = −1 (see Section 3.2) for the present EA4T steel together with the R-curve fit using Equation (3) and a 90% scatter band. To validate the Modified Nasgro model, a comparison of the experimental tests on SEB specimens and predicted crack growth rates for a short crack starting from a notch of depth of 6.5 mm at R = −1 was performed and is presented in Figure 6b. As it can be seen, in the first two load steps the crack growth rate slows down and the crack stops due to the build-up of crack closure. Only after a further increase of the load amplitude, the crack grows reaching the threshold of long cracks.
Additionally, the build-up of crack closure can be assumed to be similar at different R [33]. Following this concept, for modeling the R-curve at different stress ratios, the term (∆K th,lc − ∆K th,e f f ) in Equation (3) was described through Nasgro dependence of ∆K th,lc on R (see Equation (2) for a >> a 0 ).

Residual Lifetime Predictions
For a precise calculation of the residual lifetime of an axle, the knowledge of several factors is required: the load spectrum acting on the axle during service, the crack growth behavior of the adopted steel grade, as well as the precise evaluation of the stress intensity factors along the crack front. In particular, for analyzing semi-elliptical crack propagation, it is necessary to perform separate analyses for the surface point a and for the deepest point c of the crack front.
The Wang-Lambert [51] weight function adopted is valid for a simple plate geometry which can still be used here as it has been proved to give reasonable results when applied to full-scale axles containing semi-elliptical fatigue cracks [55] compared to other solutions [56][57][58][59][60][61]. As shown in Figure 7, the longitudinal normal stresses acting on the crack plane associated with both the rotating bending and the press-fit stresses were ob-tained by a dedicated finite element model (more details about the finite element analysis can be found in [55]). Due to the press-fit, the wheel squeezes the axle under its seat, and the surrounding regions stretch inducing local tensile stresses which are superimposed with the rotating bending stresses. As the tensile residual stresses are constant, they determine an increment of the local R which will be higher than −1 [29].
(a) (b) Figure 7. Stress analysis in the notch section for applied nominal stress of 130 MPa (rearranged from the work in [55]): (a) longitudinal stresses induced by the press-fit and (b) calculated in-depth stress profiles due to the rotating bending and press-fit and manufacturing induced compressive residual stresses (the value of the measured scatter marked with grey is acquired from [62]; the value of residual stress marked as black line is acquired from [62]).
There are additional compressive residual stresses acting on the crack plane induced by the manufacturing process. As it is shown in [63], the residual stresses determine a positive or a negative effect on the residual lifetime depending on the sign. The majority of the residual lifetime estimations studies consider the axle loadings as determined by the rotating bending (train weight + dynamic effects of moving train) and the press fit, by neglecting the compressive residual stresses effect at the axle surface [22,29,64,65]. These stresses are observed to have a positive effect determining a residual fatigue life substantially longer [30,[66][67][68]. In the recent study by P. Pokorný et al. [69], it was demonstrated that in the case of a standardly manufactured railway axle adopting EA4T steel with consequent heat treatment and machining, it is possible to develop negative residual stresses up to depth of 5 mm inducing a significant effect on the residual lifetime.
The residual stresses were measured in correspondence of the section containing the artificial notch up to depth of 2 mm, using the hole drilling method [70]. The measurements were derived for three surface points (respectively at 0 • , 200 • and 300 • ). The average of experimental residual stresses is indicated in Figure 7b with green triangle dots, while the light green shaded area indicates the dispersion of the performed measurements. The obtained residual stresses were compared with the available data by Pokorny [69] and a fair agreement was found with the present results. Therefore, for the sake of simulation of test results, the lower bound profile measured by Pokorny is considered.

Crack Growth Algorithm
The life predictions are based upon a crack growth algorithm which considers the discretized load spectrum shown in the Figure 3. A summary of the algorithm steps is given in the following: 1. Given a block load level, the local load ratio R acting at the points a and c of the crack is calculated by using Equation (8).
In the case of superposition of rotating bending and residual stresses, the effective R is changing between −1.75 < R < −0.7 along the crack propagation path. On the other hand, for the analysis carried out considering only the rotating bending stresses, the load ratio is always R = −1. 2. Evaluation of the SIFs acting at the two points a and c adopting the Wang-Lambert weight function and the normal stress distribution perpendicular to the crack plane for the un-cracked axle calculated with the FE model. 3. Given the stress intensity factor ranges ∆K and the R, for both points a and c, evaluation of the crack growth rates using the conventional Nasgro equation and Modified Nasgro equation (see Section 4). The evaluated crack growth rates for axle S067 has been shown in Figure 8. According to the modified Nasgro equation, at the beginning of the simulation the crack growth rates tend to decrease with crack propagation due to crack closure development. Finally, when the crack closure fully develops and the driving force ∆K increases due to the crack advancement, the crack growth rates increase and converge along the Paris region. 4. Updating the actual crack size and shape, which acted as an input for the following load block. 5. Simulations are run until a crack depth of 20 mm is reached.

Results and Discussion
The residual lifetime predictions for each axle are presented in Figure 9 and Table 4. The experimental data are reported with blue dots, while the predictions performed with the Nasgro equation and the Modified Nasgro equation are indicated, respectively, with the green and red colors. In addition, the solid lines indicate the predictions performed considering only the rotating bending stresses, while the dashed lines indicate the predictions considering the superposition of rotating bending stresses and residual stresses. The first comparison with experiments can be made considering the life predictions based on the stresses induced by the rotating bending (solid lines) and the ones with bending plus residual stresses (dashed lines). It can be seen, as already noted in the literature, that residual stresses have a significant effect on life prediction. In terms of comparison between the algorithms, the simulations performed with the Modified Nasgro equation are in good agreement with the experimental results. Even if they are slightly conservative (maximum underestimation 30%), it is clear that the Modified Nasgro equation produces predictions with a curvature very close to the experimental results. Remarkably, the predictions are very precise under load spectra with low stresses (Figure 9c,d) that are the ones closer to the real service conditions. By contrast, it is evident that the adoption of the conventional Nasgro equation generally leads to a significant overestimation of the propagation lifetime of the order of 120% (see Section 6.2 for an explanation of this effect).
Even if recent papers have shown the possibility to include interaction effects on propagation simulations based on the R-curve [40], the comparison with the present experiments show that, at the stress levels experienced by axles in service, neglecting those effects does not significantly affects the axle life predictions. This is surely in line with the modest interaction effects already observed for EA4T [71]. In fact, analyses based on the Strip-Yield model [62,72] showed that this effect was due to a small increase of S op under variable amplitude loads.

Crack Shape Evolution
Another interesting analysis is shown in Figure 10, where the experimental crack shape evolution characterized by the shape factor a/c is compared with the numerical predictions (the results refer to the Axle S064). The experimental a/c values show small fluctuations around an average value of approximately 0.8. From the crack propagation simulation results, it can be observed that neglecting the residual stresses determines a continuous decrease of the aspect ratio a/c from the beginning of the simulations, either performed with the Nasgro or Modified Nasgro equations. This tendency is also captured for a/D > 0.05 when considering the residual stresses in the simulation (green and dashed lines in the Figure 10). However, the residual stresses have a paramount effect on the shape factor evolution in the first stages of propagation. In fact, due to the trend of the compressive residual stresses, which displays higher absolute values at 1 mm depth compared to the free surface, the initial trend of the simulated aspect ratio shows a quick drop of a/c (crack forced to propagate faster at the free surface) and a successive increment. Interestingly, both the simulations performed with the Nasgro and the Modified Nasgro equations indicate that a maximum a/c is reached at a crack depth of approximately a/D = 0.05. Successively, the simulations display a trend of a/c which decreases according the same rate as it is observed for the simulations performed neglecting the residual stresses. From this analysis, it can be concluded that the crack shape evolution is also well captured with the crack propagation algorithm based on the Modified Nasgro equation as evidenced comparing the experimental a/c values (blue dots) and the simulated ones (dashed red line) in the Figure 10a,b.
(a) (b) Figure 10. Crack shape evolution analysis of Axle S067: (a) Crack shape aspect ratio a/c comparison between numerical prediction and experimental and (b) contour plot of the crack shape evaluated by considering residual stresses the comparison with the experimental evidence.

Effect of Closure Development
The significant overestimation of the residual lifetime obtained adopting the conventional Nasgro algorithm in the axles tested at low stresses can be clearly understood by analyzing the stress spectrum in terms of ∆K for the initial crack. Figure 11a illustrates the calculated driving forces at the surface point of the initial defect (a 0 = 1.5 mm and c 0 = 1.87 mm) with a ∆a = 200 µm crack extension ahead of the notch for axle 067. Figure 11 clearly highlights the difference between the two models. In particular, Figure 11a shows the load blocks involved in crack propagation according to the Nasgro equation (green bars) and the Modified Nasgro equation (red and green bars combined), which correspond to 1.9% and 100% of the overall spectrum, respectively.
Considering that the initial crack depth is 1.7 mm (a 0 + ∆a), according to the El-Haddad model, this crack can already be considered as a long crack with ∆K th ≈ ∆K th,LC . On the other side, according to the R-curve, at the beginning of the crack propagation the threshold ∆K th is much lower with respect to the long crack threshold ∆K th,LC . Consequently, for the El-Haddad model, only the loading cycles highlighted with the green bars (see Figure 11b) are above the threshold condition, while, according to the R-curve, all the loading cycles contribute to the early stages of crack propagation. Even if this difference was already discussed in terms of fatigue limit for small cracks/defects [73], it is the first time that its effect has been discussed and experimentally validated through components' experimental data. The dramatic consequence is that the lifetime predictionbased on the traditional Nasgro algorithm is potentially unsafe at stress levels close to real service conditions.

Concluding Remarks
Wheelset requires optimization of the inspection intervals to reduce the service life maintenance costs. The primary tool needed for reaching this objective is developing a robust residual lifetime prediction model which is applicable to realistic operational conditions of the wheelset. In this paper, series of experimental investigations were carried out on full-scale railway axles to analyze the fatigue crack growth behavior of the EA4T steel subjected to real operational load spectrums. Two crack propagation algorithms (Nasgro equation and Modified Nasgro equation) have been applied to estimate the propagation lifetime of a fatigue crack in railway axles and a comparison of predictions with a set of 4 full-scale block-loading experiments was presented. The main results can be summarized as follows. • Adopting the concept of cyclic R-curve through the Modified Nasgro equation leads to axle life predictions that are in good agreement with the full-scale experimental results (max underestimation less than 30%) if all the contribution to mean stresses (press-fit stresses + residual stresses) are included in the analysis. • The conventional Nasgro equation leads to significant overestimation of propagation lifetime, approximately by a factor 2.2, in the real service scenario for railway axles when the majority of the load classes in the load spectrum corresponds to ∆K which are below the long crack threshold.

Data Availability Statement:
The data published in this paper can be requested to the corresponding author.

Acknowledgments:
The present results were obtained with a research contract between LucchiniRS and Politecnico di Milano, Dept. Mechanical Engineering. The authors express gratitude to Luc-chiniRS, especially S. Cantini and S. Cervello, for permission to publish the results here contained.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data. The funders accepted authors' decision to publish the results in the present paper.

Abbreviations
The following abbreviations are used in this manuscript: SIF Stress intensity factor R Load (or stress) ratio EDM Electrical discharge machining EA4T 25CrMo4 steel CPCA Compression pre-cracking constant-amplitude CPLR Compression pre-cracking Load reduction