Fatigue degradation strategies to simulate crack propagation using peridynamic based computational methods

*Corresponding author http://dx.doi.org/10.1590/1679-78255022 Abstract The aim of this paper is to develop new computational tools to study fatigue crack propagation in structural materials. In particular we compare the performance of different degradation strategies to study fatigue crack propagation phenomena adopting peridynamic based computational methods. Three fatigue degradation laws are proposed. Two of them are original. Initially a cylinder model is used to compare the computational performance of the three fatigue laws and to study their robustness with respect to variations of discretization parameters. Then the fatigue degradation strategies are implemented in a peridynamic framework for fatigue crack propagation analyses. Both cylinder model and peridynamic simulations show that the third proposed degradation law is unique in its combination of high accuracy, high stability and


INTRODUCTION
Fatigue is a major cause of structural failure in many structures. Prediction of fracture propagation under cyclic loading is a challenging problem for two main reasons: from a theoretical point of view the discontinuous nature of fracture phenomena conflicts with the underlying mathematical structure of classical continuum mechanics, and from a more practical side the slow damaging effect of a large number of load cycles (high cycle fatigue) is not easy to model (Weinert and Gergely 2011). Traditional finite element method (FEM), or various modified versions of the FEM, have been used to simulate fatigue crack growth (Allix and Corigliano 1996;Branco et al. 2015;Point and Sacco 1996;Schellekens and Borst 1993). FEM models are based on the classical continuum mechanics theory, which cannot directly be used to simulate problems with discontinuities, such as cracks, requiring extra relationships that control the crack growth speed, direction and remeshing techniques. For example, the FEM with a remeshing technique is widely used to model fatigue crack growth (Branco et al. 2015). However, implementation of automatic remeshing techniques is difficult in 3D cases (Baydoun and Fries 2012;Boroomand et al. 2016;Boroomand et al. 2017), and is computationally costly to remesh geometries when multiple cracks exist. Interface elements equipped with cohesive zone models have been used Alfano and Crisfield 2003;Blackman et al. 2003;Chen et al. 1999;Crisfield and Alfano 2002;Ostergaard et al. 2011;Qiu et al. 2001) in the cases in which the crack path is predictable. The extended finite element method (XFEM) (Belytschko and Black 1999;Bidokhti and Shahani 2015;Dolbow and Belytschko 1999;Shi et al. 2010) allows cracks to pass through the elements leading to better approximations of crack paths without remeshing.
Nevertheless, one common problem for XFEM is that additional criteria are needed to study crack coalescence or branching.
Peridynamics, a newly presented continuum theory, has been developed to deal mainly with dynamic fracture problems by Silling in (Silling 2000). It is a non-local continuum theory proposed with the specific goal to tackle the shortcomings of the classical theory. Peridynamics consists of integral equations and does not need any spatial differentiation, for these reasons it is proper to describe problems affected by discontinuities (Bazazzadeh et al. 2018;Bobaru et al. 2016;Madenci et al. 2016;Shojaei et al. 2019). One of the significant features of the peridynamic approach is that it allows for spontaneous formation, interaction, and growth of discontinuities in a consistent framework (Xu et al. 2008). Using this theory, promising results have been attained in elastic bar (Aguiar et al. 2018) dynamic brittle fracture in glass Bobaru 2010, 2011;Bobaru and Ha 2011;Hu et al. 2013), fiber-reinforced composites (Ghajari et al. 2014;Hu et al. 2012), functionally graded materials (Cheng et al. 2015) and corrosion damage (Chen and Bobaru 2015;Chen et al. 2016).
A peridynamic model for fatigue cracking has been proposed in the recent past by Silling and Askari (Silling and Askari 2014), which enables the simulation, in a single model, of the three phases of fatigue failure: crack initiation, fatigue growth and final failure controlled by quasi-static crack growth. In (Silling and Askari 2014), the dynamic relaxation method is used to obtain the static solution, and fatigue crack growth. References Zhang et al. 2016) study fatigue crack propagation obtaining the static solution adopting a conjugate gradient energy minimization method, which is matrix-free and much faster than dynamic relaxation method, see (Le 2014). In a recent work (Hu et al. 2015;Hu and Madenci 2017) peridynamics is applied to simulate the fatigue life and residual strength in composite laminates. Moreover, by using peridynamic approach, mixed-mode fatigue crack propagation has been investigated in (Jung and Seok 2017). Other studies to simulate fatigue phenomena based on the peridynamic theory can be found in (Baber and Guven 2017;Jung and Seok 2016;Oterkus et al. 2010).
In the present paper a cylinder model (Galvanetto et al. 2009) is used to efficiently compare the computational performance of three fatigue laws. Three fatigue degradation strategies have been implemented in conjunction with a bilinear constitutive law. The computational performance of the three fatigue degradation strategies is examined by making use of various discretization parameters. Then the same constitutive laws and fatigue degradation strategies are implemented in a peridynamic based computational code. Due to the non-local nature of peridynamic methods, they are computationally more expensive than those based on classical continuum theories. In the conventional peridynamic approach, a constant horizon size and a uniform grid are assumed all over the domain. Accordingly, the efficiency of the peridynamic approach is strongly affected when a fine grid spacing is used. There have been many attempts, so far, to maintain the computational accuracy and reduce the cost by coupling peridynamics with local theories Seleson et al. 2013;Shojaei et al. 2016;Galvanetto et al. 2016;Shojaei et al. 2017;Zaccariotto et al. 2017;Zaccariotto et al. 2018), or by using parallel computing techniques (Lee et al. 2017;Mossaiby et al. 2017) and adaptive refinement algorithms (Dipasquale et al. 2014;Gu et al. 2017;Ren et al. 2016;Shojaei et al. 2018). Among the aforementioned techniques the FEM-PD coupling enhanced with adaptivity represents an efficient approach the use of which will be explored in future works.
The paper is organized as follows: in section 2 kinematics and statics of the cylinder model are defined. Then, the bilinear interface constitutive law is recalled in section 3. In section 4, three fatigue degradation strategies are described. In section 5 the computational performances of the three fatigue laws are compared. Finally the implementation of the three fatigue laws in a peridynamic code is illustrated in section 6 where the results they provide are also compared. Section 7 gives the conclusions.

THE CYLINDER MODEL
The mathematical model consists of a cylinder of radius R which rolls on a horizontal surface initially composed of two coinciding zero-thickness layers, see Figure 1(a) and ref. (Galvanetto et al. 2009). One of the layers is fixed and the other adheres to the surface of the cylinder at the contact point C. The upper layer remains adhered while the cylinder rotates and the contact point moves to the right (see Figure 1(b)). In the initial state the contact point is at 0 x = (see Figure 1), the rotation, θ, of the cylinder is zero and the line OA, on the cylinder cross section, is vertical. The value of the rotation of the cylinder is defined as the angle between the line OA and the vertical line connecting O to the current contact point. If a clockwise moment is applied, the cylinder will rotate in the same direction and it will move to the right of its initial position and one layer will be lifted and separated from the other, as shown in Figure-1 (b). An adhesion stress acting between the two layers opposes the motion of the cylinder. The model represents one of the two symmetric halves of the tip of a crack opening in mode I as depicted in Figure

Continuum equations of the model
The model has only one degree of freedom, the rotation of the cylinder θ . The coordinate of the contact point con X represents the crack length, which is linked to the rotation θ as The applied moment, a M , is in equilibrium with the reaction moment, t M , due to the stresses in the interface 0 ( )( )

Discretization of the model
The motion of the cylinder detaches the two layers. The adhesion between them can be represented as a discrete system of vertical springs, with uniform spacing Δl, that oppose the separation, Figure 1.b. The springs are initially unstressed but become stressed as the two layers separate behind the advancing contact point.
A constitutive law for the springs has to be defined to describe the mechanical behavior of the interface which is progressively damaged: the force in the spring eventually reduces to zero (representing complete failure) as the layer separation increases.
The value of R is to be chosen as very large compared to the displacement at which the stretched springs fail so that the springs can be considered to remain vertical while they are stressed. is known and given the spring spacing Δl, it is easy to identify the number of springs that become active (n) and their coordinates ( 1, 2, 3, Figure 3 shows that simple geometry is sufficient to define the elongation δi of the i-th spring and its distance i d from the contact point. Since the springs have a zero initial length, the elongation will be defined as the vertical displacement ( i δ ) of the top of the spring. Moreover it was previously mentioned that the spring eventually degrades as the displacement is increased and the stress finally reduces to zero. If c δ defines the critical displacement at which failure occurs, then the angle α at which such a displacement is achieved is given by (Figure 3):

Equilibrium
As the cylinder advances, a number of springs become active. The interface stress ( i σ ), which is a function of the opening displacement ( where Gc is the critical energy release rate. The moment Mc of the cylinder model can be found by imposing static equilibrium for a rotation of magnitude α.

INTERFACE CONSTITUTIVE LAWS
As it was mentioned previously, a constitutive law is required to describe the behavior of σ as a function of δ . The constitutive law should be defined in such a way that for every value of δ there is a unique value of σ . In this research the bilinear interface constitutive law has been considered Alfano and Crisfield 2003;Crisfield and Alfano 2002;Galvanetto et al. 2009;Peerlings et al. 2000;Qiu et al. 2001;Robinson et al. 2005). There are two parameters, additional to the critical displacement value, c δ , that are needed to define an interface constitutive law: the critical energy release rate ( c G ) and the maximum stress value ( 0 σ ). In the present case of the bilinear law the elastic limit relative displacement ( 0 δ ) is also required (Galvanetto et al. 2009;Peerlings et al. 2000;Robinson et al. 2005). The critical energy release rate is a physical quantity used in Linear Elastic Fracture Mechanics (LEFM) for defining a criterion of propagation of the delamination in the interface. Its value is measurable in lab tests. It is indirectly used in the definition of the constitutive laws by equating the areas under the stress-displacement curves to the critical energy release rates. Therefore, the total energy dissipated during the delamination is correctly computed at each point although it is not released instantaneously as is assumed in LEFM .
δ0 is the maximum displacement for which the spring behavior is still linear elastic, whereas the ratio ( 0 0 / σ δ ) is clearly linked to the stiffness of the interface.
The constitutive behavior for monotonic increasing relative displacement δ will be described below.
The bilinear interface law, as seen in Figure 4, is divided into three parts: linear elastic, stiffness degradation and failure. In the first part a linear-elastic behavior is assumed for a relative displacement that is increased from zero, where there is no stress ( 0 σ = ), until it reaches the elastic limit 0 δ . At the elastic limit the stress assumes its maximum value, 0 σ . The linear-elastic part of the law is described by the following expression: where K is the initial stiffness of the interface and is given by The second part of the interface constitutive law defines the behavior between the elastic limit and the critical displacement c δ . Once the elastic limit is exceeded the interface starts to be degraded. This degradation is described by the following relationship: Fatigue degradation strategies to simulate crack propagation using peridynamic based computational methods The third and last part of the interface constitutive law is for relative displacement values that are equal or larger than the critical displacement value. When this critical value is reached or exceeded the element fails irreversibly so that no stress can be carried by the interface: As it was stated previously the critical energy release rate c G is equal to the area under the δ σ − curve, so in this Therefore of the four parameters c G , c δ , 0 δ , 0 σ only three are independent. The constitutive law of the interface can be written as well in terms of a damage variable D as: where the damage variable is: The increment of the opening displacement from 1 δ to 2 δ ( 2 1 δ δ > ) causes an increment of the damage variable:

FATIGUE DAMAGE
If the fatigue phenomenon is involved a time-like independent variable has to be considered, the number of cycles N of the external load. The applied load history consists of two components: quasi-static ramping phase and fatigue phase as shown in Figure 5. In order to simplify the calculation process some assumptions are made (Peerlings et al. 2000). The first one is related to the cyclic load which is assumed to be sinusoidal and oscillating between zero and a maximum value. Since it would be very computationally expensive to simulate the relative displacement history for each spring, the load numerically applied to the structure, or to the cylinder model in this paper, is taken as constant and equal to the maximum value of the actual cyclic load as shown in Figure 5. Consequently, the relative displacement calculated at the crack tip will be the envelope of its true cyclic variation with time, as shown in Figure 5(b) (Galvanetto et al. 2009;Peerlings et al. 2000;Robinson et al. 2005).
We assume that the total damage rate in a spring can be obtained as the sum of the static and the fatigue damage rates as follows s D  is the rate of static damage due only to the increase in magnitude of δ ; f D  represent the fatigue damage increment due to the repetition of the load cycles and is greater than zero even if δ is constant. As the novelty of this work is in the fatigue part we use the same static damage rate based on the bilinear constitutive law for all different fatigue approaches.
Three strategies are considered to estimate the fatigue part of the damage.

Static damage rate
The static damage rate with respect to time can be obtained from Equation (14) as One may integrate Equation (17) over an increment of the number of load cycles ( N ∆ ) as . Equation (18) provides the same expression as Equation (15) and therefore clarifies that the static damage component is due only to the increment of the opening displacement.

First fatigue degradation strategy
In the first strategy, originally presented in (Galvanetto et al. 2009;Robinson et al. 2005) the fatigue rate is assumed to have the following form: C , λ and β are parameters to be determined by fitting of experimental data. The physical meaning of Equation (19) is that the fatigue damage rate is bigger for bigger relative displacement δ and for higher values of the current total damage.
Integrating Equation (19) over a number of loading cycles ( N ∆ ), one may write the increment of fatigue damage as Fatigue degradation strategies to simulate crack propagation using peridynamic based computational methods Moreover, based on the reference (Galvanetto et al. 2009), D µ and µ δ are taken as Where µ value is 0.7. Additionally, by putting (18) and (20) into Equation (16), the total damage increment, based on the bilinear constitutive law and stiffness degradation damage (see Figure 6), can be obtained as Note that in usual numerical procedures the only unknown variable in Equation (22) is N N D +∆ , it appears on both sides of Equation (22) which makes it implicit and nonlinear.

Second fatigue degradation strategy
The second fatigue degradation strategy, originally proposed by the authors in (Zaccariotto et al. 2015b), is based on the reduction of the value of the failure displacement, due to the fatigue phenomenon ( Figure 7). It can be assumed that the failure elongation of a spring in each increment of load cycle, N ∆ , is a function of displacement as follows.
Where B and n are parameters related to the material behavior when it is subjected to fatigue load cycles and c N N δ +∆ is the failure displacement of a spring after N N + ∆ loading cycles.
The total damage at the end of each load cycle increment can be obtained as Accordingly, the total damage increment in this approach is given by Latin American Journal of Solids and Structures, 2019, 16(2), e163 9/31 In the above formula if 0 N ∆ → , only the static damage takes place due to the increment of external load so that Equation (25) is simplified to which is consistent with Equation (18). In Equation (24), the dependence on the number of cycles is implicit in c N N δ +∆ (see Equation (23)) and in the value of the current displacement. Since N N D +∆ doesn't appear on both sides of Equation (25), this numerical approach, differently from the first strategy, is explicit.

Third fatigue degradation strategy
In the third law, which has not presented before, it is assumed that the derivative of fatigue damage with regard to N is given in Equation (27) In the above formula, A′ and m′ are constants related to the material behavior while δ represents the elongation of the spring. Equation (27) differs from equation (19) because in (27) the damage rate does not depend on the current value of total damage D. Updated fatigue damage for crack growth can be determined as (see references (Bak et al. 2016;Peerlings et al. 2000)) Using a mid-point rule, one may approximately calculate the integral of the fatigue damage rate as Latin American Journal of Solids and Structures, 2019, 16(2), e163 10/31 Accordingly, by combining Equations (29) and (18) in Equation (16), the increment of total damage due to an increment N ∆ of the number of cycles is given by the following expression, explicit in the unknown N N D +∆ : Note that in the third law, as in the first one, damage directly affects the stiffness of the springs as shown in Figure  6.
Once the sum of static and fatigue damage increments has been computed, Equation (13) is used to compute the stress in the spring.

COMPARISON OF THE PERFORMANCE OF THE THREE FATIGUE DEGRADATION STRATEGIES APPLIED TO THE CYLINDER MODEL
In this section we compare the computational performances of the three fatigue degradation strategies proposed in section 4. The fatigue phenomena considered are described by the following Paris law, based on the experimental data gathered from (Asp et al. 2001), (the specimen is under pure mode I loading): where I C and I m are material parameters, con X and N represent the crack length and number of load cycle respectively. G , c G and 0 G express the maximum, critical and threshold energy release rates respectively. For the cylinder model, using Equation (7), one may rewrite Equation (31) Table 1, provided in reference (Bak et al. 2017). In a log-log graph representing / con dX dN vs / c G G Paris law is a straight line, as shown by the red line in Figure 8.   . Given the initial static ramp of the load phase, the crack fatigue propagation starts from a value different from zero.

Fatigue laws
Model specific parameters As depicted in Figure 9, the fatigue crack propagation is divided into two phases. At the beginning of crack propagation, one may observe a slow increase of the crack length, and then the crack grows more quickly with a constant slope. It is also apparent that the slope . All three proposed approaches are able to produce the same slope, . The three graphs show the stress-elongation curve (σ−δ curve) for the same spring in the constant slope propagation phase. The three approaches generate a similar curve due to the simultaneous action of increasing displacement and fatigue cycles. For all cases the area under the σ−δ fatigue curve approximates well the green area which represents half of the area under the bi-linear law. In many practical problems, the applied G is not constant. So, following what was done also in (Bak et al. 2016), we compare the performance of the three fatigue degradation strategies when used with a variable G (i.e a a M in the case of the cylinder model). We alter the applied moment instantly as shown in Figure 11. The response to this energy release change can be observed in the crack length vs. the number of load cycles as depicted Figure 11. It can be noticed that there is a good agreement between all proposed approaches and the Paris law. The results of the Paris law are obtained by integrating numerically Equation (31) to evaluate the crack length as a function of energy release rate as (see Reference (Bak et al. 2016 Moreover, in many engineering problems gradual energy release changes occur. Assuming the energy release rate as a function of crack length, one may change the loading condition in a way to obtain an increasing/decreasing energy release rate as depicted in Figure 12.

Stability of the results with respect to variations of discretization parameters
The results presented in the previous section seem to suggest that the performance of the three degradation strategies are rather similar; the second law has a worse performance in the case of increasing Ma; and the first law requires a longer computational time than the other two, but the overall behavior of the three laws does not exhibit apparent differences. In this section we want to compare the stability of the results provided by the three laws with respect to variations of two important discretization parameters: l ∆ , the spring spacing, and N ∆ , the increment of load cycles in one step. The best approach is the one for which the slope / con dX dN is affected the least by changes in l ∆ and N ∆ . All parameters are given in Table 2. The applied moment is set to be  Figure 14 represent different percent errors in the computed slope with respect to 'exact slope'. Analyzing these error-plots, we observe that as l ∆ and N ∆ increase the error in general rises. The region of small error is much wider for the third law and much smaller for the second. Therefore, the third law, not only is computationally cheaper than the other methods for the same values of Δl and ΔN, but also can be used with larger values of l ∆ and N ∆ . Therefore, the cylinder model suggests that the third fatigue law is preferable. In the next section, we verify such a suggestion by implementing the three laws in a peridynamic-based code. In particular, the bilinear constitutive law equipped with the three fatigue degradation strategies will be adopted as constitutive law that defines the force in the bond as a function of the bond stretch.

PERIDYNAMICS FATIGUE MODEL
Bond based peridynamics (BBPD) is a non-local theory of continuum proposed in Reference (Silling 2000). Since the formulation of the theory avoids differentiations and is only based on spatial integration BBPD is particularly suited for describing crack propagation problems. For a point, x , and time, t , the peridynamic static equilibrium equation is written as (see Figure 15) In which, Ω is the body domain, f is the pairwise force function in each peridynamic bond between node x and x , u , is the displacement vector field, x H is a local spherical region of radius H δ , in which the integral is defined. ξ and η are the initial relative position and current relative displacement respectively: The stretch of the bond can be expressed as

Spatial Discretization of the model
We discretize the solution domain with a grid of points called nodes (see Figure 16), to approximate the peridynamic equation.
By adapting the one-point Gauss quadrature rule to proceed with the spatial integration, one can discretize Equation (34) as Equation (38) expresses that the volume of nodes close to the boundary of the neighborhood falls only partially within the horizon of the source node. Additionally, the distance between two adjacet nodes is referred to as the grid spacing.
As proposed in (Silling 2000), in the prototype microelastic brittle material (PMB), the pairwise force function f and the stretch s are related according to the following equation in which c is the micro-modulus that can be interpreted as the bond stiffness, ( ) µ ξ is a history-dependent function which takes the value of 1 for active bonds and 0 for broken bonds. Furthermore, for "small deformations" the magnitude of the bond force, f , is approximated as the second part of the Equation (39) in which, h t represents the thickness. The value for Poisson's ratio for the 3D and plane strain cases is fixed to 0.25 ν = ; however, for the plane stress cases this value is restricted to The limitation of Poisson's ratio is resolved in (Silling et al. 2007) and (Warren et al. 2009) by using the state-based formulation of the peridynamic theory.

Pairwise force function
Based on Reference (Silling and Askari 2005), for a PMB material the stretch in a bond is related to the pairwise function f by Equation (39). When the stretch reaches its limit value, 0 s , the bond breaks so that the global behavior of the system is nonlinear. In this section the constitutive law of the bond ( ) f f s = is assumed to be bilinear of the type presented in Section 3. The stretch of the bond plays the same role as the elongation of the spring. 0 s is the elastic limit, and c s is the failure limit of the stretch, as shown in Figure 17. According to reference (Ha and Bobaru 2010), following expressions are assumed: Plane strain 12 where, 0 G , is the energy required to start the damaging process. The ratio between c G and 0 G , called r k , is defined as (see Reference (Zaccariotto et al. 2015a)): Moreover, similar to the cylinder model, the increment of total damage in each bond for the first fatigue strategy is given by (see Equation (22) Similarly, for the second fatigue law, the total damage of each bond at the end of a load cycle can be obtained as (see Equation (24) Finally, for the third fatigue law, the increment of total damage in each bond, can be written as (see Equation (30) The damage level in a point x , at cycle number N , can be expressed by In Equation (51), φ , indicates the ratio of the number of broken bonds to the total number of bonds originally connected to the point x . The damage level, φ , takes a value between 0 and 1, 0 1 φ ≤ ≤ . 1 φ = means there is no interaction between the point and all surrounding points within its horizon, while 0 φ = denotes an undamaged state.

Nonlinearity of the system
The problem becomes nonlinear after the bonds get damaged, therefore, the stiffness matrix of the structure is not constant. Various techniques have been used to satisfy the equilibrium in the whole structure. In this paper, the Newton-Raphson (N-R) method is adapted to solve Equation (37) and update the stiffness matrix in each iteration.
Using the conventional finite element concept, one is able to write the relationship between external and internal forces as follows: in which, K represents the global stiffness matrix of the structure, u and ext F stand for the nodal displacements and external force vector respectively.
An out of balance force in the whole system is generated due to the nonlinearity of the problem. Therefore, out of balance force vector is represented by The equilibrium is satisfied when the norm of g is lower than a given tolerance. According to the Newton-Raphson solution, the derivatives of g with respect to the nodal displacement vector, u , have to be calculated as follows

Loading and boundary conditions
The application of boundary conditions in peridynamics is a significant issue: since integral operators, rather than partial differential operators, are used to calculate the peridynamic equilibrium equations, the traction boundary conditions should be applied over a volume (Madenci and Oterkus 2014).
In order to converge to the equilibrium solution, static forces are applied gradually. After static equilibrium is satisfied, the fatigue phase initiates in the model, as shown in figure 5(a).
Numerical models with static damage have been comprehensively investigated in (Zaccariotto et al. 2015a), where deformations of the structures have been compared to the exact and FEM numerical solutions. Since the main contribution of this work is on the fatigue damage, the authors just present the examples with static and fatigue damage to further validate the fatigue degradation strategies in a peridynamic framework.

Double cantilever beam peridynamic example
In order to verify the indications of the cylinder model, a double cantilever beam is simulated with a BBPD code as depicted in Figure 18.  ; and different increments of load cycles, ΔN=1000, 10000 and 50000 for the first and the third proposed fatigue laws, in contrast, due to convergence problems we used smaller values of 500, 1000 and 10000 for the second fatigue law. (See Table 3).  In this example, the crack length with respect to the number of load cycles, N , is computed by using the three fatigue strategies (See Figure 21).   Figure 21 shows the second law to be much more sensitive to parameter variations than the other two: it does not work for ΔN= 50000 and for smaller values of ΔN provides results that are clearly more scattered than those produced Latin American Journal of Solids and Structures, 2019, 16(2), e163 25/31 by third law. Moreover the comparison between Figure 21(a) and 21(b) shows that if the lines in Figure 21(a) obtained for ΔN= 50000 are neglected, since there is no equivalent counterpart in Figure 21(b), the results of the second law are worse even than those of the first law. Finally comparing Figure 21(a) and 21(c) confirms the findings of the cylinder model: the third law is much more stable than the other two and provides very similar values of the slope of the line 'crack-length vs number of cycles' in a reasonably wide range of values of H δ and ΔN. Additionally, Figure 21 shows that the third fatigue law is less sensitive than the other two to the horizon size, H δ . Both the cylinder model and the DCB discretized with the peridynamic based code confirm that the third law is the best of the three considered in the present work because it is more stable with respect to the variations of the discretization parameters and it is cheaper from a computational point of view.

Experimental validation of the model
To further validate the presented numerical methods, we simulate a DCB, which is taken from an experimental study (Juntti et al. 1999), using BBPD code. Due to the good convergence and high computational speed of the third fatigue law which represented in previous sections, hereinafter we only use the third law for the rest of the examples. In this example, we investigate the effectiveness of the selected fatigue degradation strategy by comparing with the experimental data (Juntti et al. 1999). The specimen used for this simulation is shown in Figure 22 Where, I , is the second moment of area of the specimen arm.
In a log-log graph representing

Center cracked beam with a hole
In this section, fatigue crack propagation is simulated in a single notched beam with a hole by using the third fatigue law (see Figure 24). Experimental data on this test can be found in (Miranda et al. 2003). As shown in Figure 24, the beam is loaded in two points at the top and supported in two other points at the bottom. This problem has been solved under plain strain conditions with the physical parameters of:  The solution domain is discretized with a 106×438 uniform grid of nodes which results in 45,390 nodes. This example has been solved up to 5 6.7 10 × load cycles using 50000 N ∆ = cycles. As shown in Figure 25, the crack path grows toward the center hole. Moreover, the crack path is in a very good agreement with experimental data in (Miranda et al. 2003

CONCLUSIONS
The paper compares three fatigue degradation strategies to be used in the simulation of fatigue crack propagation. The adopted constitutive law of the material is bilinear. The increment of damage due to fatigue is added to that due to the static increment of the opening displacement and is defined in three different ways. A mathematical cylinder model is employed to carry out a reasonable comparison of the computational efficiency of the three fatigue degradation strategies.
Fatigue degradation strategy 1 assumes that the fatigue damage rate is a function of the opening displacement δ and of the current value of the total damage D. In fatigue degradation strategy 2 the damage directly affects the critical value of the opening displacement δc. Finally, fatigue degradation strategy 3 assumes that the fatigue damage rate is a function only of the opening displacement δ.
Then the three laws are implemented in a code using bond based peridynamics to simulate fatigue crack propagation, and the sensitivity of the results to variations of the main discretization parameters is assessed. Both the cylinder model and the BBPD code provide the same assessment of the three fatigue degradation strategies. Two main conclusions can be obtained by the current work: 1. The third fatigue degradation strategy is the best, among those investigated, for being used in BBPD codes.
2. The cylinder model appears to provide reliable indications about the computational performance of the fatigue laws.

ACKNOWLEDGMENT
The financial support of the University of Padova BIRD2017 NR.175705/17 is gratefully acknowledged.