Adaptive Finite Element Prediction of Fatigue Life and Crack Path in 2D Structural Components

: The existence of a hole near a growing fatigue crack can cause the crack trajectory to deviate. Unless the hole is too close to the crack, the crack is arrested at the edge of the hole and does not progress further. The purpose of this paper was to predict the crack propagation and lifetime of two ‐ dimension geometries for linear elastic materials in mixed ‐ mode loading using a finite element source code program written in Visual Fortran language. The finite element mesh is generated using the advancing front method. The onset criterion of crack propagation was based on the equivalent stress intensity factor which provides the most important parameter that must be accurately estimated for the mixed ‐ mode loading condition. The maximum circumferential stress theory was used as a direction criterion. The modified compact tension (MCTS) was studied to demonstrate the influence of the hole’s presence on the direction of crack growth and fatigue life for different configurations. The Paris’ law model has been employed to evaluate the mixed ‐ mode fatigue life for MCTS in different configurations under the linear elastic fracture mechanics (LEFMs) assumption. The framework involves a progressive crack extension study of stress intensity factors (SIFs), crack growth direction, and fatigue life estimation. The results show that the fatigue growth was attracted to the hole either changes its direction to reach the hole or floats by the hole and grows as the hole is missed. The results of the study agree with several crack propagation experiments in the literature revealing similar crack propagation trajectory observations.


Introduction
To maintain an acceptable level of durability over the entire operational time, the fatigue crack growth in engineering structures should be monitored. Many experimental methods are used today to study the propagation of cracks in steel plates. Nevertheless, most of the experimental approaches to find an effective solution are destructive, expensive, and time-consuming. To overcome these challenges, scientists have proposed several analytical and finite element techniques to quickly find final solutions [1][2][3]. The structures and components are fractured in many fields of engineering due to the accumulation of fatigue failure. Fatigue has therefore been widely studied [4][5][6]. Several studies have been conducted to date to investigate the interaction between cracks and holes. The finite element method (FEM) research has been shown to offer a reliable justification for experimental dynamic crack propagation in two-dimensional situations. If a crack trajectory is associated with a hole, it has first been found that it avoids the hole with a dynamic load and falls into the hole loaded by a similar quasi-static one [7]. The stress intensity factors (SIFs) describe the intensity of the singular stress fields and are, at the same time, a measure of the amount of the displacements in the crack area, i.e., they are also a measure for crack opening or the relative displacements of the crack surfaces. The SIFs depend on many parameters which are the applied loading of the component, the crack geometry or crack length/depth, the crack position, crack arrangement, the geometry of the component, and the type and location of the load application. Numerical analyses, such as finite element analysis, have become an incredibly useful method for modeling a crack's propagation and determining fracture parameters such as the stress intensity factors and energy release rate. Over the last few decades, improvements to the finite element method have been continuously made in order to improve the handling of complex material models and experimental conditions. In the meantime, new approaches and methods in several areas of study have been suggested and developed rapidly, including the Finite Element Method (FEM), Discrete Element Method (DEM) [8][9][10], Element-Free Galerkin (EFG) method [11], Extended Finite Element Method (XFEM) [12,13], cohesive element method [14] and phase field method [15].
The practical importance of fracture mechanics applications relies on SIFs. Numerical analyses using the FEM usually computes first-order energy derivatives related to crack length or equivalent SIFs to propagate a crack gradually in small, controlled straight-line growth increments. The advantages and disadvantages of using finite elements in computational mechanics have been extensively addressed [16]. De Oliveira Miranda and Meggiolaro et al. study the problem of mesh refining and related errors in computing SIFs using the FEM [16,17]. It was noticed that too much mesh refinement in crack problems could significantly improve the accuracy of the calculation when increasing the computational time. The development of the present software began in 2004 and has continued to add a range of features from version to version for the simulation of two-dimensional fatigue crack growth with the under the LEFM assumptions [18][19][20][21][22][23][24][25]. LEFMs apply when the materialʹs nonlinear deformation is limited to a small area at the tip of the crack, whereas Elastic Plastic Fracture Mechanics (EPFMs) is proposed to analyze the relatively large plastic zones and assumes isotropic and elastic-plastic materials. The energy fields of the strain or the opening displacement near the crack tips are determined based on the assumption. If the energy or opening exceeds the critical value, the crack grows. The fracture is modeled by the splitting node approach [3] and the trajectory follows the successive linear extensions of each crack increment. The propagation process is driven by a linear elastic fracture mechanics (LEFMs) approach with minimum user interaction. In 3D modeling, the re-meshing strategy presents a number of important additional computational challenges compared to a 2D analysis. However, the predicted values of stress intensity factors in both 2D and 3D were almost identical with the same degree of mesh refinement [26]. The developed source code is a finite element software that yields results comparable with the results obtained from the currently available commercial software in fracture analysis. Concerning knowledge, commercial software is not suitable for at least two reasons: first, it does not completely understand the basic algorithm that it uses, and second there is a complete lack of the state-of-the-art programming skills in its development. The efficiency of the proposed program was demonstrated in many cases such as a compact tension specimen under mixed-mode loading [3], specimen with a hole at its center and three interacting cracks [25], modified four-points bend specimen [22], modified compact tension specimen [6], two internal noncolinear cracks [18] and polymethyl methacrylate (PMMA) specimen with three holes [18].

Program Development Procedure
The finite element process usually comprises three main elements-pre-processing, processing, and post-processing. The processing is the most time-consuming part of computing because it involves the evaluation and assembly processes of the stiffness matrices equations and the solver method. The developed code is a simulation program for the assessment of 2D crack propagation processes under conditions of LEFMs. Using the finite element (FE) process, this program predicts the growth of quasi-static fatigue crack in 2D components taking into consideration the mechanical boundary conditions of the fracture. For the crack trajectory simulation, four essential components are used for the adaptive mesh FE analysis, namely, the mesh optimization algorithm, the crack criterion, the direction criterion, and the crack propagation methodology. The mesh refining can be managed by each element's characteristic scale, predicted according to the error estimator. An incremental theory with the Von Mises yield criterion is applied to this initial model. After each load stage is over, the solution errors are calculated. When the error reaches a stated cumulative error at some point in the course of the experiment, the incremental analysis is disrupted and a new FE model is created. The system decides whether to configure the mesh automatically. The system refines the mesh when necessary, within the original boundary conditions. The solution variables (displacement, stresses, strains, etc.) are mapped from the old mesh into the new mesh after it is created. The analysis is then restarted and continues until the errors are again greater than the pre-decided number.
Users may display their responses through a graphical post-processor in the final analysis or at each stage. The crack criterion is employed to assess the start of the cracks. In LEFM, the SIFs are typically used as a criterion for fracture. The directional criteria are the third part of a crack growth simulation using an adaptive mesh in the FE strategy. The direction criterion is also used to determine the direction of the crack's growth. Several methods are used for predicting the direction of a crack's path, including maximum circumferential stress theory, maximum energy release theory, and minimum strain-energy density theory. An FE model is defined at each stage of crack growth. In the first step, the model is provided as an input for the simulation. The output of the algorithm itself is then generated by the models provided in the previous stages. At each step, as the crack grows, the elements within the geometry are deleted and rebuilt with an adaptive strategy and modified for the next propagation process. Figure 1 displays the computational procedure that is used to model fatigue crack growth. The main steps of this procedure are explained in detail by [22] and [23].

Adaptive Mesh Refinements and Crack Growth Criteria
The computational efficiency is highly dependent on numerous parameters, such as the mesh density, mesh refinement, and the number of cracks. During each modeling stage, conventional progressive fracture simulation algorithms update the model geometry with a short straight-line crack forward segment. The direction of each straight, linear growth segment can be determined with a linear elastic 2D model, implying a guideline of growth direction and defining two parameters that characterize a Williams stress field of the first-order near the tip of the crack. These two parameters are the symmetrical first-order mode I and the anti-symmetrical mode II-SIFs (KI, KII). Adaptive remeshing strategies for subsets add each segment of the crack's advance by modifying the existing mesh of FEs. It takes place in the region before the present crack tip to remove elements and nodes and replace them with smaller elements along the crack faces that must be added. Usually, these updated elements provide a small standardized element rosette in the current crack tip position. During each modeling stage, only a geometrical representation of a crack path must be modified because the cracks are not specifically depicted in an adaptive network. This direct method has been extended to represent multiple cracks without refining the meshes between successive crack tips [27]. The maximum circumferential stress theory [28] states that cracks grow in a direction normal to the maximum tangential tensile stress for isotropic materials under mixed-mode loading. The tangential stress in polar coordinates is given by: 3 cos cos sin 2 2 2 2 By solving  for  to get the direction normal to the maximum tangential stress, the solution is given by:   sin 3cos 1 0 which can be solved as: The sign 0  must be reversed to the sign of II K to guarantee the maximum stress associated with the crack increment [29]. The two alternatives are seen in Figure 2.

Stress Intensity Factors Method
The most successful engineering application of fracture mechanics is in the use of crack propagation laws based on stress intensity factor ranges [30]. It was developed by Irwin in the mid-1950s [31]. The stress intensity factor, denoted by K, was defined by Irwin as a mechanical parameter that characterizes the state of stress at a crack point. Numerous analytical solutions for the SIFs have been developed for a variety of common crack configurations. However, these analytical solutions are limited to simple crack geometries and loading conditions [32,33]. The implementation of smallsized yield criteria [34], which established a J-integral approach for studying the nonlinear behavior of the material. For FE analyses, the equivalent domain integral technique replaces the finite-size domain with a divergence theorem by integration along the contour.
The trajectory integral independent contour is as follows: where W is the strain-energy density, i j  are the stresses, ui is the local i axis displacement, s is the contour arc length, nj is the outward unit normal to the contour C around the crack tip (Figure 3a).
(a) (b) For two-dimensional structures, the contour integral i is replaced by the area integral ( Figure  3b). Then, Equation (4) is modified as: For both modes I and II, Equation (5) simplifIES the calculation of mode I and II of SIFs but fails to calculate KI and KII separately for multiple load conditions. The integral to be considered in these situations is the following [35]: where k denotes the index of the crack tip. Two possible methods can be used to calculate the SIFs. The first method is through interactions between the J-integral and the SIFs: The second strategy uses the relation between SIFs and J1 and J2:

Adaptive Mesh Refinement
The optimization scheme in the field of FE mesh is the adaptive mesh refinement. This scheme is reliant on an a posteriori error estimator that is derived from the previous mesh generation. The estimator of the error in the mesh refinement is the stress error norm which is known as relative stress norm error. The h-type adaptive mesh optimization was used to derive the ratio of element norm stress error with the whole area's average standard stress error. The size of the mesh for each element in this scheme is presented as: e e h A  (9) where Ae is the triangle element area. For each element, the norm stress error is represented by: Although, the average norm stress error over the whole domain is represented as: where m is the overall number of elements in the entire domain and *  is the smoothed stress vector.
The integration with the triangular isoparametric dimension will be converted in the FEM by the summation of quadratics according to the Radau rules as follows: and similarly where e t is the thickness of the element for a plane stress case and 1 e t  is the plane strain case.
Consequently, the relative stress norm error e  for each item is considerably less than 5%, which is acceptable for many applications in engineering [23]. Thus, h  (16) where p is the approximation of polynomial order. In this study, P = 2 since the quadratic polynomial is used for the calculation of FEs. Thus, the new element estimated size is: The existing mesh will be considered as the new background mesh and the advancing front method will be replicated depending on the user's set number of mesh refinement.

Fatigue Crack Growth Analysis
As fatigue is a cyclical energy dissipation in the form of hysteretic loops, it is a cyclical damage process. The time elapsed for damage is represented by the number of fatigue cycles (N). This process is defined by crack growth per cycle (da/dN) which depends on the applied stress intensity factor range. In order to achieve fatigue crack growth, the resulting stress intensity range must exceed the threshold stress intensity factor at each crack tip, defined as: where f is a function of geometry and loading and th  is the stress range limit. Equation (18)  indicates that fatigue crack grows will happen otherwise there is no fatigue crack propagation. Table 1 summarizes some of the common models used (proposed by the authors).
Predicting fatigue crack growth by using equivalent SIF ( eq K  ) is still the most commonly used method for dynamic loading structures. Using a modified form of Paris' law, a power law is introduced for the association of the rate of propagation of fatigue crack and the corresponding SIF [36], which is given as: The fatigue life cycles can be estimated with Equation (19) for a crack increment Δa as:

Numerical Results and Discussion
To demonstrate the ability of the developed program to predict fatigue crack growth and fatigue life, an MCTS with different hole positions is used. The dimensions of this geometry with various configurations is shown in Figure 4. The diameter of the third hole is 7 mm, located from the crack tip as seen in Figure 4, at horizontal and vertical distances K and C, respectively. The material properties are elastic modulus E = 205 GPa, Poisson's ratio υ = 0.29, yield strength σy = 285 MPa, and ultimate strength σu = 491 MPa. The Paris' law coefficient C = 8.59 × 10 −14 , Paris' law exponent m = 4.26, and load ratio R = 0.1. Figure 5 shows the initial adaptive finite element mesh for this geometry, as seen in the enlargement of the rosette elements around the crack tip which was refined with a uniform element size in order to capture the severity of the plastic deformation near the crack tip.  For CT01, CT02, CT03, and CT04 in Figures 6-9, respectively, comparisons are made between the results of the present study and the experimental and numerical fatigue crack growth paths conducted by [40,41]. As seen, the predicted crack growth directions in the present study were almost similar to the paths that were predicted experimentally and numerically using the boundary element method (BEM) with BemCracker 2D and finite element method with Quebra2, D (FEM) [40,41]. Additionally, the maximum principal stresses for the four different MCTS configurations are shown in Figure 10.    and numerical results of [24], (c) BemCracker2D [41] and (d) Quebra2D [41]. of [24], (c) BemCracker2D [41] and (d) Quebra2D [41]. The stress intensity factor is the essential criterion for fatigue life assessment. Numerous handbooks can also provide analytical calculations of the SIF for the standard CT geometry. The SIF solution is formulated as follows for the standard CT specimen [40][41][42].
where P is the load applied, t is the thickness of the geometry and f(a/w) refers to the dimensionless stress intensity factor which is denoted as: In this modified specimen, the presence of the hole curved the path of the crack. Because of the curved crack of the path, the solution given in Equation (21) is no longer valid. At this point, we can see the main advantage of mixed-mode crack propagation using numerical methods. In every simulation stage, Mode I SIFs (KI) are extracted from the developed program solutions and substituted through Equation (21) to generate the dimensionless stress intensity factor f(a/w). From the extracted values of KI, a fourth-degree polynomial is fitted for dimensional stress intensity factor indications in CT01, CT02, CT03, and CT04, respectively, as defined in the following equations.
The collected data set can be used in general linear regression methods to create an easy-to-use model to represent the SIF according to the desired parameters of crack length and geometry parameters to help determine the propagation behavior of the crack. These representations of the dimensionless stress intensity factor are among the foremost independent contributions of the present work compared to the previous work [3] besides the variety of the geometries and loading conditions, which confirm the programʹs ability to predict fatigue life and crack growth path reliably. The numerical results for the dimensionless SIF for the MCTS specimen were compared with the analytical solution represented in Equation (23) for the CTS without a hole for the four different configurations of the MCTS as shown in Figures 11-14. The f(a/w) patterns vary from each other, as shown, as the curved crack path is formed. Furthermore, the results of this analysis for the dimensionless SIF f(a/w) were compared with the dimensionless stress factor values calculated by [41] using BEM with BemCracker2D software and the FEM with Quebra2D (FEM) for the four different configurations as seen in Figures 10-14. It can be seen from these figures that the results of this analysis are closely associated with Quera2D rather than BemCracker2D.
In both cases, the predicted crack growth path precisely matches the results of [41]. The propagation of the fatigue crack is always attracted to the hole and can either change its direction and grow through the hole (sink in the holeʹs behavior) or can only be redirected into the hole and propagate as long as it misses (missing the holeʹs behavior). Even if the position of the hole is just slightly changed, the fatigue life cycles of each geometry will vary dramatically depending upon the variations in the path of the crack. It also demonstrates the effectiveness of numerical analysis in predicting such unexpected fatigue crack growth. The simulations showed that the fatigue crack was still attracted to the hole, so it could either curve its path and grow toward the hole, or simply be deflected by the hole as it propagates.
For the first two configurations of the MCTS, i.e., CTS01 and CTS03, the dimensionless SIFs affecting the curves have approximately the same values up to 0.5 of (a/w). Beyond that, the curves begin to diverge depending on the different locations of the third hole. Likewise, the dimensionless SIF value for CTS04 reached the highest value of 10.7 with a crack length of 21 mm directly before the crack sank in the hole, while in CTS01 the f(a/w) value is 8.97 with the same crack length of 21 mm. Therefore, the difference between the two f(a/w) values for CTS03 and CT01 is 1.81 with the same crack length of 21 mm as shown in Figure 15. The difference is caused by the third hole's different position. Therefore, it is concluded that the position of the third hole plays a major role in calculating the dimensionless SIF.     The life of fatigue crack growth is predicted through the mixed-mode equivalent SIF provided by [36] as represented in Equation (20). Comparisons between the present study simulated fatigue life results and the experimental and numerical results using the BemCracker2D (BC2D) and ViDa software [41] for four different MCTS geometries are shown in Figures 16-19 for CT01, CT02, CT03, and CT04, respectively. The simulated fatigue crack growth (FCG) life using the developed program has excellent agreement with the experimental and numerical results of [41] as can be seen in these figures.

Conclusions
In this work, the developed FE source code program has been applied to two-dimensional fatigue crack growth on an MCTS specimen with configurations differentiated by various positions of the third hole relative to the crack tip. The structure and configuration of the geometry play a crucial role in the acquisition of higher values of SIFs in mixed modes, which demonstrated the crack growth trajectory. The presence of a hole on the plate affects the crack and deflects the crack in the direction of the hole, depending on the location of the crack, so that the crack will change or even pass into the hole and grow until the hole is lost. Comparisons with experimental studies show that the developed program can predict, efficiently and economically, the propagation of the crack and fatigue life of arbitrary two-dimensional structural components.