Peridynamic numerical investigation of asymmetric strain-controlled fatigue behaviour using the kinetic theory of fracture

Numerical fatigue process modelling is complex and still an open task. Discontinuity caused by fatigue cracks requires special finite element techniques based on additional parameters, the selection of which has a strong effect on the simulation results. Moreover, the calculation of fatigue life according to empirical material coefficients (e.g., Paris law) does not explain the process, and coefficients should be set from experimental testing, which is not always possible. A new nonlocal continuum mechanics formulation without spatial derivative of coordinates, namely, peridynamics (PD), which was created 20 y ago, provides new opportunities for modelling discontinuities, such as fatigue cracks. The fatigue process can be better described by using the atomistic approach-based kinetic theory of fracture (KTF), which includes the process temperature, maximum and minimum stresses, and loading frequency in its differential fatigue damage equation. Standard 316L stainless steel specimens are tested, and then the KTF-PD fatigue simulation is run in this study. In-house MATLAB code, calibrated from the material S ‒ N curve, is used for the KTF-PD simulation. A novel KTF equation based on the cycle stress ‒ strain hysteresis loop is proposed and applied to predict fatigue life. The simulation results are compared with the experimental results, and good agreement is observed for both symmetric and asymmetric cyclic loading.


Introduction
Fatigue life prediction is very important and is still one of the most challenging engineering tasks.For example, the fatigue life of 316L stainless steel constructions used in reactors from nuclear power plants is estimated to be 60 years at certain cyclic loadings.Newer studies suggest that 60 years of exploitation time is too conservative; the same structures can last up to 80 years, and more materials, energy and manufacturing time can be saved.To predict the fatigue life of such responsible structures more precisely, reliable fatigue damage detection methods and numerical fatigue models validated against experiments are necessary.Numerical fatigue process modelling faces several issues.Most of them are related to the finite element (FE) formulation itself.The FE governing equations of motion are cast in terms of spatial derivatives of displacement, which is not valid for discontinuities, e.g., fatigue cracks.Although special FE methods, such as the extended finite element method (XFEM) [1], virtual crack closure technique (VCCT) [2], cohesive zone model (CZM) [3] or even a combination of these methods [4], can be applied to address the discontinuity problem caused by fatigue cracking, the complexity of the numerical fatigue model and possible uncertainty increase, especially for 3D cracking.In this aforementioned problem, the fatigue crack in the FE model cannot be simulated directly, e.g., without previously defining the crack path.The crack path definition has a strong effect on the simulation results; for example, if true crack paths are unknown, then numerical fatigue failure prediction may be inaccurate.The second issue of numerical fatigue modelling comes from the fatigue damage definition based on empirical material coefficients [5] (e.g., Paris's law [6] for crack growth and Miner's rule for fatigue damage summation).Despite the fact that the numerical model can be calibrated and coefficients can be chosen according to the material S-N curve and crack growth rate, for other testing conditions (different stress ratios and environmental temperatures), the material coefficient values are different and cannot be calculated without experimental testing under those conditions.Although the Goodman and Gerber relations [7] can be helpful in evaluating fatigue life at different stress ratios, for more precise fatigue life prediction, experimental validation is still needed.
The alternative theory of classical continuum mechanics (CCM), that is, peridynamics (PD), reformulates CCM equations in terms of spatial integrals (without spatial derivatives), and this approach is very useful for the case of (fatigue) crack modelling [8].PD was originally created by Silling [9] in 2000 and applied only for modelling brittle materials, but later PD models were further developed for both brittle and ductile materials by Silling et al. [10], Madenci and Oterkus [11], Bobaru and Hu [12] and others, as described in detail in the PD review given by Javili and Morasata [13].Currently, PD models, which can be grouped into bond-based PD (BBPD, the simplest ones with a fixed material Poisson's ratio) [9,14] and state-based PD (SBPD, free of BBPD limitations and more advanced ones that allow complex material deformations) [15,16] models are used to model various mechanical (static and dynamic crack propagation, impact damage, wave propagation) [17]) and physical (thermal, diffusion, etc.) [18,19].However, CCM reformulation in PD theory creates new BBPD and SBPD material stiffness and strength parameters [20], but their values can be derived from the deformation energy equilibriums of the CCM and PD theories, as opposed to the CZM, XFEM, and VCCT methods, where FE mesh and crack parameters have to be adjusted by the user to ensure model correlation to experiment, stability, etc. Silling and Askari [21], Oterkus et al. [22] and others [23][24][25] successfully applied the PD theory for the fatigue modelling of brittle and ductile materials using cyclic failure criteria based on the Paris law [6].Silling and Askari [21] demonstrated that the single PD model can predict the initiation of fatigue cracks, crack growth, and final failure phases without predefined crack paths.Pashazad and Khrazi [15] showed that even kinematic, isotropic, and mixed hardening effects can be captured in the PD model.The Freimanis and Kaewunruen study of fatigue damage in rail squats [25], Zhang fatigue analysis of isotropic and polycrystalline materials [26] and other works [17] demonstrated excellent PD simulation agreement with experimental testing results (correlation coefficients R 2 [17] between the PD simulation and the experimental curves up to 0.97 ÷ 1.00).All abovementioned PD fatigue studies, including recent PD applications for fatigue crack growth in hydrogels [27], are still based on empirical fatigue process coefficients, which can be found only from experiments and are not valid for different testing conditions.
The issue of purely empirical cyclic material empirical coefficients can be addressed by taking the kinetic theory of fracture (KTF) instead of the Paris law approach.KTF was created separately by Coleman [28] and Zhurkov [29] in the middle of the XX century and treats fatigue damage as a thermally activated process that is dependent on the maximum and minimum cycle stresses, process temperature, and cyclic loading frequency.Although the KTF fatigue model parameters still must be calibrated based on the material S-N curve, compared with the Paris law, KTF equations better describe the fatigue process; thus, the same model can be easily implemented to simulate the process under different conditions (stress ratio, temperature, and loading frequency).KTF-based fatigue life predictions [29,30] are in agreement with experimental results for various materials, such as composites and metals, over a wide range of temperatures and stress ranges.
Madenci et al. first applied KTF combined with PD [5] and further developed it for modelling constant amplitude aluminium [14] and constant, variable amplitude and stress ratio composite cyclic loading [31].Although KTF models correlate well with experimental results, for asymmetric cyclic loading at low stress amplitudes, the KTF fatigue damage equation does not provide the correct fatigue life [32].To ensure agreement with the experimental data, Fertig and Kenik [32,33] modified this equation for composite materials, including energy losses proportional to the square of the cyclic loading stress amplitude.To date, no works that continue to combine the research of Madenci et al. [5,14,31] and Fertig and Kenik [32,33] have been conducted; therefore, the development and establishment of KTF-PD models for symmetric loading is an important topic, especially for asymmetric cyclic loading.The aim of this study is to develop a KTF-PD model for fatigue life prediction in metallic materials based only on symmetric cyclic loading S-N data and to apply this model to predict fatigue life under symmetric and asymmetric cyclic loading.The prepared model is validated using experimental data.Standard cylindrical 316L stainless steel specimens are tested, and then a KTF-PD fatigue simulation is carried out in this study.In-house PD MATLAB code, calibrated from the material S-N curve, given in the datasheets [34], is used for the KTF-PD simulation.A novel KTF equation that extends energy loss based on a cyclic stress-strain hysteresis loop is proposed and applied to predict fatigue life under asymmetric cyclic loading.The simulation results are compared with the experimental results, and good agreement is observed.The manuscript is organized as follows: PD and KTF theories and their numerical implementation details are provided in "Theoretical background"; the experimental testing and KTF-PD fatigue modelling in MATLAB details are presented in "Experiments and modelling"; the fatigue simulation results and their comparison to the experimental results are discussed in "Results and discussion"; and the summary and conclusions are provided in "Conclusion".

Peridynamics theory
Peridynamics (PD) theory due to its integral equations of motion does not contain spatial derivatives.Thus, the PD approach is different than the CCM equations, which are expressed in terms of the partial derivatives of displacement.According to PD theory, it is assumed that each material point x interacts with the material points x' in the interaction range H x (Fig. 1), which is usually a sphere with radius δ ≈ 3Δx in the 3D PD model.The PD material domain Ω consists of a grid of equally spaced material points, where Δx is the PD grid spacing.The interaction between points x and x' is defined as the PD bond.
Then, the equation of motion in the simplest case of PD theory (bond-based (BBPD)) can be written as [35]: where x, u and x', u' are position and displacement vectors (Fig. 1) associated with the material points x and x', respectively; b(x, t) represents the body forces; f represents the PD bond force, ρ is the density of the material and V x is the volume of the material point x.Conventional PD is a macroscale model because the number of PD material points in the macroscopic body numerical model is finite.In the BBPD model (Fig. 1), the PD bond forces f and f' are equal in magnitude and parallel in direction, resulting in a Poisson's ratio of only 0.25 for the 3D and 2D plane strain BBPD models and 0.33 for the 2D plane stress model [13].Additionally, plastic deformations in BBPD can be modelled only as the permanent deformation of a material undergoing a volumetric strain without shear [10].To circumvent these previously mentioned BBPD limitations, more advanced SBPD models have been developed, namely, ordinary and non-ordinary state-based (1)  PD (OSBPD and NSBPD) models [13].However, the computational time of the SBPD models significantly increases compared to that of BBPD.Based on the computational resource limitations and acceptable accuracy, the BBPD model is used in this study.Moreover, the complexity of SBPD models renders their implementation in commercial FE software not possible for existing FE formulations [36].In contrast, the BBPD fatigue model has been successfully implemented in the MATRIX 27 elements in the ANSYS framework by Zhang and Madenci [14].Therefore, there is incentive to use and develop the BBPD fatigue model to increase the availability of the PD method.
The additional difficulty in both the BBPD and SBPD models is related to the reformulation of the CCM, which results in new PD material quantities, such as PD bond strains (stretches).Bond, connecting points x and x', stretch is calculated as follows: Similarly, for the stress-strain relation in CCM, the PD bond forces f (also see Eq. ( 1)) are related to the PD bond stretch (Fig. 2a), and for the BBPD formulation, the bond forces are computed as follows: where c is the PD bond stiffness for the 3D BBPD elastic material model expressed as follows: where E is the elastic modulus and δ is the PD horizon size (Fig. 1).
The PD bond force-stretch relation, defined in Eqs.(3)(4) and presented in Fig. 2a, is not valid for the elasto-plastic material model and thus should be updated by changing (2) ( the PD bond stiffness c after the material proportional limit is reached.Similar to the CCM tensile diagram, bilinear [37] or power law [38] approximation, the same type of force-stretch relation (Fig. 2b) can be used.The proportional or yield limit stretch value is derived by equalizing the CCM and PD strain energy densities at the proportional or yield limit deformation [11]: where σ ij and ε ij are the CCM stress and strain tensor components, respectively.Special calculations are required to find the CCM strain and stress tensor of the material PD point x.The deformation gradient F x of the PD point x, which is valid even for a discontinuous displacement field, is also calculated without spatial derivatives of the displacement by using the peridynamic differential operator [10,39] and expressed as follows: The term ω(x ′ − x) in Eq. ( 6) is the PD bond length-dependent influence function, the selection of which was analysed in several papers [40,41], and it was found that the linear function for the 3D model is sufficient.Once the strain gradient of point x is computed by the PDDO, further CCM equations are applied to obtain the Lagrange strain tensor E of point x: Then, the stress tensor of point x components can be found according to the 3D generalized Hooke's law relation: where C ijkl are components of the CCM material stiffness tensor and e kl are compo- nents of the Lagrange strain tensor from Eq. ( 7).When plasticity occurs, the material stiffness tensor components should be updated to follow the elastoplastic behaviour of the material.
Static failure of the PD bond is considered when its stretch exceeds the defined critical value s c .The critical stretch cannot be directly related to the CCM failure strain or stress because the stretch of the PD bond depends on the PD interaction range (horizon) size δ and is thus not equal to the strain used in the CCM.The critical stretch s c value used in this study is calculated according to the relationship between PD s c and the material critical energy release rate G c , derived by Madenci and Oterkus [35].For the 3D BBPD case, it is calculated as follows: (5 C ijkl e kl i, j = 1 . . .3, .Fatigue modelling under PD still consists of static simulations, but additional cyclic PD bond failure criteria are needed.Although the critical stretch s c value of the PD bond can be reduced with an increasing number of cycles that result in cyclic failure, such a model [22] cannot capture the fatigue initiation process.The PD bond "remaining life" variable q approach, proposed by Silling and Askari [21], is better because it allows modelling of both the crack initiation and growth phases.q depends on the loading of the PD bond and the number of cycles, when its initial value q(x, x'-x, 0) = 1, and the cyclic failure of the bond is considered when q(x, x'-x, N) ≤ 0.Then, both the static and cyclic PD bond failures in the PD fatigue model can be evaluated by multiplying the PD bond stiffness c by the historydependent scalar value function µ(x'-x), expressed as follows: Finally, the PD damage at point x as a ratio of broken PD bonds N µ(x ′ −x)=0 to the number N of total PD bonds at point x horizon H x is defined as follows:

Kinetic theory of fracture
The evolution of fatigue damage n (n varies from 0 to 1; n = 0 indicates no damage) in the KTF is expressed by the following differential equation: where n is the fatigue damage, n 0 is the initial constant, λ is the shape parameter ( [14,31] suggested value of λ = 9), and K b is the fatigue crack formation and growth rate, which is dependent on the durability equation [30]: where τ is the failure time and τ 0 = 10 -13 s is the characteristic period of atom oscillation in the solid body.Equation ( 12) can be easily adopted for the PD model by calculating the fatigue damage variables n for each PD bond.Then, the rate K b , defined as the PD bond breakage rate, using Planck's law, is expressed as [5]: where k and h are the Boltzmann and Planck constants, respectively; T is the temperature; U and γ are the process activation energy and activation volume, respectively ( U and γ can be found from the material S-N curve); and σ xx′ is the PD bond stress.PD (9) bond, connecting points x and x', stress is calculated as the average stress at point's x and x', namely, σ xx′ = σ x +σ x′ 2 .By assuming that the maximum PD bond cyclic stress σ max_xx′ occurs at time t and that the minimum stress σ min_xx′ occurs at time t 1 and relating σ min_xx′ = Rσ max_xx′ to t − t 1 = N f , where R is the stress ratio and f is the cyclic loading frequency, integrating Eq. ( 12) finally leads to KTF fatigue damage of the PD bond x'-x equation [14]: where n Ixx′ is the initial PD bond x'-x fatigue damage (if N = 1, n Ixx′ = 0 ).The initial constant n 0 is found from the definition of the process initial and end conditions [29], resulting in = 1 , and for λ = 9, n 0 = 1.771.Then, the PD bond "remaining life" parameter q(x, x'-x, N) can be related to the fatigue damage variable by the following equation: KTF Eq. ( 15) better describes the fatigue process, including temperature, stress ratio, stress value, and loading frequency variables, than does the empirical coefficient-based Paris [6] law approach.Nevertheless, increasing the stress ratio R in Eq. ( 15) reduced both the stress amplitude σ a = 1 2 σ max (1 − R) and the number of cycles to failure, while experimental testing [32] revealed that the number of cycles should increase.Fertig and Kenik [32,33] modified Eq. ( 15) for a composite including energy losses proportional to the square of the cyclic loading stress amplitude.Energy losses affect the process temperature, and then the total temperature term T in Eq. ( 15) is rewritten as: where T * is the ambient temperature, Δt is the cycle time and ψ is the material pro- portional constant.A similar approach described in subchapter 3.2.2based on cycle (15)  hysteresis loop areas is applied in this study to model asymmetric cyclic loading in metals.

PD-KTF model numerical implementation aspects
The PD equation of motion (Eq. 1) is solved numerically by applying spatial and time integrations.Because the number of material points is discrete in the numerical PD model, points near the surface of the PD horizon of point x are not fully embedded in the horizon and have truncated volumes (Fig. 3a).Additionally, material points at the boundaries do not have full horizons (Fig. 3b), which results in different material properties at the body boundaries (surface effect) from those in the bulk.To address these problems, volume corrections and surface corrections [35] based on introduced correction coefficients are used.
The PD surface effect takes one size of the PD horizon δ from the body surface for BBPD and twice for SBPD, and after the correction, it cannot be completely (only up to 90%) eliminated for complex loading.Notably, for the fatigue model, the remaining surface effect can be beneficial because fatigue cracks start on the surface of the body due to reduced body stiffness at the boundaries in the PD model and stretch; thus, the strains and stresses are higher, leading to earlier cyclic failure.
Time integration can be performed by explicit and implicit solutions [42].The implementation of explicit time integration is simpler, although the algorithm is only conditionally stable, and the time step value should be calculated according to the PD model discretization (PD grid spacing and PD horizon) and material (density and stiffness) parameters [35].For the BBPD 3D model, the stable time step is equal to the following: where V x = v cx V is corrected volume of the PD point x (Fig. 2) and f s ≤ 1 is the safety factor.To adapt the PD equation of motion to a quasistatic solution, fictitious damping is introduced into the explicit integration of Eq. (1).
Fatigue numerical modelling consists of a set of static simulations.First, the PD model is run for static simulation to achieve the maximum positive loading σ max_xx′ . .Then, (18)  using the relation σ min_xx′ = Rσ max_xx′ , the number of cycles N to failure of the maxi- mum loaded bond is calculated from Eq. ( 15).After N cycles, this maximum loaded PD bond is broken, while the fatigue damage of the other PD bonds is calculated according to Eq. ( 15).Because after N cycles there are broken PD bonds, the static simulation should be run again, and then the process is repeated.
The fatigue process consists of fatigue crack initiation and growth phases, which should be replicated in the model.The process activation energy U and activation vol- ume γ for the fatigue crack initiation phase are selected by comparing the number of cycles to failure calculated according to Eq. ( 15) and predicted according to the material S-N curve, as shown in Fig. 4. Usually, it is very difficult to find U and γ values for all stress levels; therefore, U and γ are selected for specific stress ranges (1, 2, …, M in Fig. 4) only.Although the material S-N curve corresponds to macroscopic fatigue damage, including both initiation and crack growth phases, macrocrack initiation is treated as material fatigue failure at stress concentrator PD point x and can then be predicted from the material S-N curve [21].Once the macrocrack is initiated, it starts to grow.The transition criterion between the crack initiation and growth phases is selected as the PD bond damage ratio φ x .According to the literature [14], if φ x ≥ 0.5 is achieved at point x, then the fatigue crack growth phase can start at point x.Because fatigue crack growth is faster than crack initiation, U and γ for the crack growth phase should be updated as U = k U •U and γ = k γ •γ.The coefficients k U and k γ are usually selected from the experimental relationship of the crack length and number of cycles by trial and error to ensure the same fatigue crack growth in the PD model (Fig. 4b prediction P) as in the experiment.Final failure is considered to be the maximum stress or reaction force in the model drop by a defined ratio.The testing protocol [43] suggests that fatigue experimental testing should end when the maximum stress decreases by more than 25%; thus, the same criteria are applied to PD fatigue process simulation.
The generalized KTF-PD fatigue modelling scheme (based on previous studies) [14,31]) can be explained by the following steps: 1. Run a static PD simulation at a given maximum positive loading; 2. Calculate the strain tensor of each point in the PD model; 3. Calculate the stress tensor of each point in the PD model; 4. Determine an equivalent stress of each point in the PD model.For uniaxial loading, the absolute longitudinal stress can be taken as the equivalent stress; 5. Calculate each PD bond stress as σ xx′ = σ x +σ x′ 2 ; 6. Calculate the broken PD bond ratio φ x .If φ x ≤ 0.5, then it is the fatigue crack initiation phase; otherwise, it is the fatigue crack growth phase.In the fatigue crack growth phase, update U and γ as k U •U and k γ •γ, respectively; 7. Calculate the number of cycles N according to Eq. ( 15) to the maximum stressed bond failure.Calculate the fatigue damage parameter n of each PD bond by taking this N value.Increase the number of cycles N and remove bonds for which n ≥ 1. 8. Go to step 1 and repeat steps 1-7 until the maximum stress reduction is greater than 25%.
testing machine with a maximum force of 10 kN.According to the standard ASTM E1012 [45], the bending strain in the specimen gauge area cannot exceed 5% of the axial strain (percent bending ratio).The INSTRON E10000 construction not allowed to achieve required percent bending, resulting more than the 10 times shorter specimen fatigue life than predicted from testing guidelines [43].To solve this issue, a special alignment flange, shown in Fig. 5a, was used.It consists of four alignment adjustment screws (the front side screws are marked "2" in Fig. 5a) and is mounted on the machine base with 3 screws (marked "1" in Fig. 5a) before the force cell (Fig. 5a).The specimen centring procedure was performed in the following order: 1) screws "1" and "2" and the bottom grip of the machine were released.(2) The specimen was inserted into the top machine grip, and the top grip was clamped.
(3) The bottom machine grip was clamped.
(4) The machine was set for a very small displacement of approximately 50 µm.This results in the effect that the entire bottom part of the machine stud hangs on the specimen and gravity aligns the specimen.( 5) All 4 alignment adjustment screws "2" were gently tightened until they touched the internal part ("3" in Fig. 5a) of the adjustment flange.( 6) The upper machine grip was slowly returned to the initial position, and all the mounting screws "1" were tightened.With the help of such an adjustment flange, early specimen fatigue failure did not occur; thus, the required specimen alignment was achieved.Strains were measured using an external extensometer (Fig. 5a) mounted on the specimen.The testing temperature was 20 °C (293 K), and the cyclic loading frequency was 0.56 Hz.Strain-controlled fatigue tests were run at symmetric and asymmetric cycles (shown in Fig. 6) and are as further defined as follows: (1) symmetric constant amplitude cycles with ε max = 0.18% strain amplitude; (2) POL-variable amplitude loading, which consists of blocks of one symmetric ε max = 0.6% cycle and 200 asymmetric cycles of strain amplitude ε = 0.18% when the mean strain ε = − 0.42%; and (3) PUL-vari- able amplitude loading, which consists of blocks of one symmetric ε max = 0.6% cycle that starts with compression and then 200 asymmetric cycles of strain amplitude ε = 0.18% when the mean strain ε = 0.42%.A 0.18% strain results in an approximately 210 MPa stress (slightly higher than the material yield limit); thus, the fatigue process can be identified as an intermediate for symmetric cyclic loading and a low cycle for asymmetric cyclic loading.

PD model
The PD fatigue model is implemented in the MATLAB software framework by using the principles of open-source PD MATLAB code, available in [46].The model of the standard specimen is shown in Fig. 7a.It consists of 726 material points (17 points per specimen diameter), and the PD horizon size is selected δ = 3.1�x .To significantly reduce the computational time, only the gauge area of the specimen is modelled.To prevent false failure at the boundaries, the PD bond failure of the material is active only in the middle of the model (Fig. 7a).The model is fixed at the bottom by the δ layer of nodes and loaded by prescribed displacement at the top, resulting in the required strain value (Fig. 7b; 0.0018 for symmetric cyclic loading) given in the testing conditions [43].First, static simulations were run to test the PD model (model validation according to the displacement and strain fields at static loading is shown in Fig. 7b).To obtain the quasistatic solution from the explicit time integration of the discretized PD Eq. ( 1), an artefactual damping of 1.6•10 4 kg/s was introduced.Selected 150 times steps result in the fastest possible static solution, ensuring correct displacement and strain fields at a given loading.The 316L steel tensile curve segment was taken from the experimental first cycle hysteresis loop, as shown in Fig. 7c.The material elastoplastic behaviour was evaluated by using the linear and power law tensile curve approximation implemented in the force-stretch and stress-strain relations (Fig. 7c, d), where the stress plastic region was described as σ = -0.242•ε−0.891 + 295 MPa if ε > σ y E = 6.6 • 10 −4 (Fig. 6c).The equiva- lent material plastic stretch values s p = 6.98 • 10 −4 and 316L stainless steel critical stretch value s c = 8.23 • 10 −3 calculated from Eqs. ( 5) and ( 7) were used (Fig. 7d).

KTF symmetric cyclic loading model
A fatigue model is created based on the static PD model of the standard specimen by introducing the fatigue crack initiation and crack growth phases.Because the specimen in the model is perfect, e.g., without any defects, fatigue failure would occur as predicted from the material S-N curve at a given stress in all the models without any initiated cracks.To prevent this issue, the elastic modulus Gaussian distribution (E = 193 ± 3 GPa, values selected from experimental observations) is used to start the fatigue crack (Fig. 8a) by creating a point in the specimen where the deformation is highest; thus, the lowest number of cycles is necessary for local material failure (PD bond failure) predicted from the S-N curve to initiate the fatigue crack there.The fatigue crack imitation phase material activation energy U and activation volume γ are calibrated based on the ε-N NUREG/CR-6909 curve [34], as shown in Fig. 8b.The σ-N curve, necessary for KTF model calibration, is derived from the NUREG/ CR-6909 ε-N curve by calculating stresses from strains according to the material inhouse experimental tensile curve, as presented in Fig. 7c.
For this case, the material S-N curve is approximated with three line segments, namely, σ ≤ 250 MPa, 250 MPa < σ ≤ 290 MPa, and σ > 290 MPa, by selecting material activation energy U and activation volume γ values that result in minimum differences between the fatigue life calculated according to Eq. ( 4) and that predicted from the material S-N curve (Fig. 8b) at the same stress level.More line segments can be used for a better material S-N curve approximation, but for the NUREG/CR-6909 curve, the three-line segment approximation results in a 5.02• × 10 -5 remaining area, expressed in relative units, between the NUREG/CR-6909 curve and its approximation; thus, the approximation accuracy is sufficient.
Once a fatigue crack is initiated on the virtual specimen surface, the fatigue crack growth phase starts.The switching from crack initiation to crack growth phase criteria is based on the PD point x PD damage ratio φ x .According to the [14], the value φ x ≥ 0.5 is selected to determine the initiated fatigue crack.For the fatigue crack growth phase, the activation energy U and activation volume γ are updated as k U •U and k γ •γ, respectively.This can be carried out by ensuring the same fatigue crack length during the same number of cycles (the same crack growth rate as shown in Fig. 4 b) in the PD model and experiment.Because crack growth experiments were not performed, such in-house experimental data were not available.Alternatively, the experimental maximum stress vs. linear cycle plot, presented in Fig. 9, is used to Fig. 9 Determination of the coefficients k U and k γ from maximum stress vs. linear cycle plot determine the coefficients k U and k γ in the PD model.The increase in the maximum stress reduction rate shows the start of the fatigue crack growth phase, whose duration is ΔN cycles, as shown in Fig. 9.The experimental and simulation failure criteria (at the end of the fatigue crack growth phase) are considered to be greater than a 25% reduction in the maximum stress [43].
Selecting k U = 1.2 and k γ = 1 results in the same number of cycles ΔN (Fig. 9) in the PD model as in the experiment when a 25% maximum stress reduction is achieved.
The effect of material hardening according to the number of cycles in the PD model is evaluated in a manner similar to that in finite element software [47] by using the maximum stress versus linear cycles experimental curve approximation dependent on the material plastic strain.For this purpose, the polynomial logarithmic function σ = f (ε pl , N ) can be expressed as is applied in the PD fatigue model.Here, N is the number of cycles, and a 0 − a 5 are coef- ficients that depend on the plastic strain of the 316L stainless steel material, namely,  The selection of function σ = f (ε pl , N ) is based on an in-house experimental max- imum stress curve approximation only, and material coefficients, e.g., the Lemaitre and Chaboche model [48], were not taken from material datasheets.
Following the modelling order mentioned in "PD-KTF model numerical implementation aspects", the general PD fatigue simulation scheme is presented in Fig. 10.

Extended KTF model for asymmetric cyclic loading
PD models with the same static material parameters were applied to simulate fatigue failure under POL and PUL cyclic loading.Because the POL and PUL loading schemes (Fig. 6) consist of variable amplitude blocks (1 cycle of 0.6% strain and 200 cycles of asymmetric loading with a strain amplitude of 0.18%), once static and cyclic calculations were performed in the PD model for the ε max = 0.6% cycle, the next modelling step was a static deformation of ε = − 0.6% for POL loading (for PUL ε = 0.6%).This was necessary to evaluate the hardening effect of the material.When analysing the experimental hysteresis loops, pure kinematic hardening was assumed to be accurate enough in the PD model.Then, the strains and stresses at ε = − 0.6% for POL and ε = 0.6% for the PUL loading of each PD point were saved, a new yield limit was calculated according to the kinematic hardening law, and the PD model was run from these saved strainstress values for static and cyclic analysis for 200 cycles of asymmetric loading with a strain amplitude of 0.18%.The static simulation and cyclic bond failure evaluation were repeated not only for each PD bond failure but also after each loading block, as shown in the POL cycling loading model scheme presented in Fig. 11.PUL modelling principles are the same and thus are not given in separate schemes.
Because KTF Eq. ( 15) with calibrated U and γ values for the symmetric cycle results in a significantly shorter fatigue life in cycles for R > − 1, especially when R ≈ 0.5…1, it should be corrected.Similar to the composite, the missing term in Eq. ( 15) can be treated as the energy loss during cyclic loading, which increases the process temperature [32,33,49].When term kT denotes the thermal energy in Eq. ( 15), the cycle energy loss can be expressed proportionally to each cycle hysteresis loop area A, and then the term kT in Eq. ( 15) is replaced by kT + ψA , resulting in the following KTF-modified low-cycle fatigue damage equation: (23) a 3 = 687.1εpl + 0.539; (24) a 4 = −54.71εpl − 0.0239; (25) a 5 = 1.559ε pl + 0.0002654.(26) where ψ is a proportional constant and A σ y is the hysteresis loop area at the material yield limit.The fatigue process activation energy U and activation volume γ in Eq. ( 26) are calibrated by applying the same principles as in Fig. 8.The material S-N curve from Fig. 8 is approximated with three line segments.Each segment is considered to be accurate over the stress range of 40 MPa (segments from 210 to 250, 250 to 290 and 290 to 330 MPa, respectively), and then the hysteresis loop areas for the Eq. ( 26) term kT + ψA are calculated for each segment middle point, e.g., the line segment from 250 to 290 MPa → A 270 MPa , as shown in Fig. 12a.Measurements of the experimental hysteresis loop areas at the first cycle showed that the POL asymmetric loading loop area (orange colour in Fig. 12b) is 22% smaller than the symmetric loading loop area (A 230 MPa ) with the same maximum stress of 230 MPa.This also applies for PUL asymmetric loading, but in this case, the maximum stress is 270 MPa, and the PUL asymmetric loading hysteresis loop area (blue in Fig. 12b) is only 45% of the symmetric loading loop area when the maximum stress is 270 MPa.As shown in Fig. 12b, 0.78A 230 MPa and 0.45A 270 MPa hysteresis loop area values and stress ratios of R = − 1.5 and − 0.7 were used to evaluate the fatigue damage and number of cycles to failure according to Eq. ( 26) for asymmetric POL and PUL loading, respectively.A

Symmetric loading
The KTF-PD model for symmetric cyclic loading with a strain amplitude of 0.18% was first calibrated according to the NUREG/CR-6909 [34] curve without evaluating cyclic energy losses (Fig. 8).PD damage (ratios of the broken PD bonds to total PD bonds of each point) plots showing simulated fatigue crack paths on this model are presented in Fig. 13.The comparison of the simulation results to those of in-house and other laboratories [43] experimental tests is shown in Fig. 14.
According to the results in Figs. 13 and 14, the PD model, which was calibrated from the NUREG/CR-6909 curve, shows good agreement with the FrFr laboratory experimental results in terms of the number of cycles up to 25% of the maximum stress reduction (N25 in Figs. 13 and 14).The NUREG/CR-6909 curve predicts 171,517 cycles, while the PD model predicts N25 = 183,410.During the experiment, which was performed in an FrFr laboratory, a 25% reduction in the maximum stress was achieved after a total of 205,500 cycles.However, in-house experimental testing showed N25 = 342,870 cycles, Fig. 14 Comparison of the experimental and simulation results for symmetric ε = ± 0.18% cyclic loading which is far from the NUREG/CR-6909 curve.Nevertheless, experimental testing performed by the Japan Nuclear Energy Safety Organization (JNES) [34] predicted an average fatigue life of 321,085 cycles at ε = ± 0.18% cyclic loading, which correlates well with the in-house experimental result of 342,870 cycles.Thus, more experimental testing should be performed to possibly extend the average fatigue life of 316L steel nuclear power plant structures and establish a safe and prolonged exploitation period.
A comparison of the maximum experimental and simulation stress versus linear cycle curves is presented in Fig. 15.Here, the material kinematic hardening effect, simulated according to Eqs. ( 19)- (25), can be seen.Figure 16 presents the hysteresis loops of the 10, 1000, 100,000, 181,970 (simulation), and 341,970 (experiment) comparisons.Both the experimental and simulated hysteresis loops are in agreement according to their shape and area.

Asymmetric loading
Comparisons of the in-house and other laboratories [43] experimental and KTF-PD simulation results of both the conventional (without hysteresis loop energy loss evaluation) and improved (with loop energy loss evaluation) KTF-PD models for all loading types (symmetric, POL and PUL) are presented in Fig. 17.Each KTF-PD simulation took ~ 8 h of computational time when the simulations were performed on a computer with an Intel(R) Core(TM) 3.2 GHz processor and 16 Gb of RAM memory.
The newly proposed KTF method equation with material energy losses proportional to its loading hysteresis loop area validation is shown in Fig. 17a for the symmetric cycle.Compared with the conventional KTF approach, the modified KTF method results in a 3.8% longer fatigue life (190,170 versus 183,410), but both models are in very good agreement with those predicted from the N25 fatigue life (171,517 cycles) of the S-N curve and the FrFr experiment [43] (205,500 cycles).However, for POL loading (Fig. 17b), the conventional KTF-PD model better estimates the N25 fatigue life (147,060 cycles and the modified KTF method 185,280 cycles versus the predicted 158,719 cycles), but on the other hand, the modified KTF method result is closer to the VTT experiment [43] result (197,460 cycles) and still correlates well with the predicted N25 fatigue life [43] of 158,719 cycles.Moreover, the simulated N25 fatigue lives for POL and symmetric cyclic loading with both conventional and modified methods show that the POL N25 fatigue life tends to be shorter than that for symmetric loading, as also observed in the experimental results.The effect of the modified KTF method is mostly visible for PUL loading (Fig. 17c), where the simulated N25 fatigue life of 74,675 cycles correlates well with the predicted N25 fatigue life [43] of 70,029 cycles and the VTT experiment result of 63,980 cycles [43], while the conventional KTF equation estimates an extremely shortened N25 fatigue life equal to only 3185 cycles.In-house experimental testing for PUL loading resulted in significantly earlier failure after 23,614 cycles, similar to the results of the IRSN experiment [43].According to the results presented in Fig. 17, it is simple to conclude that the conventional KTF calibrated for symmetric cycles lacks accuracy for asymmetric cyclic loading when R > − 1.After the material energy losses proportional to its loading hysteresis loop area evaluation, the KTF method results in acceptable agreement between the simulated and experimental fatigue lives for all types of loading: constant amplitude symmetric and variable amplitude asymmetric with R < − 1 (POL) and even for R > − 1 (PUL).The average differences between the KTF-simulated fatigue lives for all types of loading are 10.8% and 40.2% for unmodified KTF.Nevertheless, more experimental testing is necessary to obtain the fatigue life variation and then better select the numerical fatigue model parameters, especially ψ.The simulated and experimental hysteresis loops for POL loading are plotted in Fig. 18 and for PUL loading-in Fig. 19.Each loop presents ε = 0.6% symmetric cycle, and the next asymmetric cycle has a strain amplitude of 0.18%.
Figures 18 and 19 show agreement between the simulated and experimental loops in terms of the maximum strain-stress values for symmetric cycles, but the selected stress-strain approximation leads to differences in terms of the shape of the curves.

Summary and conclusions
The kinetic theory of fracture under the peridynamic material model is able to predict fatigue failure with very good agreement with the experimental results.Moreover, with the help of integral (without spatial derivatives) peridynamics theory formulation, fatigue crack initiation and growth phases can be simulated without predefining the crack path.The cyclic material model parameters can be found only from a single Fig. 19 Comparison of experimental and KTF-PD simulated hysteresis loops for PUL loading experimental S-N curve of the material for symmetric loading (R = − 1).However, the kinetic theory of fracture fatigue damage equation (Eq.15) calibrated on the symmetric cycle S-N curve results in an incorrect fatigue life for low-amplitude asymmetric cyclic loading (R > − 1, especially R ≈ 0.5…1), and variable amplitude and asymmetric cyclic loading can still be simulated based on the same S-N curve of the material by extending the KTF-PD method with material energy losses during cyclic loading evaluation.It is assumed that energy losses are proportional to the area of the loading hysteresis loop.The modified KTF-PD method is accurate for different types of symmetric and asymmetric cyclic loading.The material kinematic hardening effect can also be captured in the KTF-PD simulation.

Fig. 2
Fig. 2 PD bond force-stretch relation (without failure of the material): (a) for purely linear elastic material (original Silling PD model) and (b) linear elastic PD bond force stretch relation

Fig. 3
Fig. 3 Numerical spatial integration issues: (a) truncated PD point volumes, requiring volume correction V i = v ci V and (b) not full PD horizon at the boundaries, requiring PD bonds stiffness corrections

Fig. 4
Fig. 4 Definition of activation energy U and activation volume γ : (a) from the S-N curve for the fatigue crack initiation phase and (b) updating U = k U •U and γ = k γ •γ for the fatigue crack growth phase

Fig. 8
Fig. 8 Modelling fatigue crack initiation: (a) probabilistic material elastic modulus distribution to define the crack start position and (b) material activation energy U and activation volume γ calibration from the S-N curve a 0 = 5457ε pl + 227.1;

a 1 =Fig. 10
Fig. 10 Schematics of the general PD fatigue model for symmetric cyclic loading

Fig. 11
Fig. 11 Schematics of the general PD fatigue model for POL cyclic loading

Fig. 12
Fig. 12 Modified KTF approach: (a) U and γ calibration from the S-N curve when evaluating cycle energy losses and (b) evaluating the experimental hysteresis loop areas for asymmetric cycles

Fig. 13
Fig. 13 Fatigue crack initiation and growth in the KTF-PD model

Fig. 15 Fig. 16
Fig.15 Comparison of the maximum stress vs. linear cycles and the simulated kinematic hardening effect

Fig. 17
Fig. 17 Comparison of experimental and KTF-PD simulation (conventional and extended with hysteresis loop energy losses) results: a symmetric cycle, b POL loading, and c PUL loading

Fig. 18
Fig. 18 Comparison of experimental and KTF-PD simulated hysteresis loops for POL loading Page 7 of 27 Vaitkunas et al.Adv.Model.and Simul. in Eng.Sci.