A cohesive XFEM model for simulating fatigue crack growth under mixed-mode loading and overloading

Structures are subjected to cyclic loads that can vary in direction and magnitude, causing constant amplitude mode I simulations to be too simplistic. This study presents a new approach for fatigue crack propagation in ductile materials that can capture mixed-mode loading and overloading. The extended finite element method is used to deal with arbitrary crack paths. Furthermore, adaptive meshing is applied to minimize computation time. A fracture process zone ahead of the physical crack tip is represented by means of cohesive tractions from which the energy release rate, and thus the stress intensity factor can be extracted for an elastic-plastic material. The approach is therefore compatible with the Paris equation, which is an empirical relation to compute the fatigue crack growth rate. Two different models to compute the cohesive tractions are compared. First, a cohesive zone model with a static cohesive law is used. The second model is based on the interfacial thick level set method in which tractions follow from a given damage profile. Both models show good agreement with a mode I analytical relation and a mixed-mode experiment. Furthermore, it is shown that the presented models can capture crack growth retardation as a result of an overload.


INTRODUCTION
Numerous structures, such as wind turbines, bridges, and cars, are exposed to cyclic loading. The fatigue crack propagation behavior of these structures is commonly simulated with models that are only valid for mode I constant amplitude loading. However, in real-life applications, the applied loads can vary in direction and magnitude. For example, a change in loading direction may cause a mixed-mode stress field around the crack tip, changing the crack growth direction. 1 Furthermore, a change in loading magnitude, such as an overload, can create a significant crack growth retardation effect by means of plasticity induced crack closure. 2,3 This study presents a new fatigue crack growth model that can capture the effects of mixed-mode loading and overloading.
Traditionally, fatigue crack propagation is predicted by means of the Paris equation, 4 which links the cyclic change in stress intensity factor (SIF) to a crack growth rate. A drawback of this method is that the SIF is only a valid measure for materials that show small scale yielding around the crack tip, otherwise linear elastic fracture mechanics (LEFM) cannot be used. A significant amount of yielding can be found for instance when considering overloading in ductile materials, making the SIF lose its validity. Furthermore, Paris parameters are generally determined for mode I loading and thus cannot readily be used for mixed-mode loading. Modifications to the Paris equation exist to extend its applicability to mixed-mode loading 5 and overloading. 6,7 However, these methods are only valid for specific load cases.
To deal with mixed-mode scenarios for which the crack growth direction is not known a priori, the extended finite element method (XFEM) has been used in fatigue crack propagation models. 8,9 In these models, crack tip enrichments following from LEFM are applied to capture the strain field around the crack tip. Furthermore, the J-integral 10 is used to calculate the SIF, which is coupled to the crack growth rate by means of a modified Paris equation. 5 However, the LEFM crack tip enrichments and the J-integral are invalid for elastic-plastic materials. Therefore, the models presented in the works of Singh et al 8 and Pathak et al 9 cannot readily be extended by simply using an elastic-plastic material around the crack to capture the effects of overloading. The XFEM can also be used with a cohesive zone (CZ) approach instead of with crack tip enrichment. 11,12 This is more suitable in combination with plasticity, 13 but has not yet been used for fatigue problems.
The CZ models were originally developed for quasi-static crack growth predictions 14,15 and have been successfully used in the context of interface elements for cases where the crack path is known a priori. [16][17][18] In a CZ model, the loss of bonding strength in the fracture process zone (FPZ) ahead of the macro crack tip is captured by means of tractions. These tractions are calculated with a traction-separation or cohesive law. However, these CZ models cannot readily be used to model crack growth due to cyclic loading. In fatigue crack growth cases, the maximum applied load is smaller than required for quasi-static crack propagation. Consequently, under constant amplitude cyclic loading, the crack will simply not grow. This problem can be solved by means of two different approaches.
In the first approach, the crack tip is forced forward and the amount of cycles required for this jump is computed by means of the Paris equation. In these models, the energy release rate (ERR) is computed by the J-integral around the zero thickness interface elements 19 or by using the local ERR extracted from the cohesive law of the most damaged point in the FPZ. [20][21][22] Since the bulk material is not considered in this J-integral, it can also be used to compute the ERR as a local driving force for propagation of a crack embedded in an elastic-plastic bulk material.
In the second approach, dissipative mechanisms are added to the traction separation law of the CZ model such that the crack can keep on growing in fatigue loading. [23][24][25][26][27][28] This circumvents the use of the Paris equation and its modifications for different loading scenarios. In some of these modified CZ models, the bulk plastic behavior is separated from the creation of crack faces in which case it was shown that an overload causes the plastic zone of the bulk material around the crack tip to increase in size. 24,26,27 This increase in plastic zone size effectively closes the crack before the minimum load level is reached and consequently causes crack growth retardation, as is observed in experiments. The downside of these modified CZ models is that the dissipative equations lack physical meaning. Recently, the interfacial thick level set (ITLS) method was introduced as an alternative for CZ models. [29][30][31] The ITLS method is derived from the thick level set method, 32,33 which is a regularized continuum damage model. The ITLS method itself is similar to a CZ model in the sense that it provides a traction-separation relation. The main difference between the ITLS method and conventional CZ models is that the damage profile over the FPZ for the ITLS method is given instead of being computed from a cohesive law. The shape of the damage profile and the length, which is equal to the length of the FPZ, should be known or estimated a priori in the ITLS method. The method allows for straightforward evaluation of the ERR, which makes the ITLS particularly suitable for combining it with the Paris equation.
In this study, a cohesive XFEM model that can capture fatigue crack growth in arbitrary direction under mixed-mode loading and overloading is presented. Two different models to define cohesive tractions are compared. The first model is a CZ model that uses a static cohesive law as starting point. The CZ model is implemented with the two different ERR extraction methods given by Bak et al 19 and by Harper and Hallet. 20 The second model uses the ITLS method where the ERR extraction method follows from the work of Voormeeren et al. 30 Unlike current ITLS models, the length of the FPZ is not a user input, but will follow automatically from the simulation. For both models, the Paris equation is used for relating the SIF, and thus the ERR, to crack growth. By separating the plastic bulk material from the creation of crack faces and by including a mixed-mode description of the FPZ, Paris' equation does not need any adjustments for overloading or mixed-mode loading. Furthermore, by using XFEM, the crack can continuously change direction, depending on the stress field around the crack tip. Finally, adaptive meshing is used to minimize computation time.
This paper starts with the numerical framework followed by a description of the behavior of the bulk material and the FPZ. The FPZ behavior for both the CZ and ITLS models is presented together with their ERR extraction methods.
Subsequently, the crack tip propagation criteria are discussed. Finally, three numerical examples are presented to verify, validate, and compare the accuracy of the models.

NUMERICAL FRAMEWORK
A cracked medium, as shown in Figure 1, is considered. The specimen is subjected to cyclic loading, causing the crack to propagate. The crack can grow in an arbitrary direction depending on the specimen geometry and the applied loading.
Physically, an FPZ is present ahead of the macro or physical crack tip, which is indicated as the black area in Figure 2. For the numerical models, it is assumed that this zone can be compressed into a single line. The process zone is then represented numerically by means of traction forces across a displacement discontinuity, shown at the bottom of Figure 2. As a result, in the numerical model, the crack tip is located ahead of the FPZ, which is called the fictitious or numerical crack tip. Material ahead of the numerical crack tip does not have any damage. On the other hand, the material in the wake of the physical crack tip has maximum damage indicating traction-free separation of the crack faces.
In order to accurately and efficiently capture fatigue crack growth under general load conditions, the numerical framework needs to have two important characteristics. First of all, the crack path is not known a priori, which means that the crack needs to have the freedom to grow in any arbitrary direction. The standard finite element method will not allow this because displacement jumps can only be present along element boundaries. For this reason, this study uses XFEM, which enables a crack to grow through the elements.
The second characteristic is also related to the unknown crack path. The FPZ for fatigue crack growth in metals is very local. Therefore, a high mesh density is required in the area around the crack tip to correctly capture the strain field. However, as the crack path is unknown, it means that the whole domain should consist of small elements, increasing the computation time significantly. In this study, an adaptive meshing strategy is used, which ensures a high mesh density around the FPZ and a low mesh density far away from this zone to minimize the computation time without significant loss of accuracy.

Phantom node technique
The XFEM was initially developed by Belytschko and Black 34 and Moës et al 35 as a method in which a discontinuity like a crack does not necessarily have to be located along an element boundary as is the case for the standard finite element method. Node enrichments are used to capture a discontinuity in an element without the need of remeshing. Therefore, XFEM can be used to predict crack propagation for cases where the crack growth direction is not known a priori. 11,36,37 For ease of implementation, an alternative formulation of XFEM has been used in this study, which is called the phantom-node method. 12,[38][39][40][41] The phantom node method is visualized in Figure 3. Four phantom nodes are introduced if a crack crosses an element, which initially have the same location as the original regular nodes. With these eight nodes in total, the single element is split into two elements. The first element has regular nodes above the crack and phantom nodes below, and the second element the other way around. The gray areas indicate which part of the element is active, which means that the internal force and stiffness matrix contribution of the first element are evaluated by integrating over the element subdomain above the crack line. A displacement jump between the two elements can be present because the two elements do not share any nodes.

Adaptive meshing
In this study, four noded quadrilateral elements are considered, enabling the use of a quad-tree algorithm for the adaptive meshing process. [42][43][44] A simple visualization of a quad-tree algorithm is shown in Figure 4. It starts with the original element on the left, which is defined as level zero. If refinement of the element is required, it is divided into four new elements, all having the same aspect ratio as the original one. These new elements are defined to be one refinement level higher than the element from which they were created, which makes them level one. This process can then be repeated until the required refinement has been achieved.
The application of the quad-tree algorithm results in nodes that end up in the middle of an element boundary, instead of only at element corners. These nodes are called hanging nodes. A hanging node is constrained to the adjacent corner nodes to maintain a continuous displacement field. The maximum allowable difference in refinement level between adjacent elements is one to achieve better efficiency, following the work of Palle and Dantzig. 42 Refinement is performed when an element is within a certain distance from the numerical crack tip. De-refinement is applied when the physical crack tip has run past a refined element and is a certain distance away from it. The extent of the refinement zone and the minimal required element size are thus problem specific and based on user experience. Other more general refinement criteria that use a finite element error estimation are treated in the works of Palle and Dantzig 42

Regular node
Hanging node FIGURE 4 Mesh refinement by means of a quad-tree algorithm and Tabarraei and Sukumar. 44 It should be noted that elements crossed by a crack are not de-refined to keep the geometry of the crack. Consequently, the total number of elements increases for increasing crack length.

MATERIAL BEHAVIOR
The bulk material behavior is considered separately from the process zone behavior, having both their own set of constitutive relations. As a result, the crack growth process is driven by a combination of the two. First, the bulk material behavior is discussed. Second, the FPZ behavior for both the CZ model and the ITLS method is treated.

Bulk material
An elastic-plastic bulk material is considered. The elastic behavior follows from Hooke's law. The plasticity model uses both isotropic and kinematic hardening to capture plastic flow under cyclic loading. The Von Mises criterion is used to describe the yield surface where s and are the deviatoric stresses and deviatoric back-stresses, respectively. The yield stress y is defined with a nonlinear isotropic hardening rule 45 where 0 is the initial yield stress, Q ∞ the limit value for the yield stress increase, and b is a measure for the rate of change of the yield surface. The increment in equivalent plastic strain d̄p is given as follows: where d p is the plastic strain increment. The kinematic hardening rule is given by 46 where C i is the linear kinematic hardening coefficient and i the nonlinear one. The phantom node technique and the adaptive meshing algorithm require mapping of the plasticity history variables. This mapping is performed based on the distance between the new and old integration points, following the work of Wells et al. 47

Cohesive zone model
The CZ model uses a static cohesive law for each integration point on the crack, from which the traction, and thus the amount of damage can be computed. A simple bilinear cohesive law is used, which is given in Figure 5A. The cohesive law uses a mixed-mode formulation in which the effective traction is a function of the effective displacement jump , which is given by 48 where n and s are the displacement jumps in the normal and tangential direction, respectively. Furthermore, Macaulay brackets are denoted by ⟨·⟩. Parameter is equal to the ratio between the tensile and shear strength of the material. The area below the curve of the cohesive law is equal to G c , which is the critical ERR for mode I or cleavage failure. Parameter max is the maximum effective traction stress and is equal to the tensile strength of the material. Parameter m is defined as the maximum effective displacement jump that has been reached for a specific material point, which has a corresponding damage d that is defined as follows: where it can be seen that the damage d increases for increasing m . Here, i and f are the fracture initiation and final failure displacement jumps, respectively, which follow from max , G c , and K c . Damage only starts to accumulate once the effective displacement jump becomes larger than i . The damage reaches its maximum value of one when f has been reached. The initial cohesive stiffness K c is adjusted for the accumulated damage by multiplying it with (1 − d). The effective traction-displacement relations can therefore be written as follows: From the effective traction, the tractions in the normal direction n and sliding direction s can be calculated using the following equations: Note that, for a negative displacement jump, the normal traction is calculated with the initial cohesive stiffness, irrespective of the amount of accumulated damage. The large initial stiffness ensures contact between the two crack faces by preventing any large negative normal displacement jumps. Therefore, K c should be given a sufficiently large value. Under that condition, its exact value does not influence the global response. The maximum local ERR of a material point having a certain combination of m and d can be calculated as follows: which is equal to the gray area in Figure 5A. Upon unloading and reloading, the local ERR of a single material point is given by which comes from the consideration that G SP is a quadratic function of the applied load and is equal to the gray area in Figure 5B. Notice how the outer static curve moves inwards during unloading and reloading. Harper and Hallett 20 showed that the ERR of a specimen is equal to the local ERR G sp of a single material point at the physical crack tip. Instead of using the single point (SP) estimation, the ERR can also be extracted by computing the J-integral around the interface elements as done by Bak et al. 19 This is equal to integrating the traction-separation behavior of the complete FPZ. The J-integral over the FPZ can be defined as follows: which is computed by means of the Riemann sum.

Interfacial thick level set
The main difference between the ITLS method and the CZ model is in the definition of the damage. For the ITLS method, the damage profile over the FPZ should be provided, instead of having a cohesive law from which the damage is calculated.
An example of such a damage profile is given in Figure 6.
Here, Equation (7) is combined with the damage distribution over the FPZ from the work of Voormeeren et al, 30 which is given by in which has a value of zero at the physical crack tip and is equal to the FPZ length at the numerical crack tip. Furthermore, c d determines the slope of the function and l p is the FPZ length. Only c d needs to be given as an input because in the present XFEM formulation l p automatically follows from the simulation. The ERR is determined globally from integrating the change of the local interfacial free energy over the FPZ for an increase in crack length, which gives the following relation 30 : Here, d ′ is the derivative of the damage with respect to the level set function . The effective displacement jump is the same as given in Equation (5), which means that the ITLS method can be used for mixed-mode loading as well. It should be noted that, for a given damage distribution and opening profile, G IE is equivalent to G J from Equation (12), but written in a form more suitable for the ITLS method.

Fatigue crack relation
The Paris equation 4 is used to compute the crack growth rate da∕dN, which is defined as follows: where a is the crack length, N is the number of cycles, ΔK is the SIF range, c is the Paris constant, and m is the Paris exponent. Both the CZ model and the ITLS method compute the ERR for a given state of the model, and not the SIF. However, for elastic materials, the SIF can simply be computed from the ERR with the following equation: where G is the ERR and E is Young's modulus. Note that the aforementioned equation is only valid for plane stress, but a similar relation exists for plane strain. The effect of different mode-mixities is captured in G by means of the effective displacement jump given in Equation (5). Therefore, even for mixed-mode loading, mode I Paris parameters are used to determine the fatigue crack growth rate.
It should be noted that Equation (16) cannot give a real SIF when considering an elastic-plastic material. However, the SIF obtained from the FPZ inside the plastic bulk can still be used as an indication of the magnitude of the crack driving force and will therefore be considered as an effective K.

CRACK TIP PROPAGATION
The physical and the numerical crack tip both have their own propagation criterion similar to the work of Yang et al. 23 As a result, the size of the FPZ l p is variable because it is defined by the positions of the physical and the numerical crack tip. The CZ model will determine the damage over the whole process zone by means of its traction-separation law. For the ITLS model, the damage distribution is therefore recalculated every time l p is changed by means of Equation (13). In both cases, irreversibility of damage growth is ensured.
In this study, both the physical and numerical crack tip are only allowed to grow from element boundary to element boundary. Figure 7 shows an illustration of a crack with an FPZ in a finite element framework in which the numerical and physical crack tip are indicated. The FPZ is indicated with the solid line of which the arc length is equal to l p . The dashed line indicates the part where the crack faces have been completely separated.

Numerical crack tip
Numerical crack tip propagation, illustrated in Figure 8, is the same for both quasi-static and fatigue loading. It occurs when the maximum principal stress in an integration point located in a small region around the numerical crack tip exceeds the failure stress, which has the same value as max in Figure 5A. This region is indicated in Figure 8 with a circle, which generally has a radius of six times the element size.
The direction of crack growth is determined by means of the maximum principal stress criterion, 50 which is valid for low and moderately high mode mixities, 51 using a nonlocal approach. 11,12 The nonlocal stress is calculated by taking the stress states of the integration points in the neighborhood of the numerical crack tip and weighing them with the following function: Here, r w is the distance from the numerical crack tip to the specific integration point and l w defines the rate of decay of the weight function, generally taken equal to three times a typical element size. The same region as for the maximum principle stress criterion is considered when computing the crack growth direction. The crack extends to the next element boundary once the numerical crack tip propagation criterion is met. As a result, the FPZ size will increase, as shown in Figure 8. After numerical crack tip propagation, a new equilibrium solution is found for the updated discretization. The increased process zone size results in smaller stresses ahead of the numerical crack tip. Numerical crack tip propagation is repeated until max is no longer exceeded at the crack tip.
Generalization of the element-by-element crack growth algorithm to a 3D situation is possible as well following the work of Moës et al 52 and Gravouil et al. 53

Physical crack tip
In simulations for static loading, the physical crack tip moves forward when the local cohesive damage reaches its critical value of one. This means there is no need to keep track of the physical crack tip location. In the case of fatigue crack growth, the physical crack tip for the CZ model does not move forward without help because a static damage of one will not be reached. In fatigue crack growth casesm, the maximum applied load is smaller than required for quasi-static crack propagation. Furthermore, the physical crack tip for the ITLS method does, by definition, not move forward by itself.
For both models, this is solved through forcing the physical crack tip forward by setting the damage of the process zone element that is adjacent to the current physical crack tip equal to one. 22 As a result, the FPZ decreases in size, as is shown in Figure 9. Consequently, the stresses ahead of the numerical crack tip increase such that numerical crack tip propagation could occur again. Physical crack tip propagation is done after every single simulated loading cycle, indicated with the dots in Figure 10.
During postprocessing, the real amount of fatigue cycles dN that a single simulated loading cycle represents is calculated by means of Equation (15). The increase in crack growth length da is known, which is the distance the physical crack tip shifted forward. The ERR is extracted at every time step and, thus, the SIF range ΔK can be computed for every cycle with Equation (16).

NUMERICAL EXAMPLES
In this section, three different numerical examples are presented. Firstly, a mode I linear elastic case is considered in which the presented models are compared against an analytical solution. An overload is included to investigate the accuracy of the different ERR extraction methods for a scenario that does not have constant amplitude loading. Secondly, a mixed-mode linear elastic case is shown where both models are compared against an experiment. Finally, an elastic-plastic mode I case is considered, which is subjected to both a constant amplitude loading and an overload to demonstrate the suitability of the method for modeling crack retardation.
For all cases, the material is aluminium 2024-T4, of which the properties are given in Table 1. The plasticity parameters are taken from the work of Abdollahi and Chakherlou 54 and the Paris parameters from the work of Jeong. 55 The parameter is based on the quasi-static failure stress taken from the work of Davis 56 and the critical ERR is taken from the work of Dursun and Soutis. 57 The cohesive stiffness, which is a numerical parameter, should be chosen sufficiently stiff. All the given numerical examples are in a state of plane-stress.

Mode I linear elastic
The centre crack specimen (CCS) given in Figure 11 is considered for the mode I linear elastic case. Here, the half crack length a has an initial value of 5 mm. The applied cyclic stress has a maximum value of 100 MPa and the load ratio R is equal to 0.1. The main point of interest is the SIF. If both models can correctly capture the SIF, then their fatigue behavior, for which the Paris equation is used, should be accurate as well. In this numerical example, an elastic bulk material is considered.
The models are compared against an analytical relation from LEFM, 58 indicated with "Ana" in the upcoming figures, which is given by where W is the width of the specimen.  The SIF versus crack length for the constant amplitude case is given in Figure 12A. All three lines for K max and K min are on top of each other, showing that both the CZ model and the ITLS method are accurate. Here, the ERR, and thus the SIF for the CZ model, is calculated using the SP approach.
The overload case is similar to the constant amplitude example except for the application of an overload of 1.5 times the constant amplitude load applied at a crack length of 8 mm. The results are given in Figure 12B. Both models and the analytical relation show an increase in the SIF when the overload is applied. After the overload, the SIF is expected to go back to the constant amplitude line, as is seen for the analytical relation and the ITLS method. However, the CZ model shows different behavior with a temporary increase in SIF after the application of the overload.
The cause for this discrepancy is visualized in Figure 13, which shows the traction-separation law at the physical crack tip, where the area below the curve is defined as G SP in the SP method. The second curve gives the traction versus displacement jump of the whole FPZ, the area below which is the physically correct ERR G J . The SP method assumes that both curves are identical.
Before the application of the overload, both curves are exactly on top of each other (see Figure 13A). The SP method is correct because all the material points in the FPZ are on the static envelope of the cohesive law when considering the maximum applied load. All material points are on the traction-separation law of the physical crack tip because, there, the displacement jump is largest. Therefore, only the material point at the physical crack tip has to be considered for determining the ERR. This also holds true for unloading, in which all the material points enter their unloading/reloading branch at the same time. During unloading, all the material points in the FPZ are in the same unloading/reloading state of the traction-separation law as the physical crack tip point. Thus, all the material points will lie on the inner contour, as shown in Figure 5B. During the overload, the entire CZ remains in a critical state and the SP approach still works. After returning to the original load level, all points are unloading. However, when the physical crack tip is then shifted forward, points close to the physical crack tip start to experience an increase in damage, whereas other points remain in the unloading stage. Consequently, the traction-separation law of the physical crack tip no longer represents the whole FPZ, which can be seen in Figure 13B. Therefore, the SP method looses its accuracy when overloads are considered. As soon as the distribution of damage is back to the constant amplitude state, which is reached after multiple physical crack tip propagations, the SP method is working correctly again.
As explained in Section 3.2, the J-integral can also be used to compute the ERR for the CZ model. The J-integral is the area under the traction-separation curve of the whole FPZ in Figure 13. Figure 14 shows that the J-integral does work correctly in the presence of an overload.

Mixed-mode linear elastic
In the second numerical example, the CZ and ITLS models are compared with a mixed-mode experiment done by Jeong. 55 The specimen with geometry, as given in Figure 15, is subjected to a uni-axial constant cyclic loading with a maximum applied stress of 110 MPa and a load ratio of 0.1. The inclined initial crack leads to a mixed-mode stress field at the crack tip, causing the crack to grow along a nontrivial path. The initial crack angle is equal to 22.2 • and its initial length a 0 is equal to 4.22 mm. The bulk material behavior is modeled as purely elastic.  Figure 16 shows that there is a good agreement between the experiment and the models in term of crack growth rate and crack path. Following the work of Jeong, the crack length is defined as the shortest distance between the physical crack tip and the circumference of the hole. The J-integral ERR and the SP method result, in this case, in the same curve because of the constant amplitude loading. The mode I Paris parameters are used, which are fitted for a load ratio of 0.1. With these effective properties and constant amplitude loading, there is no need to include plasticity. It is concluded that the maximum principal stress criterion is appropriate for determining the crack growth direction and that the mixed-mode cohesive law in combination with mode I Paris parameters gives the correct crack growth rate. Figure 17 shows the Von Mises stress field around the crack for the mixed-mode specimen for two different magnifications. The left figure shows that the mesh refinement is only present around the FPZ and close to the fully separated crack faces. The right figure shows that the crack can grow through the elements because of the phantom node technique.

Mode I elastic plastic
Elastic-plastic material behavior is considered for this numerical example. The CCS geometry from Figure 11 is considered again. The maximum applied cyclic stress is equal to 100 MPa and the load ratio R is equal to 0.1. The initial crack length is 5 mm. A constant amplitude and overload case are considered. The overload is applied again at a crack length of 8 mm. As stated in Section 3.4, for an elastic-plastic material, the SIF is considered to be an effective K. However, the Paris parameters c and m are not calibrated for it and the values given in Table 1 are used. The results are therefore qualitative and not quantitative.
The SIFs for a given crack length under constant amplitude loading are given in Figure 18A. The CZ model and the ITLS method both show a lower K max , K min , and ΔK than the elastic solution. Furthermore, due to plasticity, the crack faces are compressed together at the minimum applied stress resulting in an ERR, and thus an effective SIF, of zero. The SP and the J-integral method for the CZ model give the same SIF values for this case. This means that, even for an elastic-plastic material, the SP method works correctly.
It can also be seen that the CZ model and the ITLS method provide different SIFs. This can be explained by inspecting the actual traction distributions in Figure 18B, where the tractions over the FPZ are given for a specific crack length. In addition to having a different traction distribution, the maximum traction for the ITLS is slightly higher than that of the CZ model. Unlike the CZ model, the ITLS is not limited by a maximum stress from a traction-separation law. A larger maximum traction means larger stresses in the neighboring elements resulting in more plasticity and thus a lower SIF. Note that, by changing the damage profile parameter c p and the cohesive stiffness K ITLS c , the traction profile can be changed and the maximum traction could be lowered or increased if required. However, obtaining a specific maximum traction is not as straightforward as for the CZ model.
The SIF and the FPZ size versus crack length for the overload case are given in Figure 19. The SIF for the CZ model is now computed by means of the J-integral. At the application of the overload, the SIF increases in magnitude, leading to an increase in plastic zone size and an increase in length of the FPZ. After the overload, the physical crack tip has to grow through this zone of increased plasticity, which is pushing the crack faces together more severely than before the overload was applied. Consequently, a reduction in SIF is observed as well as a decrease in the FPZ size.
The stress concentration at the numerical crack tip is decreased, due to the increased plasticity and process zone length. As a result, no numerical crack tip propagation occurs. However, after a certain amount of physical crack tip propagation, the stress is large enough again to obtain numerical crack tip propagation. From this point on, the SIF and FPZ size  Figure 20 shows a comparison for both models with and without an overload. It can be seen that the crack growth rate is slowed down after the application of an overload, which is beneficial for the fatigue life. Furthermore, there is a difference between the CZ model and the ITLS method, which is attributed to both models using the same Paris parameters. These parameters should actually be tuned separately for both models such that their model specific plastic behavior is taken into account. However, here, the comparison is purely qualitative and therefore calibration is not performed.

CONCLUSIONS
This study presented a new approach for fatigue crack propagation in ductile materials. The approach is built in a phantom-node framework enabling arbitrary crack growth direction. Adaptive meshing is applied to keep the simulations efficient. The FPZ ahead of the physical crack tip is captured by means of cohesive tractions. As a result, there are two different crack tips, being the macro or physical crack tip and the fictitious or numerical crack tip. The two crack tips have different criteria for propagation. Consequently, there is no need to define an FPZ size, as it automatically follows from the simulation.
The tractions in the FPZ are computed with two different models, being a CZ model and an ITLS model. The main difference between the two is on how the damage is defined. For the CZ model, the damage follows from a traction-separation law, whereas, for the ITLS, the damage profile over the FPZ is predefined.
The two models result in different traction distributions over the FPZ and, therefore, in the case of an elastic-plastic material, a different distribution of plasticity. As a result, the SIFs for the two models are not the same for a given crack length and both models require a different set of Paris parameters. Furthermore, it is observed that it is not straightforward to control the maximum cohesive traction and, thus, the amount of plasticity, for the ITLS method in comparison to the CZ model.
Both models show good agreement with a mode I analytical relation and a mixed-mode experiment. Furthermore, it is shown that the presented models can also capture crack retardation due to an overload when the J-integral is employed for the ERR extraction. However, extracting the ERR from the traction-separation law of the physical crack tip for the CZ model only gives correct results in the case of constant amplitude loading, irrespective of using an elastic or an elastic-plastic material.