## 2 Implementation of the Algorithms of Fatigue and J-Integral in FEM
### Implementation of the Surface Correction Factor  
Implementation of the peridynamic theory in an FEM software package does not allow for the micromodulus, or the corrected micromodulus, to be input directly. The solution is formulated in terms of the stiffness of truss elements representing bonds with characteristics of length, area, and elastic modulus. The micromodulus of bonds near a free surface may be described as the product of \(c\) and a correction factor \(\lambda\) to obtain \(c_{corrected}\). Substituting \(c_{corrected}\) into Equation (6) gives a corrected [PERSON]'s modulus, \(E_{FE,corrected}\), in terms of a corrected micromodulus. The correction factor for each load case can be evaluated via a slight modification to the Volume Method for surface correction. Noting that for a square grid the volume of each collocation point is the same, the correction factor may be written in terms of the number of nodes contained in the family. By evaluating the number of  
Figure 1: Illustration of bonds and horizons for the central node with a crack surface in peridynamics. Not all bonds are drawn for the central node.  
nodes in the family of \(x\) and \(q\) as \(N_{x}\) and \(N_{q}\), and \(N_{0}\) as the number of nodes in the family of a node in the bulk, the correction factor can be written as  
\[\lambda=\frac{2\Delta x^{3}N_{0}}{\Delta x^{3}(N_{x}+N_{q})}=\frac{2N_{0}}{N_{x }+N_{q}} \tag{7}\]  
This procedure does not apply a correction for material volumes bisected by the horizon which have only partial volumes contained inside the horizon. Methods for such volume correction have been established in the book of [PERSON] and [PERSON] [17].  
In this study, a procedure is outlined, as the flow chart shown in Figure 2, by which the remaining life fatigue damage model can be paired directly with an existing FEM software via use of user subroutines. Firstly, the whole domain is divided into FEM region (continuum) and peridynamics region (bonds). Next, the global stiffness matrix is assembled by adding element matrix obtained from continuum mechanics and from peridynamics. Then boundary conditions are applied before solving the linear system. After that, the strain for each bond in peridynamic region is checked to see if any bond will be broken. The updated bond connection will lead to an updated stiffness matrix that will be assembled to the global stiffness so that it will move to the next time step. Interested readers are referred to literature [18; 19] for more details. ABAQUS was used as the platform for this study. The ABAQUS subroutines USDFLD (user-defined field variables) and URDFIL (user-defined subroutine for reading the results.fil file) are capable of satisfying the requirements of the fatigue damage model. The function of the URDFIL subroutine allows access to nodal and element data at the end of every solution step. The USDFLD subroutine allows access to element integration point variables and modification of solution dependent variables during a solution step.  
ABAQUS allows for unique variables to be defined via the solution-dependent variables (SDV) functionality. This method allows variables to be stored in the ABAQUS framework rather than in an external array allowing for reduced computation time and enabling results to be plotted readily. In this study, two SDVs were defined for each element with one retaining the remaining life of the particular element and the other maintaining  
Figure 2: Algorithm for Fatigue Damage Model via User Subroutine.  
element status (intact or severed). Remaining life was initialized as one for all elements in the model in the first step of the analysis, and was updated in each subsequent step. The element status parameter was initially a uniform value of one for all elements in the model. Once the remaining life of a particular element was less than or equal to zero, the status was changed to zero and the element was inactivated (deleted from the solution phase).