Fatigue Analysis and Defect Size Evaluation of Filled NBR including Temperature Influence

The fatigue behavior of a filled non-crystallizing elastomer was investigated on axisymmetric dumbbell specimens. By plotting relevant Wöhler curves, a power law behavior was found. In addition, temperature increases due to heat build-up were monitored. In order to distinguish between initiation and crack growth regimes, hysteresis curves, secant and dynamic moduli, dissipated and stored energies, and normalized minimum and maximum forces were analyzed. Even though indications related to material damaging were observed, a clear trend to recognize the initiation was not evident. Further details were revealed by considering a fracture mechanics. The analysis of the fracture surfaces evidenced the presence of three regions, associated to initiation, fatigue striation, and catastrophic failure. Additional fatigue tests were performed with samples in which a radial notch was introduced. This resulted in a reduction in lifetime by four orders of magnitude; nevertheless, the fracture surfaces revealed similar failure mechanisms. A fracture mechanics approach, which considered the effect of temperature, was adopted to calculate the critical defect size for fatigue, which was found to be approximately 9 μm. This value was then compared with the particle size distribution obtained through X-ray microcomputed tomography (μ-CT) of undamaged samples and it was found that the majority of the initial defects were indeed smaller than the calculated one. Finally, the evaluation of J-integral for both unnotched and notched dumbbells enabled the assessment of a geometry-independent correlation with fatigue life.


Introduction
Fatigue is one of the most common causes of the failure of elastomer components such as tires, dampers, seals, hydraulic hoses, and conveyor belts. In fact, rubber articles are commonly subjected to cyclic loads that may lead to the accumulation of damage and to the development and growth of cracks, which then lead to failure. Generally speaking, materials are assumed to naturally contain flaws and inhomogeneities that can lead, under specific loading conditions, to localized stress concentrations promoting the nucleation of cracks [1]. Indeed, fracture occurs due to defects, which may nucleate from inhomogeneities that grow when subjected to repeated loads and subsequently coalesce up to failure. Similarly, cyclic tests were performed on notched circular dumbbell specimens, denoted as Crack Round Bars-CRB (Figure 1b). According to ISO 18489 [36], an initial circumferential notch of approximately 1 mm depth (Figure 1c) was inserted at midheight using a razor blade (thickness = 0.1 mm, tip-radius < 5 μm) mounted on a lathe rotating at 80 rpm. The specimens were cyclically loaded at maximum strains of 20, 30, and 50%, using the same frequency and load ratio used for the unnotched samples. Furthermore, the tests were monitored using a camera system CV-5701P by Keyence (Osaka, Japan) and pictures were recorded at every cycle.

Fracture Surface Analysis
An optical microscope Zeiss Stemi 2000C (Carl Zeiss AG, Oberkochen, Germany) was used for the fracture surface analysis.

X-ray Microtomography (μ-CT)
μ-CT scans of the dumbbell samples were acquired on a NSI X-25 inspection system (North Star Imaging, High Wycombe, UK). The X-ray source was operated at 40 kV and 100 μA, granting a magnification factor of 3.88 and a voxel resolution of 19 μm. It should be considered that the apparatus is capable of a nominal certified resolution of 2 μm [37]. The set of images was processed with a custom-made Matlab code (Matlab R2020b), performing binarization based on a combined global/local algorithm procedure [38] as a first step. The global algorithm adopted is the one developed by Otsu [39,40], which is able to distinguish between measurement artefacts and particles that are actually present inside the inspected specimen. Then, a local gradient-based watershed algorithm [41] was employed to better define the contour of the inspected particles. Following successful segmentation (i.e., separation of the voxels representing rubber and void, with Boolean values of 1 and 0, respectively), the code allowed particle detection by exploiting a spatial connectivity equal to 26, which means that two or more voxels belong to the same particle if they are in contact with at least one vertex. Each particle was then labelled, and its volume was computed as the sum of the connected full voxels (Boolean value equal to 1). Similarly, cyclic tests were performed on notched circular dumbbell specimens, denoted as Crack Round Bars-CRB (Figure 1b). According to ISO 18489 [36], an initial circumferential notch of approximately 1 mm depth (Figure 1c) was inserted at mid-height using a razor blade (thickness = 0.1 mm, tip-radius < 5 µm) mounted on a lathe rotating at 80 rpm. The specimens were cyclically loaded at maximum strains of 20, 30, and 50%, using the same frequency and load ratio used for the unnotched samples. Furthermore, the tests were monitored using a camera system CV-5701P by Keyence (Osaka, Japan) and pictures were recorded at every cycle.

Fracture Surface Analysis
An optical microscope Zeiss Stemi 2000C (Carl Zeiss AG, Oberkochen, Germany) was used for the fracture surface analysis.

X-ray Microtomography (µ-CT)
µ-CT scans of the dumbbell samples were acquired on a NSI X-25 inspection system (North Star Imaging, High Wycombe, UK). The X-ray source was operated at 40 kV and 100 µA, granting a magnification factor of 3.88 and a voxel resolution of 19 µm. It should be considered that the apparatus is capable of a nominal certified resolution of 2 µm [37]. The set of images was processed with a custom-made Matlab code (Matlab R2020b), performing binarization based on a combined global/local algorithm procedure [38] as a first step. The global algorithm adopted is the one developed by Otsu [39,40], which is able to distinguish between measurement artefacts and particles that are actually present inside the inspected specimen. Then, a local gradient-based watershed algorithm [41] was employed to better define the contour of the inspected particles. Following successful segmentation (i.e., separation of the voxels representing rubber and void, with Boolean values of 1 and 0, respectively), the code allowed particle detection by exploiting a spatial connectivity equal to 26, which means that two or more voxels belong to the same particle if they are in contact with at least one vertex. Each particle was then labelled, and its volume was computed as the sum of the connected full voxels (Boolean value equal to 1).

Evalaution of J-Integral Using Finite Element Method (FEM)
FEM simulation of the CRB specimen behavior was carried out using a 2D axisymmetric model created in Ansys Parametric Design Language (APDL). The element type PLANE 183 was used, which is a higher order 2D, 8-node, or 6-node element [42]. Mapped mesh was used where possible. The mesh was made very fine in the vicinity of the crack tip with elements as small as 0.01 mm in the area closest to the crack tip. Moreover, special wedge-shaped crack tip elements were used at the crack tip. The rest of the model was meshed with a coarser mesh, with the largest elements up to 1.2 mm. The total number of elements was approximately 7000 for each of the modelled crack lengths. See Figure 2 where the mesh of the model for the crack length of 2.5 mm is depicted. A neo-Hookean incompressible hyperelastic solid was chosen for the material model; the shear modulusobtained by fitting the uniaxial tensile test data-was 12.3 MPa, whereas the bulk modulus was not acquired as the material was assumed to be incompressible. Even though the rubbers are not strictly incompressible, this simplification does not make a big difference for the intended analysis. The goal of the simulation was to calculate values of the J-integral for various crack lengths and various loads. The calculation of the non-linear J-integral was performed through closed path integral calculation (CINT function in APDL).

Evalaution of J-Integral Using Finite Element Method (FEM)
FEM simulation of the CRB specimen behavior was carried out using a 2D axisymmetric model created in Ansys Parametric Design Language (APDL). The element type PLANE 183 was used, which is a higher order 2D, 8-node, or 6-node element [42]. Mapped mesh was used where possible. The mesh was made very fine in the vicinity of the crack tip with elements as small as 0.01 mm in the area closest to the crack tip. Moreover, special wedge-shaped crack tip elements were used at the crack tip. The rest of the model was meshed with a coarser mesh, with the largest elements up to 1.2 mm. The total number of elements was approximately 7000 for each of the modelled crack lengths. See Figure 2 where the mesh of the model for the crack length of 2.5 mm is depicted. A neo-Hookean incompressible hyperelastic solid was chosen for the material model; the shear modulus-obtained by fitting the uniaxial tensile test data-was 12.3 MPa, whereas the bulk modulus was not acquired as the material was assumed to be incompressible. Even though the rubbers are not strictly incompressible, this simplification does not make a big difference for the intended analysis. The goal of the simulation was to calculate values of the J-integral for various crack lengths and various loads. The calculation of the non-linear J-integral was performed through closed path integral calculation (CINT function in APDL).

Fatigue of the Unnotched Dumbbell Specimens
The fatigue lifetime of the specimens subjected to cyclic loading is reported in Figure  3, in which the relevant ε-N or Wöhler curve [43] is represented. In particular, the maximum engineering strain εmax was selected as the representative loading parameter and plotted in double logarithmic scale against the number of cycles to failure, Nf.

Fatigue of the Unnotched Dumbbell Specimens
The fatigue lifetime of the specimens subjected to cyclic loading is reported in Figure 3, in which the relevant ε-N or Wöhler curve [43] is represented. In particular, the maximum engineering strain ε max was selected as the representative loading parameter and plotted in double logarithmic scale against the number of cycles to failure, N f .
As depicted in the plot, a certain scattering of the results was found, which is typical for the fatigue of polymers and rubbers. The statistical analysis of the results was performed through the ASTM E739-10 [44], and the resulting fitting power law was found to be: Furthermore, two interrupted tests, for which no failure was recorded, are shown in Figure 3, and they correspond to ε max of 5% and 20%. In these two cases, the loads were so low that the failure would occur at a number of cycles that is well beyond the experimental window.
An alternative criterion for studying the fatigue nucleation of rubbers is to consider N i , corresponding to the number of cycles for the transition between the phase of initiation and crack growth. This is normally related to an abrupt variation of force (or stiffness) [45]. By applying this concept to the results previously shown, it was found that qualitatively, the plot was almost identical to the one reported in Figure 3 (provided in the Supplementary Material, Figure S1). In fact, the abrupt variation of force took place only in the last cycles (approximately 10 to 20 cycles). In actuality, this crack propagation phase corresponds to the macroscopic crack growth, which is visible to the naked eye. Generally, after the initial formation of cracks from inhomogeneities in the material, a smallscale growth of microcracks is present until macroscopic crack growth occurs, leading to a reduction of the ligament area and rapid failure.  (1)) was obtained using ASTM E739-10.
As depicted in the plot, a certain scattering of the results was found, which is typical for the fatigue of polymers and rubbers. The statistical analysis of the results was performed through the ASTM E739-10 [44], and the resulting fitting power law was found to be: Furthermore, two interrupted tests, for which no failure was recorded, are shown in Figure 3, and they correspond to εmax of 5% and 20%. In these two cases, the loads were so low that the failure would occur at a number of cycles that is well beyond the experimental window.
An alternative criterion for studying the fatigue nucleation of rubbers is to consider Ni, corresponding to the number of cycles for the transition between the phase of initiation and crack growth. This is normally related to an abrupt variation of force (or stiffness) [45]. By applying this concept to the results previously shown, it was found that qualitatively, the plot was almost identical to the one reported in Figure 3 (provided in the Supplementary Material, Figure S1). In fact, the abrupt variation of force took place only in the last cycles (approximately 10 to 20 cycles). In actuality, this crack propagation phase corresponds to the macroscopic crack growth, which is visible to the naked eye. Generally, after the initial formation of cracks from inhomogeneities in the material, a small-scale growth of microcracks is present until macroscopic crack growth occurs, leading to a reduction of the ligament area and rapid failure.
Further analysis of the fatigue tests led to the inspection of temperature evolution upon cyclic loading. In Figure 4a, the variation of the surface temperature for the first 15,000 cycles is shown. The temperature variation is a consequence of heat build-up; during every loading cycle, part of the mechanical energy is dissipated in the specimen volume and converted into heat. Due to the low thermal conductivity of rubbers [46,47], the internally generated heat slowly flows to the surrounding environment, and therefore, the temperature of the specimen significantly increases. This phenomenon strongly depends on the material compliance, on the frequency of oscillation, on the loading amplitude, and on the thickness of the component [48]. In fact, the rate of temperature variation ∆ ̇ can be described as [49,50]: Further analysis of the fatigue tests led to the inspection of temperature evolution upon cyclic loading. In Figure 4a, the variation of the surface temperature for the first 15,000 cycles is shown. The temperature variation is a consequence of heat build-up; during every loading cycle, part of the mechanical energy is dissipated in the specimen volume and converted into heat. Due to the low thermal conductivity of rubbers [46,47], the internally generated heat slowly flows to the surrounding environment, and therefore, the temperature of the specimen significantly increases. This phenomenon strongly depends on the material compliance, on the frequency of oscillation, on the loading amplitude, and on the thickness of the component [48]. In fact, the rate of temperature variation . ∆T can be described as [49,50]: where f is the frequency, ∆σ is the stress amplitude, D is the loss compliance, c p is the specific heat and ρ is the density; therefore, larger load amplitudes and frequencies will result in larger temperature increases. Indeed, the effect of increasing the loading amplitude is confirmed by the results reported in Figure 4a,b. The monitored temperature increases are consequences of the geometry of the specimen (i.e., the thickness) and the used frequency (4 Hz). In actuality, a decrease of these two variables could lead to the reduction of the temperature increase; however, thicker specimens are more representative for certain applications, and the axisymmetric dumbbell samples (or slight variations of it) represent a well-established geometry for fatigue characterization of rubbers [14,15,18,[51][52][53][54][55][56][57][58][59][60][61][62]. The choice of the frequency has to match the requirement of being representative for the material application, and it is a tradeoff between testing time and temperature increases. Overall, with the used loading conditions, the surface temperatures increased in a range between 20 and 40 • C. Generally, the temperatures showed fluctuations related to the cooling system of the laboratory. Nevertheless, after approximately 2500-3000 cycles, the surface temperature became stable. As depicted in Figure 4b, higher maximum strains corresponded to higher temperatures in an almost linear relationship. The monitored temperature increases were so relevant that they must be considered during the analysis of fatigue results since such temperatures would affect the mechanical behavior of the material. In general, due to the entropic elasticity of the rubber matrix, an increase of stiffness would be expected; however, due to the presence of fillers and due to viscoelastic effects, the material would become softer at higher temperatures. Thus, the overall trend would be a trade-off between these opposite trends, and it would depend on the material composition. A reduction of stiffness at a higher temperature was previously reported for the material under investigation [35]. Nevertheless, the mechanical behavior is expected to be constant after reaching the plateau temperature.
The choice of the frequency has to match the requirement of being representative for the material application, and it is a tradeoff between testing time and temperature increases. Overall, with the used loading conditions, the surface temperatures increased in a range between 20 and 40 °C. Generally, the temperatures showed fluctuations related to the cooling system of the laboratory. Nevertheless, after approximately 2500-3000 cycles, the surface temperature became stable. As depicted in Figure 4b, higher maximum strains corresponded to higher temperatures in an almost linear relationship. The monitored temperature increases were so relevant that they must be considered during the analysis of fatigue results since such temperatures would affect the mechanical behavior of the material. In general, due to the entropic elasticity of the rubber matrix, an increase of stiffness would be expected; however, due to the presence of fillers and due to viscoelastic effects, the material would become softer at higher temperatures. Thus, the overall trend would be a trade-off between these opposite trends, and it would depend on the material composition. A reduction of stiffness at a higher temperature was previously reported for the material under investigation [35]. Nevertheless, the mechanical behavior is expected to be constant after reaching the plateau temperature.  Further analyses of the fatigue data were implemented to investigate the distinct regimes of fatigue, namely, crack initiation and crack propagation. The hysteresis curves during fatigue tests for different strain levels are reported in Figure 5a, in terms of engineering stress and engineering strain at 3000 cycles. As depicted, higher strains corresponded to larger hysteresis. In Figure 5b, the evolution of the hysteresis curves during the tests is shown for a particular specimen (εmax = 55%) for a different number of cycles. Further analyses of the fatigue data were implemented to investigate the distinct regimes of fatigue, namely, crack initiation and crack propagation. The hysteresis curves during fatigue tests for different strain levels are reported in Figure 5a, in terms of engineering stress and engineering strain at 3000 cycles. As depicted, higher strains corresponded to larger hysteresis. In Figure 5b, the evolution of the hysteresis curves during the tests is shown for a particular specimen (ε max = 55%) for a different number of cycles. A big variation of the hysteresis is present in the beginning of the test as a consequence of stress relaxation. The hysteresis seems to stabilize after approximately 3000 cycles, which corresponds to the stabilization of the surface temperatures. During the successive lifetime, the hysteresis variation was much smaller, which was a consequence of the combination of both stress relaxation and material damage accumulation. cycles. A big variation of the hysteresis is present in the beginning of the test as a consequence of stress relaxation. The hysteresis seems to stabilize after approximately 3000 cycles, which corresponds to the stabilization of the surface temperatures. During the successive lifetime, the hysteresis variation was much smaller, which was a consequence of the combination of both stress relaxation and material damage accumulation. For a better understanding of the results, the hysteresis curves were analyzed considering two parameters, the secant Esec and the dynamic Edyn modulus that can be calculated according to Equations (3) and (4) where σmax and σmin represent the maximum and minimum stresses and εmax and εmin are the maximum and minimum strain, respectively. Although the dynamic modulus is correlated with the instantaneous material response, the variation of the secant modulus is related to viscoelastic effects and to the accumulation of damage in the material [64,65].
Equations (3) and (4) are probably not able to fully describe the highly non-linear mechanical behavior of elastomers; however, for a first approximation, their trends for the number of cycles may provide interesting insights about fatigue. Due to the large variation in temperature at the beginning of the tests, both moduli initially showed large drops. For a better comparison, the values were normalized with respect to those at 2500 cycles (i.e., once temperature and mechanical behavior became stable). The results at the different strains are reported in Figure 6. In all cases, it is possible to observe a similar initial drop for both moduli. Subsequently, Edyn showed a lower rate of reduction and became more stable-even slight increases were monitored; therefore, the material stiffness remained essentially constant in the last stages of fatigue. The initial decrease and the successive stabilization were associated with cyclic stress softening [66], which seems to be connected to the breakdown and reformation of filler-molecule interactions [67]. On the other hand, For a better understanding of the results, the hysteresis curves were analyzed considering two parameters, the secant E sec and the dynamic E dyn modulus that can be calculated according to Equations (3) and (4) [63,64]: where σ max and σ min represent the maximum and minimum stresses and ε max and ε min are the maximum and minimum strain, respectively. Although the dynamic modulus is correlated with the instantaneous material response, the variation of the secant modulus is related to viscoelastic effects and to the accumulation of damage in the material [64,65].
Equations (3) and (4) are probably not able to fully describe the highly non-linear mechanical behavior of elastomers; however, for a first approximation, their trends for the number of cycles may provide interesting insights about fatigue. Due to the large variation in temperature at the beginning of the tests, both moduli initially showed large drops. For a better comparison, the values were normalized with respect to those at 2500 cycles (i.e., once temperature and mechanical behavior became stable). The results at the different strains are reported in Figure 6. In all cases, it is possible to observe a similar initial drop for both moduli. Subsequently, E dyn showed a lower rate of reduction and became more stable-even slight increases were monitored; therefore, the material stiffness remained essentially constant in the last stages of fatigue. The initial decrease and the successive stabilization were associated with cyclic stress softening [66], which seems to be connected to the breakdown and reformation of filler-molecule interactions [67]. On the other hand, E sec shows a steady drop for the entire test duration, as a consequence of the continuous stress relaxation and accumulation of damage during the test. Esec shows a steady drop for the entire test duration, as a consequence of the continuous stress relaxation and accumulation of damage during the test. Figure 6. Dynamic (Edyn) and secant (Esec) moduli evaluated from hysteresis at maximum strain of (a) 50%, (b) 55%, (c) 60%, and (d) 65%, respectively. The fatigue tests were implemented on axisymmetric dumbbell specimens of CB filled NBR at a frequency of 4 Hz and a load ratio Rε of 0.5.
Subsequently, the energies were analyzed by integrating the area of the force-displacement curves from the hysteresis curves previously discussed. In this sense, it was possible to distinguish between the dissipated energy, Udis, as the area inside the hysteresis loop, and the stored energy, Usto, as the area below the unloading curve. For a proper comparison, the evolution of the energies was plotted as a function of the normalized number of cycles (currently elapsed number of cycles, N, divided by the number of cycles to failure, Nf) as reported in Figure 7a. These results originate from the behavior already described for the hysteresis curves: in all the testing configurations, a drop in both energies was observed during the first cycles. As previously discussed, this effect is mainly a consequence of stress relaxation. After the drop, the values of Udis are then substantially stabilized for plateau values. The stored energy Usto instead kept decreasing with the number of cycles, albeit slightly, due to energy dissipation caused by material damage. On the other hand, the dissipated energies moderately increased, evidencing the increase in Figure 6. Dynamic (E dyn ) and secant (E sec ) moduli evaluated from hysteresis at maximum strain of (a) 50%, (b) 55%, (c) 60%, and (d) 65%, respectively. The fatigue tests were implemented on axisymmetric dumbbell specimens of CB filled NBR at a frequency of 4 Hz and a load ratio R ε of 0.5.
Subsequently, the energies were analyzed by integrating the area of the force-displacement curves from the hysteresis curves previously discussed. In this sense, it was possible to distinguish between the dissipated energy, U dis , as the area inside the hysteresis loop, and the stored energy, U sto , as the area below the unloading curve. For a proper comparison, the evolution of the energies was plotted as a function of the normalized number of cycles (currently elapsed number of cycles, N, divided by the number of cycles to failure, N f ) as reported in Figure 7a. These results originate from the behavior already described for the hysteresis curves: in all the testing configurations, a drop in both energies was observed during the first cycles. As previously discussed, this effect is mainly a consequence of stress relaxation. After the drop, the values of U dis are then substantially stabilized for plateau values. The stored energy U sto instead kept decreasing with the number of cycles, albeit slightly, due to energy dissipation caused by material damage. On the other hand, the dissipated energies moderately increased, evidencing the increase in damage in the material. Moreover, a clear trend with respect to the maximum strains is present: both dissipated and stored energies were higher for higher strains; however, this description is different when considering the percentage of energy, which is the normalization with respect to the total energy, U tot , calculated as the sum of U dis and U sto . In fact, as reported in Figure 7b, the percentage energies were very similar for different strain levels and for the whole duration of the tests, maintaining values of approximately 20% for the dissipated energy and 80% for the stored one.
Materials 2022, 15, x FOR PEER REVIEW 9 of 23 damage in the material. Moreover, a clear trend with respect to the maximum strains is present: both dissipated and stored energies were higher for higher strains; however, this description is different when considering the percentage of energy, which is the normalization with respect to the total energy, Utot, calculated as the sum of Udis and Usto. In fact, as reported in Figure 7b, the percentage energies were very similar for different strain levels and for the whole duration of the tests, maintaining values of approximately 20% for the dissipated energy and 80% for the stored one. Fatigue data were further analyzed by considering the peak values of the forces detected during the measurements. For a clearer comparison at different loading conditions, the maximum force, Fmax, and the minimum force, Fmin, were normalized with respect to their respective initial values and they were plotted against the normalized number of cycles (Figure 8a,b). Due to stress relaxation, the maximum and minimum forces showed a considerable variation in the first cycles, until consistent stress relaxation occurs in the specimen. After the initial drop, the values were constantly decreasing as a consequence of both stress relaxation and damage accumulation. In this region, the reduction extent depended on the loading conditions, but it was generally limited. Ultimately, a final major drop was present, which was related to the macroscopic crack growth and consequence of the reduction of the specimens' ligament section. As depicted in Figure 8a, the maximum forces decreased at the very beginning of the tests by about one third and stabilized up to the final failure of the specimen, regardless of the strain level. Conversely, a different behavior was found for the minimum force. It is evident from the plot in Figure 8b that after the initial drop, the minimum forces are different for distinct strain levels: larger drops were monitored for lower values of maximum strain. Furthermore, small peaks were observed after the initial drop, which are more evident for larger maximum strains in particular. Fatigue data were further analyzed by considering the peak values of the forces detected during the measurements. For a clearer comparison at different loading conditions, the maximum force, F max , and the minimum force, F min , were normalized with respect to their respective initial values and they were plotted against the normalized number of cycles (Figure 8a,b). Due to stress relaxation, the maximum and minimum forces showed a considerable variation in the first cycles, until consistent stress relaxation occurs in the specimen. After the initial drop, the values were constantly decreasing as a consequence of both stress relaxation and damage accumulation. In this region, the reduction extent depended on the loading conditions, but it was generally limited. Ultimately, a final major drop was present, which was related to the macroscopic crack growth and consequence of the reduction of the specimens' ligament section. As depicted in Figure 8a, the maximum forces decreased at the very beginning of the tests by about one third and stabilized up to the final failure of the specimen, regardless of the strain level. Conversely, a different behavior was found for the minimum force. It is evident from the plot in Figure 8b that after the initial drop, the minimum forces are different for distinct strain levels: larger drops were monitored for lower values of maximum strain. Furthermore, small peaks were observed after the initial drop, which are more evident for larger maximum strains in particular.
Therefore, several parameters showed a reduction during the cycles, which was related to relaxation and material damage. These changes could be observed from the early stages of fatigue, suggesting a slow crack growth process; however, clear evidence to identify a transition between the regimes of crack initiation and crack growth were not found. Thus, further investigations were performed considering a fracture mechanics-based approach. Therefore, several parameters showed a reduction during the cycles, which was related to relaxation and material damage. These changes could be observed from the early stages of fatigue, suggesting a slow crack growth process; however, clear evidence to identify a transition between the regimes of crack initiation and crack growth were not found. Thus, further investigations were performed considering a fracture mechanics-based approach.

Fracture Surface Analysis of the Unnotched Dumbbell Specimens
Subsequently, the fracture surfaces of the specimens were observed under the light microscope. The typical appearance of the fracture surfaces is shown and described in Figure 9. The most evident characteristic is the presence of fatigue striations which display an increasing distance from each other as they move farther from the initiation area. Normally, fatigue striations were observed for up to half of the surface area, whereas the remaining section corresponded to the catastrophic failure. An overview of the fracture surfaces at different strains is reported in Figure 10. Generally, failure occurred perpendicularly to the loading direction, with quite flat fracture surfaces, and the features displayed in Figure 9 were the most frequently observed with evident fatigue striation in all cases; however, in some cases, the fracture surfaces results were more complex, with evident macroscopic cracks out of the plane of the fracture surface at εmax of 50% and 60%. Such kinds of crack deviation could be related to local inhomogeneities in the material such as large agglomerates of carbon black. Moreover, at the lowest strain (εmax of 50%), fatigue striations in two opposite directions were observed (Figure 11), suggesting that the initiation is taking place independently in different areas and that at lower strains, macroscopic crack growth took place in different directions. The single crack direction at larger strains suggests that due to the larger energy and stress intensification involved, once the crack began to form, it became the predominant crack.

Fracture Surface Analysis of the Unnotched Dumbbell Specimens
Subsequently, the fracture surfaces of the specimens were observed under the light microscope. The typical appearance of the fracture surfaces is shown and described in Figure 9. The most evident characteristic is the presence of fatigue striations which display an increasing distance from each other as they move farther from the initiation area. Normally, fatigue striations were observed for up to half of the surface area, whereas the remaining section corresponded to the catastrophic failure. An overview of the fracture surfaces at different strains is reported in Figure 10. Generally, failure occurred perpendicularly to the loading direction, with quite flat fracture surfaces, and the features displayed in Figure 9 were the most frequently observed with evident fatigue striation in all cases; however, in some cases, the fracture surfaces results were more complex, with evident macroscopic cracks out of the plane of the fracture surface at ε max of 50% and 60%. Such kinds of crack deviation could be related to local inhomogeneities in the material such as large agglomerates of carbon black. Moreover, at the lowest strain (ε max of 50%), fatigue striations in two opposite directions were observed (Figure 11), suggesting that the initiation is taking place independently in different areas and that at lower strains, macroscopic crack growth took place in different directions. The single crack direction at larger strains suggests that due to the larger energy and stress intensification involved, once the crack began to form, it became the predominant crack.

Fatigue of the Cracked Round Bar Specimens
Furthermore, fatigue measurements on the Crack Round Bar (CRB) were conducted. The CRB geometry (a dumbbell with a circumferential notch in the middle) was developed for accelerated testing of quasi-brittle failure resistance of polymers, elastomers, and thermoplastic elastomers [68][69][70]. In principle, due to the presence of the notch, this specimen geometry would allow the crack length and crack growth rate to be measured using a Paris approach, focusing on the crack propagation regime; however, due to the cylindrical geometry, the crack front results are complex, and it is difficult to assess crack length measurements through optical methods. More accurate results could be obtained by using extensometers and retrieving the crack length through a calibration between Crack Opening Displacement (COD) and the actual crack length [71]; however, due to the large strain involved, this approach has limitations for elastomers.
The results of fatigue tests obtained on CRB specimens in comparison with the results of standard specimens are shown in Figure 12. As depicted, a general reduction in lifetime was found for CRB as a consequence of the introduction of the notch. A reduction of about four orders of magnitude of the fatigue lifetime at εmax of 50% was recorded. Furthermore, in contrast to the results of unnotched axisymmetric dumbbells, failure was observed even at 20% of maximum (nominal) strain. The fitting power law from the statistical analysis of the results, performed through the ASTM E739-10 was found to be: Figure 11. Light microscope picture of the fracture surface of an axisymmetric dumbbell specimen of CB filled NBR after fatigue testing in which two fronts of macroscopic fatigue crack growth are present (displayed by the arrows). The test was performed at a frequency of 4 Hz and a load ratio Rε of 0.5.

Fatigue of the Cracked Round Bar Specimens
Furthermore, fatigue measurements on the Crack Round Bar (CRB) were conducted. The CRB geometry (a dumbbell with a circumferential notch in the middle) was developed for accelerated testing of quasi-brittle failure resistance of polymers, elastomers, and thermoplastic elastomers [68][69][70]. In principle, due to the presence of the notch, this specimen geometry would allow the crack length and crack growth rate to be measured using a Paris approach, focusing on the crack propagation regime; however, due to the cylindrical geometry, the crack front results are complex, and it is difficult to assess crack length measurements through optical methods. More accurate results could be obtained by using extensometers and retrieving the crack length through a calibration between Crack Opening Displacement (COD) and the actual crack length [71]; however, due to the large strain involved, this approach has limitations for elastomers.
The results of fatigue tests obtained on CRB specimens in comparison with the results of standard specimens are shown in Figure 12. As depicted, a general reduction in lifetime was found for CRB as a consequence of the introduction of the notch. A reduction of about four orders of magnitude of the fatigue lifetime at ε max of 50% was recorded. Furthermore, in contrast to the results of unnotched axisymmetric dumbbells, failure was observed even at 20% of maximum (nominal) strain. The fitting power law from the statistical analysis of the results, performed through the ASTM E739-10 was found to be: Moreover, the variations of both the normalized maximum force and crack length during the test with CRB geometry were analyzed and reported in Figure 13a. It is worth noting that the machine reached such high strains with a delay; since a stable deformation was observed after a transient period of eight cycles, the maximum force was attained after nine cycles. Considering the cylindrical symmetry of the samples, the crack length was approximated as the projection of the actual radial crack. In addition, the reported crack length was considered as the sum of the crack lengths observed on both sides of the specimen in the camera record (see the picture in Figure 13b). To have a proper measurement of the crack evolution, more cameras would be needed; however, this approximation is considered sufficiently accurate to monitor the evolution of the specimen's ligament section. As depicted in Figure 13a, the larger crack increase was observed at about 30 cycles, a time at which an abrupt variation of the maximum force was also recorded; nevertheless, a smaller crack length increase was present even before this stage. In actuality, the crack starts to grow immediately after the application of the load-at the ninth cycle after stable deformation is reached. On the other hand, in unnotched specimens, the microscopic crack growth occurred after initiation. Moreover, the variations of both the normalized maximum force and crack length during the test with CRB geometry were analyzed and reported in Figure 13a. It is worth noting that the machine reached such high strains with a delay; since a stable deformation was observed after a transient period of eight cycles, the maximum force was attained after nine cycles. Considering the cylindrical symmetry of the samples, the crack length was approximated as the projection of the actual radial crack. In addition, the reported crack length was considered as the sum of the crack lengths observed on both sides of the specimen in the camera record (see the picture in Figure 13b). To have a proper measurement of the crack evolution, more cameras would be needed; however, this approximation is considered sufficiently accurate to monitor the evolution of the specimen's ligament section. As depicted in Figure 13a, the larger crack increase was observed at about 30 cycles, a time at which an abrupt variation of the maximum force was also recorded; nevertheless, a smaller crack length increase was present even before this stage. In actuality, the crack starts to grow immediately after the application of the load-at the ninth cycle after stable deformation is reached. On the other hand, in unnotched specimens, the microscopic crack growth occurred after initiation.

Fracture Surface Analysis of the Cracked Round Bar Specimens
Additionally, the fracture surfaces of the CRB specimen were analyzed and a representative example is reported in Figure 14. There are several similarities in the failure mechanism when compared with the unnotched kind in Figure 9. In addition to the region of radial notching, the fracture surfaces of CRBs also displayed the three aforementioned regions of initiation, fatigue striation, and catastrophic failure. Moreover, even in this case, the fatigue striations show increasing distances moving from the initiation region, and in correlation with the increasing crack length reported in Figure 13a. Overall, these results suggest that a similar damage mechanism is maintained, and that the introduction of the artificial notch is equivalent to a specimen with a larger initial defect in an unnotched specimen. Figure 13. (a) Normalized maximum force (circles) and crack length (diamonds) against number of cycles for CRB specimens obtained from axisymmetric dumbbell specimens of CB filled NBR. (b) Picture of the tested CRB at N = 30; on both sides, the crack opening can be recognized. The test was performed at the maximum strain of 50% at a frequency of 4 Hz and a load ratio R ε of 0.5.

Fracture Surface Analysis of the Cracked Round Bar Specimens
Additionally, the fracture surfaces of the CRB specimen were analyzed and a representative example is reported in Figure 14. There are several similarities in the failure mechanism when compared with the unnotched kind in Figure 9. In addition to the region of radial notching, the fracture surfaces of CRBs also displayed the three aforementioned regions of initiation, fatigue striation, and catastrophic failure. Moreover, even in this case, the fatigue striations show increasing distances moving from the initiation region, and in correlation with the increasing crack length reported in Figure 13a. Overall, these results suggest that a similar damage mechanism is maintained, and that the introduction of the artificial notch is equivalent to a specimen with a larger initial defect in an unnotched specimen.
Additionally, the fracture surfaces of the CRB specimen were analyzed and a representative example is reported in Figure 14. There are several similarities in the failure mechanism when compared with the unnotched kind in Figure 9. In addition to the region of radial notching, the fracture surfaces of CRBs also displayed the three aforementioned regions of initiation, fatigue striation, and catastrophic failure. Moreover, even in this case, the fatigue striations show increasing distances moving from the initiation region, and in correlation with the increasing crack length reported in Figure 13a. Overall, these results suggest that a similar damage mechanism is maintained, and that the introduction of the artificial notch is equivalent to a specimen with a larger initial defect in an unnotched specimen. Figure 14. Light microscope picture of the fracture surface of a notched axisymmetric dumbbell specimen (CRB) of CB filled NBR after fatigue testing. At the specimen edge the radial notch can be observed, and a further three regions are displayed: initiation region (yellow), fatigue striation Figure 14. Light microscope picture of the fracture surface of a notched axisymmetric dumbbell specimen (CRB) of CB filled NBR after fatigue testing. At the specimen edge the radial notch can be observed, and a further three regions are displayed: initiation region (yellow), fatigue striation (orange)-corresponding to macroscopic fatigue crack growth-and catastrophic failure (red). The test was performed at a frequency of 4 Hz and a load ratio R ε 0.5.

Defect Size Assesment
As already mentioned, the fatigue lifetime can be correlated with the crack growth by considering the initial size of defects that cause fatigue failure. In fact, by considering the defect size c 0 , it is possible to connect the failure of the material to defects and inhomogeneities naturally present in the material.
In order to obtain the value c 0 , the integration of the crack growth characteristic is required. Generally, the fatigue crack growth of rubbers within the stable regime can be described using the Paris law [5]: where dc/dN represents the crack growth rate, G represents the tearing energy, and C and m are materials constants. For small cracks, the energy release rate (i.e., tearing energy) can be factored into strain energy density and crack size [2] and its estimation can be given by: where W is the strain energy density and c is the crack length; k(λ) is a function of the stretch ratio [72]. In actuality, Equation (7) represents the tearing energy for single edge notch specimens [3]; nevertheless, it is a good approximation for small cracks in terms of initial defects. Moreover, the strain energy density W was evaluated considering the area under the loading curve by numerical integration of hysteresis data as: where σ is the engineering stress, ε is the engineering strain, F is the force, A the resistant area of the specimen, s is the displacement, s min and s max are the minimum and the maximum displacements, respectively, and V is the volume of the specimen. By combining Equations (6) and (7), and integrating the number of cycles for failure, N f , can be obtained: Assuming that the initial defects size c 0 is much smaller than the final crack length c f , this last contribution can be neglected, obtaining: This equation can be used to calculate the number of cycles to failure as a consequence of the growth of a pre-existing defect with the dimension c 0 . Alternatively, Equation (10) can be reversed to evaluate the effective critical defect size from which the failure originated. It is worth noting that for the evaluation of the overall lifetime, this approach for elastomers often neglects the initiation time; however, this is justified considering what was found through the analysis of the evolution of defects during fatigue, which investigated for filled polychloroprene rubber (CR) through µ-CT [31,32], and for filled natural rubber (NR) through SEM analysis [28,29]. It was shown that cracks are initiated in the early stages of fatigue tests and that the number of defects mainly increase within the first 10% of fatigue life. In the successive stages, the initiated defects slowly increased in size, whereas the number of cracks is more stable: it constantly increases but with a lower rate of nucleation. Therefore, the majority of the lifetime of rubbers is governed by crack growth and the initiation time can be neglected.
The effective critical defect size c 0 was calculated by reversing Equation (10) and using the number of cycles to failure N f , as reported in Figure 3. The values of k and W were calculated for the corresponding strain levels. In particular, the strain energy density for each specimen was evaluated through Equation (8) at the beginning of the test. As for the coefficients m and C of the Paris law, these were evaluated from the data and the fatigue master curve shifting procedure was previously reported in [35]. In that work, fatigue crack growth characteristics at different temperatures (and at the same frequency and load ratio as in the work presented here) were analyzed, and a fatigue master curve based on the temperature dependence of the loss modulus was constructed. Using this procedure, horizontal shift factors were evaluated for the temperatures measured during the testing of axisymmetric dumbbells, allowing for the calculation of m and C for each test, while accounting for the influence of temperature on these parameters. An overview of all the values used for the calculation of c 0 are reported in Table 1. Using these parameters, the effective critical defect sizes were calculated, and the values at different strains are reported in Figure 15. As depicted, the values c 0 were found to be in the range between 4 and 14 µm for the different maximum strains. The values of the critical defect size are similar, regardless of the applied strain level: this points towards the conclusion that the defect size is an intrinsic value related to the material. The reason for the different degrees of dispersion in the results at different strains was mainly related to the large scattering of the fatigue data at the same strain values (see Table 1). Due to the independence of the size on the strain, it is possible to consider the mean value for all specimens: the effective critical defect size c 0 was 9 ± 3 µm. Table 1. Overview of the values used to calculate the critical defect size with the reversion of Equation (10). The values of the number of cycles to failure are the same as reported in Figure 3. each specimen was evaluated through Equation (8) at the beginning of the test. As for the coefficients m and C of the Paris law, these were evaluated from the data and the fatigue master curve shifting procedure was previously reported in [35]. In that work, fatigue crack growth characteristics at different temperatures (and at the same frequency and load ratio as in the work presented here) were analyzed, and a fatigue master curve based on the temperature dependence of the loss modulus was constructed. Using this procedure, horizontal shift factors were evaluated for the temperatures measured during the testing of axisymmetric dumbbells, allowing for the calculation of m and C for each test, while accounting for the influence of temperature on these parameters. An overview of all the values used for the calculation of c0 are reported in Table 1. Using these parameters, the effective critical defect sizes were calculated, and the values at different strains are reported in Figure 15. As depicted, the values c0 were found to be in the range between 4 and 14 μm for the different maximum strains. The values of the critical defect size are similar, regardless of the applied strain level: this points towards the conclusion that the defect size is an intrinsic value related to the material. The reason for the different degrees of dispersion in the results at different strains was mainly related to the large scattering of the fatigue data at the same strain values (see Table 1). Due to the independence of the size on the strain, it is possible to consider the mean value for all specimens: the effective critical defect size c0 was 9 ± 3 μm. Figure 15. Critical defect size c0 at different strains evaluated from fatigue data using Equation (10). The fatigue tests were implemented on axisymmetric dumbbell specimens of CB filled NBR at a frequency of 4 Hz and a load ratio Rε of 0.5. Table 1. Overview of the values used to calculate the critical defect size with the reversion of Equation (10). The values of the number of cycles to failure are the same as reported in Figure 3.  Figure 15. Critical defect size c 0 at different strains evaluated from fatigue data using Equation (10). The fatigue tests were implemented on axisymmetric dumbbell specimens of CB filled NBR at a frequency of 4 Hz and a load ratio R ε of 0.5.
Moreover, the initial defect size of the axisymmetric dumbbell specimens was evaluated through the reconstruction of µ-CT (Figure 16a). From the reconstructed volume, naturally present defects have been identified (Figure 16b) and their volumes have been determined. In order to have a direct comparison with the critical defect size c 0 , the cavities have been approximated as spheres and their radii have been calculated. The particle size distribution is displayed in Figure 16c, in which the relative frequency for each value of the particle radius is reported. As depicted, the dimensions of detected cavities were of few µm, between 1 and 24 µm. It is important to stress that considering the temperature influence on crack growth parameters strongly affects the calculation of the defects' dimensions; however, a very good agreement between CT data and the back-calculation from CT measurements is obtained. Indeed, the defect size distribution indeed vanishes above 8-10 µm, which corresponds quite well with the value of c 0 of 9 µm previously identified (see inset of Figure 16c). the particle radius is reported. As depicted, the dimensions of detected cavities were of few μm, between 1 and 24 μm. It is important to stress that considering the temperature influence on crack growth parameters strongly affects the calculation of the defects' dimensions; however, a very good agreement between CT data and the back-calculation from CT measurements is obtained. Indeed, the defect size distribution indeed vanishes above 8-10 μm, which corresponds quite well with the value of c0 of 9 μm previously identified (see inset of Figure 16c). In addition, further observations could be made by calculating the lifetime based on the evaluated defects. In fact, Equation (10) was used to estimate the lifetime for the two interrupted tests. More specifically, the average c0 (9 μm) was used as the initial defect size, and the calculated fatigue life for the maximum strains of 5% and 20% were of 3.5·× 10 12 and 2.9·× 10 8 , respectively. The results are well above the number of cycles at which the tests were stopped (7·× 10 6 and 2·× 10 6 , respectively). The details concerning the parameters used for the calculation are reported in Table 2. In addition, further observations could be made by calculating the lifetime based on the evaluated defects. In fact, Equation (10) was used to estimate the lifetime for the two interrupted tests. More specifically, the average c 0 (9 µm) was used as the initial defect size, and the calculated fatigue life for the maximum strains of 5% and 20% were of 3.5 × 10 12 and 2.9 × 10 8 , respectively. The results are well above the number of cycles at which the tests were stopped (7 × 10 6 and 2 × 10 6 , respectively). The details concerning the parameters used for the calculation are reported in Table 2. Similarly, fatigue lifetime predictions were calculated considering the mean values of c 0 and the minimum and maximum values of particle radius found by µ-CT (so 1.7 and 22.2 µm, respectively). Some approximations have been made in the calculations. In particular, the linear relation between the strain energy density W and maximum strain (reasonable in the strain range between 50 and 65%) was considered, which resulted in the following equation: The values of k were calculated at the different strains, whereas the material coefficients m and C were estimated through the aforementioned temperature shifting procedure. A linear relation between temperature and maximum strain was considered: Therefore, the fatigue lifetime predictions were calculated using Equation (10), and these results were then compared with the experimental fatigue results described in the previous section. These results illustrate the impact of initial defect size. As depicted in Figure 17, when considering the minimum and maximum value of particle size, the fatigue predictions show a broad spectrum spanning across three orders of magnitude. Of course, the calculations using the average c 0 give values similar to the experimental data and to the power law fit (Equation (1)): these have been included to illustrate the typical scatter in the data.
The values of k were calculated at the different strains, whereas the material coefficients m and C were estimated through the aforementioned temperature shifting procedure. A linear relation between temperature and maximum strain was considered: Therefore, the fatigue lifetime predictions were calculated using Equation (10), and these results were then compared with the experimental fatigue results described in the previous section. These results illustrate the impact of initial defect size. As depicted in Figure 17, when considering the minimum and maximum value of particle size, the fatigue predictions show a broad spectrum spanning across three orders of magnitude. Of course, the calculations using the average c0 give values similar to the experimental data and to the power law fit (Equation (1)): these have been included to illustrate the typical scatter in the data. Figure 17. Comparison between fatigue data and the calculated fatigue lifetimes for the calculated defect size (c0,calc), and the maximum and minimum sizes evaluated from μ-CT (c0,max and c0,min, respectively).

Correlation of the Results Based on J-Integral
Finally, a non-linear J-integral hyperelastic model for the CRB has been developed. This energetic approach has the advantage of being applied in order to describe materials that show non-linear mechanical behavior, and it is useful for comparing different geometries. In particular, the J-integral was evaluated as [73]: Figure 17. Comparison between fatigue data and the calculated fatigue lifetimes for the calculated defect size (c 0,calc ), and the maximum and minimum sizes evaluated from µ-CT (c 0,max and c 0,min , respectively).

Correlation of the Results Based on J-Integral
Finally, a non-linear J-integral hyperelastic model for the CRB has been developed. This energetic approach has the advantage of being applied in order to describe materials that show non-linear mechanical behavior, and it is useful for comparing different geometries. In particular, the J-integral was evaluated as [73]: where F is the force, r out the external radius of the CRB, a is the size of notch in the CRB, and f (a/r out ) is a geometric factor found by the FEM simulation as: Equation (13) was used to compute J max and J min using F max and F min , respectively. The difference between J max and J min was then used to compare the CRB results with those of standard dumbbell samples. The J-integral of the axisymmetric dumbbell was evaluated using Equation (7); in fact, J can be approximated as the energy release rate G [60]. The values of J for dumbbell specimens were evaluated using the calculated c 0 and the strain energies reported in Table 1. The considered forces were those measured at the beginning of the test. An overview of the data used for the evaluation of J is reported in Table 3. The J-integral as a function of number of cycles to failure for both geometries is plotted in Figure 18. A unique fitting curve was found for both geometries; this demonstrates that independently of the specimen geometry, it is possible to evaluate the fatigue lifetime based on J-integral. Although this J-integral formulation and its applicability still need to be verified on different geometric parameters and materials, it seems to be a promising candidate for extending fracture mechanic tools for the fatigue assessment of elastomers with reduced testing times. energies reported in Table 1. The considered forces were those measured at the beginning of the test. An overview of the data used for the evaluation of J is reported in Table 3. The J-integral as a function of number of cycles to failure for both geometries is plotted in Figure 18. A unique fitting curve was found for both geometries; this demonstrates that independently of the specimen geometry, it is possible to evaluate the fatigue lifetime based on J-integral. Although this J-integral formulation and its applicability still need to be verified on different geometric parameters and materials, it seems to be a promising candidate for extending fracture mechanic tools for the fatigue assessment of elastomers with reduced testing times.  Figure 18. Wöhler curve in terms of the J-integral for axisymmetric dumbbell (blue) and CRB (red) specimens. A unique fitting line was found correlating the two geometries.

Conclusions
The fatigue behavior of a carbon black filled NBR was investigated in the presented work. Fatigue measurements were performed on axisymmetric dumbbell specimens and also on CRB specimens (i.e., axisymmetric dumbbells with a circumferential notch). The number of cycles to failure was measured for both sets, and Wöhler curves were constructed. Although the recorded points were considerably scattered, the results could be described with a power law.
The results of the axisymmetric dumbbell specimens were investigated in more detail through the analysis of temperature variations, hysteresis, secant, dynamic moduli, stored and dissipated energies, and the normalized force peaks. In general, all the observed values largely varied in the early stages of fatigue (approximately up to 2500-3000 cycles), but they stabilized with a lower degree of variation up to the failure of the specimens, suggesting a slow crack growth process. Even though these variations evidenced a continuous and progressive accumulation of relaxation and damage, a clear distinction between initiation and propagation was not observed.
Additional details could be obtained by further investigations considering a fracture mechanics approach. The fracture surfaces featured initiation areas from which evident fatigue striations originated from. Similar aspects of the fracture surfaces were found for axisymmetric dumbbells with a circumferential notch. These results suggested similarities in the damaging mechanism and that the notch can be represented as a larger defect. As expected, the introduction of a notch led to a reduction in fatigue life up to four orders of magnitude.
Moreover, the fatigue life of dumbbells was correlated to the defect size by integrating the crack growth characteristic and taking into account the effect of temperature, obtaining values that were independent on the strain level. The evaluated defect size was then compared with the defect sizes analyzed through X-ray microtomography, confirming similar dimensions and a good accuracy of the calculations.
Finally, through the evaluation of the J-integral, it was possible to link the fatigue life of unnotched and notched axisymmetric dumbbells, evidencing a geometry independency of the J-integral formulation. This formulation for notched axisymmetric dumbbells may represent a new tool for the investigation of accelerated fatigue in elastomers.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ma15113745/s1, Figure S1: Comparison of Wöhler curves considering N i and N f .

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.