This paper describes implementation of the finite element method (FEM) to investigate crack growth problems in linear elastic fracture mechanics and the correlation of results with experimental and numerical data. The approach involved using two different software to compute stress intensity factors (SIFs), the crack propagation trajectory, and fatigue life estimation in two and three dimensions. According to the software, crack modeling might be run in various ways. The first is a developed source code program written in the Visual Fortran language, while the second is the widely used ANSYS Mechanical APDL 19.2 software. The fatigue crack propagation trajectory and the corresponding SIFs were predicted using these two software programs. The crack direction was investigated using the maximum circumferential stress theory, and the finite element (FE) analysis for fatigue crack growth was done for both software based on Paris's law. The predicted results in both software demonstrated the influence of holes on the crack growth trajectory and all associated stresses and strains. The study's findings agree with other experimental and numerical crack propagation studies presented in the literature that reveal similar crack propagation trajectory observations.
Keywords
 ANSYS mechanical
 smart crack growth
 stress intensity factors
 Linear Finite Element Method (LEFM)
 source code
 hole influence
The prediction of curvilinear fatigue and fracture is a topic for ongoing research. A mathematical study has provided few curved crack path propagation solutions [1, 2], but experimental results can be accurately reproduced in computer ways. The main objective of fracture mechanics is to determine whether or not a structure will fail based on a crack. Crack analysis should start from field stresses in the crack tip, evaluated by the equivalent stress intensity factor (Keq). By comparing the Keq value to a peculiar characteristic of each material called fracture toughness or threshold stress intensity factor, it is possible to determine whether or not a crack member would fail when subjected to static or fatigue loading, respectively. In general, cracks’ initiation and propagation must be associated with stress intensity factors (SIFs) in a complicated state [3,4,5,6,7]. The prediction of curvilinear fatigue and fractures remains a topic of ongoing research. Few curved crack path propagation solutions were obtained by mathematical analysis, but computational approaches can accurately reproduce experimental results [8,9,10]. Numerical simulations typically calculate the energy derivatives in the first order concerning crack length or the equivalent SIFs. Meanwhile, 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) [11,12,13], element free Galerkin (EFG) method [14], extended finite element method (XFEM) [15,16,17,18], cohesive element method [19, 20], and phasefield method [21]. Nevertheless, as a versatile method to simulate crack propagation, the FEM is still prevalent. In general, the propagation of an internal crack in the engineering practice is essential because of its significant effect on the quality and stability of engineering structures. Thus, it is important to predict the crack propagation path for safety and reliability of engineering structures. Therefore, in many industries, the accurate estimation of the crack path and fatigue life estimation is of primary importance in terms of reliability. In various applications, such as aerospace manufacturing and the aviation industry, experimental studies are necessary for fatigue analysis. However, due to high costs, an accurate computational approach is required for crack propagation analysis to prognosticate the direction of crack growth and fatigue life in both static and dynamic loading conditions [22,23,24]. The failure is related to (a) the presence of flaws such as interfaces and cracks and (b) the nature of loads that fluctuate. Cracks tend to initiate and propagate when subjected to fluctuating loads until the structure does not bear the load contributing to complete failure. These cracks are considered to be fatigue cracks, and the expected fatigue life is one of the most significant parameters to determine the structure's safety. It can be computed by adding the number of loading cycles needed to propagate the fatigue crack and lead to failure. The calculation of the growth rate of the crack is generally based on the relationship within the range of SIFs and the crack geometry. One way to predict SIFs is using computational methods such as the XFEM. The XFEM introduced by Belytschko and Black in 1999 [25] was widely used in novel investigations. It is based on the conventional finite element (FE) framework. It uses a novel displacement feature allowing discontinuities to happen, to defeat the need to remesh continuously during the crack tip development process [26]. Extended work was offered to promote practical models to estimate the fatigue failure and fatigue growth to subdue fatigue failures. Several experimental tests are proposed [27, 28], but the performance of the experimental procedures is usually timeconsuming and costly. An effective way to reduce the laboratory work, time, and costs is to incorporate a simulation methodology to including the numerical analysis. Many fatigue crack problems identified to date by the literature use different computational approaches in simulation for simplistic and complicated geometries in 2D and 3D [10, 22, 29,30,31,32,33,34,35]. In this paper, the ANSYS APDL 19.2 is used to precisely predict the mixedmode stress intensity and correlated fatigue response for seven specimens of the modified compact tension specimen (MCTS). In particular, three methods were widely used to illustrate the fatigue estimation of different specimens: (a) the fracture mechanics based on Paris and Erdogan [36]; (b) the method of strainlife independently proposed by Coffin [37]; and (c) the method of stresslife proposed by Wöhler [38]. The first approach was employed in this study for predicting fatigue life by which the crack tip can be described separately by the SIFs. The primary motivation for this study was to perform a proper comparison of the predicted results for crack propagation parameters using the author's 2D software and the 3D ANSYS software as an alternate tool for simulating fatigue crack growth subjected to mixedmode loading. The computational time to simulate the crack propagation in the developed 2D program is approximately onefifth compared with the ANSYS modeling for the same geometry and mesh density.
There are three types of cracks that ANSYS can model: arbitrary, semielliptical, and premeshed. The Separating Morphing and Adaptive Remeshing Technology (SMART) uses the crack front in the premeshed crack approach, and the stress intensity factor is the failure criteria. The rendered node sets will distribute to the front of the crack, as well as to the top and bottom of the crack. By employing SIFs as the failure criteria, the front of the crack is represented in the premeshed crack process used by the “Smart Crack Growth” simulation tool. SMART automatically updates the mesh from crack geometry modifications based on crack propagation on each solution level, eliminating the need for long preprocessing sessions. The sphere of influence method around the crack tip that goes through the thickness can be used to improve the mesh around the crack tip. The geometric regions to be described are the crack tip, the crack's top surface, and the crack's bottom surface. Instead, each of these regions is associated with a nodeset to be used for analysis. The ANSYS MECHANICAL APDL 19.2 software considers mixedmode loading where the maximum circumferential stress criterion is implemented. The following is the formula for the direction of crack propagation [32, 33]:
In ANSYS Mechanical, the simulation of crack propagation is restricted to region II of the typical fatigue crack growth graph based on Paris’ law that can be represented as:
From Eq. (2), for a crack increment, the number of fatigue life cycles may be predicted as:
The equivalent range for the stress intensity factor formula can be found as follows [39]:
The developed program is a Visual Fortranbased simulation for analyzing two dimensional crack propagation processes under LEFM conditions. This program uses the FE approach to estimate the propagation of quasistatic fatigue cracks in 2D components, taking into account the fracture's mechanical boundary conditions. The adaptive mesh FE method includes four main components for crack trajectory simulation: the mesh optimization technique, the crack criterion, the direction criterion, and the crack propagation strategy. The mesh refinement can be controlled by the characteristic scale of each element, which the error estimator predicts. The step by step analysis is disturbed, and a new FE model is constructed when the error reaches a stated, cumulative error at some point throughout the simulation. Within the initial boundary conditions, the system refines the mesh as needed. After the new mesh is generated, the solution variables (displacement, stresses, strains) are mapped from the previous mesh into the newly created mesh. The analysis is then repeated and continued until the errors exceed the predecided number. SIFs are employed in LEFM as fracture criteria. The maximum circumferential stress theory [40], maximum energy release theory [41], and minimum strainenergy density theory [42] are just a few of the methodologies used to predict a crack's propagation direction. The maximum circumferential stress theory has been used in the developed program as a direction criterion.
The most successful engineering implementation of fracture mechanics is crack growth laws based on stress intensity factor ranges [43]. Irwin [44] defined the stress intensity factor, abbreviated as K, as a mechanical parameter that represents the state of stress at a crack point. For various common crack configurations, many analytical solutions for the SIFs were developed. On the other hand, these analytical solutions are restricted to simple crack geometries and loading conditions [2]. Smallsized yield requirements were implemented [45], which established a Jintegral methodology for investigating the nonlinear behavior of the material [46]. The equivalent domain integral methodology replaces the finitesize domain with a divergence theorem by integration along the contour for FEM. The integral independent contour of the trajectory is as follows [46]:
For 2D structures, the contour integral i is replaced by area integral (Figure 2B). Then, Eq. (5) is modified as:
The step by step procedure of the developed program software is illustrated in Figure 3. The details of these steps are explained in the literature [29, 45, 47,48,49,50]. The adaptive mesh refinement was explained in detail by Bashiri and Alshoaibi [45].
In LEFM simulation, the Jintegral by definition takes into account translational mechanical energy balance in the tip region of the crack along x axis. Eq. (6) simplifies the computation of mode I and II of SIFs but they do not compute
The second method uses the relation between SIFs and
Where
Adaptive mesh refinement is an optimization method used in the field of FE mesh. The a posteriori error estimation used in this scheme is derived from the previous mesh generation. The stress error norm, also known as relative stress norm error, is used to estimate the mesh refinement error. The ratio of element norm stress error to the average standard stress error of the entire region was calculated using htype adaptive mesh optimization. The mesh size for each element in this scheme is provided as:
However, the average norm stress error over the whole domain is represented as:
This means each element with
Depending on the user numbers of refinement, the existing mesh will be considered as the new background mesh and the advancing front method will be replicated. The element type was triangular with six nodes.
As illustrated in Figure 4, a four holes’ rectangular plate of 10 mm diameter and a 6 mm long edge crack is analyzed. The dimensions of the plate are 100 mm × 100 mm × 10 mm. The plate is supported on the bottom and is subjected to a tension force of 10 MPa. The plate is made of aluminium 7075T6, whose properties are shown in Table 1.
Mechanical properties of aluminium 7075T6 [51]
Modulus of elasticity, E  72 GPa 
Poisson's ratio, 
0.33 
Yield strength, 
469 MPa 
Ultimate strength, 
538 MPa 
Fracture toughness, 
3,288.76 MPa

The geometry is fixed in the x and y directions at the bottom face in the developed program, while in Ansys the geometry is fixed in the x, y, and z directions.
The meshed geometry created in Ansys workbench has 14,737 nodes and 7,177 elements, while the meshed geometry generated in the developed program has 3,933 nodes and 1,911 elements, as represented in Figure 5.
The predicted crack growth paths in both software is shown in Figure 6, and these were also in good agreement with that path predicted by Knowles and Sternberg [51] using the fast multipole method (FMM) boundary element method.
The maximum and minimum principal stresses were compared side by side in both programs, as represented in Figure 7. In Figure 7, the crack starts to grow in a straight line as there is no influence of any holes. As it approaches the hole, the path propagates towards the hole. The crack deviated its direction because the first upper hole was not close enough to attract the crack to sink in the hole. It continued to propagate in a straight line up until it reached near the second upper hole, close enough to the crack trajectory so that the crack sinks into this hole. The predicted estimation of the mixedmode stress intensity parameters
Figure 9 shows the comparison for the stress contour distribution in y direction between the developed program results, Ansys results, and the results obtained by Knowles and Sternberg [51] using FMM boundary element method.
A more complex mixedmode crack growth problem was investigated to illustrate the performance of the developed program to represent the propagation of cracks from beginning notches in a three holes single edge cracked plate and precisely prognosticate the values of SIFs. The geometry of three holes single edge cracked plate, and its final adaptive mesh, is presented in Figure 10. This case was investigated for the first time by Ingraffea and Grigoriu [53]. The results revealed that small changes in the crack path's initial position significantly impacted the results. As a benchmark for numerical models, this example is commonly used to assess the system's predictive abilities. As indicated in Table 2, the geometry analyzed has 2L = 508 mm, width = 203.2 mm, and thickness = 12.7 mm, with two alternative configurations depending on the initial crack length and position. The initial meshes for this geometry, generated by the developed program and Ansys before starting the simulation, are shown in Figure 11.
Configurations of the three holes’ single edge cracked plate
1  25.4  152.4 
2  38.1  127 
The length of the crack in the first configurations is a_{o} = 25.4 mm, and it is positioned at a distance of b = 152.4 mm from the midlength of the beam. The second configuration differs from the first by having an initial crack length of a_{o} = 38.1 mm, and it is positioned at a distance of b = 127 mm. The specimen is subjected to a cyclic point load at the upper midspan position with a force of 4,448 kN. The material properties are shown in Table 3.
Mechanical properties of the three holes single edge cracked plate [22]
Modulus of elasticity, 
205 GPa 
Poisson's ratio, 
0.3 
Yield strength, 
516 MPa 
Fracture toughness, 
730 MPa

Threshold SIF, Δ 
80 MPa

Paris’ law coefficient, 
1.2 × 10^{−11} 
Paris law exponent, 
3 
Specimens are constrained to all degrees of freedom from the left pin position, similar to the experimental specimen. In contrast, the correct pin only has a first degree of freedom of translation UX. For the first configuration of specimen number one, the simulated crack propagation using the developed program and Ansys software is shown in Figure 12 with an identical crack growth trajectory. The crack growth trajectory goes between the lower and middle holes in this specimen and approaches the middle hole from the right side. It implies that the shear stress intensity factor (KII) surrounding the holes has significantly increased, forcing the crack stepsizes to be decreased.
The crack trajectory using the developed program and Ansys ends up in excellent agreement with experimental data presented in the literature [47, 48], numerical results using the coupled extended meshfreesmoothed meshfree method presented by Knowles and Sternberg [51], and XFEM results obtained by Zhang and Tabiei [21]. Numerical results obtained by Kanth et al. [14] using a polygonal XFEM with numerical integration as shown in Figure 13 A, B, C, D, E, and F, respectively.
The results of this simulation for dimensionless stress factors by using the developed program and Ansys are shown in Figure 14, which are almost identical to each other. As can be seen in this figure, the SIFs increase in magnitude as the crack approaches the hole, demonstrating the hole's influence on the crack trajectory. For the second configuration of specimen number two, the simulated crack propagation using the developed program and Ansys software is shown in Figure 15 with an almost identical crack growth trajectory. The crack in this specimen propagates above the lower hole to the middle hole's left side.
Figure 16 depicts the direction of crack growth for specimen 2, which was compared to the experimental crack trajectory performed in the literature [47, 48] and the numerical findings utilizing the coupled expanded meshfreesmoothed meshfree technique provided by Knowles and Sternberg [51], as illustrated in Figure 16 A, B, C, and D respectively.
The predicted values of the two modes of SIFs in the developed program and Ansys are depicted in Figure 17 with excellent agreement. As observed in the figure, the lower hole has no influence on the crack growth direction, since the initial crack was far from the lower hole, while the middle hole has a stronger influence, attracting the crack to sink into the hole.
Regarding the simulation results of two and threedimensional analysis using the developed program and Ansys software, a satisfactory agreement was obtained between numerical and experimental crack propagation growth trajectories for the simulated specimens. However, for the first geometry there was a differential between the 2D and 3D
The threedimensional simulation is preferable to the twodimensional simulation to visualize better stresses, strains, deformation, and many other characteristics, despite the twodimensional representation being sufficient.
Configurations of the three holes’ single edge cracked plate
1  25.4  152.4 
2  38.1  127 
Mechanical properties of aluminium 7075T6 [51]
Modulus of elasticity, E  72 GPa 
Poisson's ratio, 
0.33 
Yield strength, 
469 MPa 
Ultimate strength, 
538 MPa 
Fracture toughness, 
3,288.76 MPa

Mechanical properties of the three holes single edge cracked plate [22]
Modulus of elasticity, 
205 GPa 
Poisson's ratio, 
0.3 
Yield strength, 
516 MPa 
Fracture toughness, 
730 MPa

Threshold SIF, Δ 
80 MPa

Paris’ law coefficient, 
1.2 × 10^{−11} 
Paris law exponent, 
3 
Studies on new material: carbon dotgraphene oxidezinc oxide nanocomplex Growth, spectroscopic, electrical and mechanical studies on hexamethylenetetramine succinate crystal: a new third harmonic generation material An insight on spectral, microstructural, electrical and mechanical characterization of ammonium oxalate monohydrate crystals Stereometric analysis of Ta_{2}O_{5} thin films Growth, spectral, optical, dielectric and third order nonlinear studies on 2amino 4methylpyridinium salicylate single crystals