Simulation of Quasi-Static Crack Propagation by Adaptive Finite Element Method

: The ﬁnite element method (FEM) is a widely used technique in research, including but not restricted to the growth of cracks in engineering applications. However, failure to use ﬁne meshes poses problems in modeling the singular stress ﬁeld around the crack tip in the singular element region. This work aims at using the original source code program by Visual FORTRAN language to predict the crack propagation and fatigue lifetime using the adaptive dens mesh ﬁnite element method. This developed program involves the adaptive mesh generator according to the advancing front method as well as both the pre-processing and post-processing for the crack growth simulation under linear elastic fracture mechanics theory. The stress state at a crack tip is characterized by the stress intensity factor associated with the rate of crack growth. The quarter-point singular elements are constructed around the crack tip to accurately represent the singularity of this region. Under linear elastic fracture mechanics (LEFM) with an assumption in various conﬁgurations, the Paris law model was employed to evaluate mixed-mode fatigue life for two specimens under constant amplitude loading. The framework includes a progressive analysis of the stress intensity factors (SIFs), the direction of crack growth, and the estimation of fatigue life. The results of the analysis are consistent with other experimental and numerical studies in the literature for the prediction of the fatigue crack growth trajectories as well as the calculation of stress intensity factors.


Introduction
The finite element method (FEM) is definitely the most common and effective analytical technique for analyzing the behavior of a wide variety of engineering and physical issues. One of the essential uses of FEM is the study of crack propagation. The propagation of the crack reduces components' ability to resist the external load and eventually break the components. Analyzing fatigue crack growth is necessary to ensure the stability of structures subjected to cyclic loading. Cracks begin due to the presence of plastic strain caused by cyclic tension, and they grow due to the tensile stress. However, compressive loads do not lead to fatigue cracks due to the local tensile stress [1]. A variety of software has been developed for general purposes for finite elements, verified and calibrated through the years and now available on request, the most well-known being three-dimensional, such as ANSYS [2], ABAQUS [3], NASTRAN, FRANC3D, and COMSOL. In addition, there are numerous 2D simulation software for crack propagation simulation, e.g., NASGRO, AFGROW, FRANC2D, and FASTRAN. Many researchers have also developed an effective method for estimation of fatigue breakage growth in 2D linear elastic structures with multimode loading [4][5][6][7]. Determining the accurate stress intensity factor of a cracked structure in LEFM is very crucial in accessing the integrity of the crack, especially if the calculation is carried out using the finite element technique with extremely fine mesh. The propagation of the crack can be simulated at the highest accuracy by increasing the mesh density, as well as estimating the stress intensity factors accurately. In addition, very fine mesh around the crack tip is needed for precise prediction of SIFs using a nodal displacement technique such as the displacement extrapolation technique (DET). The DET requires configuration of special elements in the vicinity of the crack tip, by correctly representing the stress field singularity at the crack tip. These special elements, known as singular elements, need to be constructed in a rosette formation around every crack tip. Very small-size elements can be optimally created around the crack tip with the use of an adaptive mesh refinement scheme. Generating overall fine mesh leads to greater computational time. This procedure was reduced by using the adaptive mesh strategy, which increases the mesh only on the required areas. The adaptive mesh refinement scheme is another method to generate the optimal mesh in a very efficient way. Many studies on mesh refinement problems and related errors in computing SIFs using the FEM were conducted [8,9]. Another study [10] was been performed to clarify the effect of the in-plane and out-of-plane constraints on the ductile fracture with different crack sizes, specimen thicknesses, and span lengths. They concluded that the lower in-plane and out-of-plane constraint levels introduce higher fracture properties. It is more challenging to combine the extreme fine mesh generation with the adaptive scheme and solve the stiffness equation matrix. The benefits sought here are both faster execution time and the ability to process larger problems. In order to simulate 2D cracks under mixed mode loading, the current developed software code is formulated to allow the researcher to estimate the fatigue life and crack trajectory using the automated adaptive mesh finite element [11][12][13][14][15]. This software was created in 2004 and continues to include several features for the simulation of two-dimensional fatigue crack growth under LEFM assumptions [12,[16][17][18][19][20][21]. The use of commercial software for engineers is not appropriate in at least two aspects: First, the basic algorithm that lies behind it is not fully understood, and second, the execution is completely apprehended throughout the programming ability. Commercial software can be used to model crack propagation as well, but such software is very expensive and can hardly get the source code to develop it.

Developed Program Framework
The code that was developed is a simulation software to assess the 2D crack propagation process under LEFM conditions. This software predicts the growth of quasi-static crack growth in 2D components using the finite element method, taking into account the mechanical parameters of the fracture. Four essential features for the adaptive mesh finite element (FE) analysis are used for the crack direction simulation, namely, the mesh optimization algorithm, the crack criteria, the criterion of direction, and the methodology of crack propagation. The mesh refining can be controlled by the characteristic scale of each element predicted, based on the current error estimator. An incremental principle with the von Mises yield criterion is applied to this initial model. The solution errors are computed after each load stage is over. The incremental analysis is interrupted when the error exceeds a specified cumulative error at some stage and a new FE plan is generated. The program automatically configures the mesh with a new mesh refinement. After it is generated, the solution variables (displacement, stresses, strains, etc.) are transferred from the old mesh to the new mesh. The analysis is then resumed and progresses until the errors are again higher than the pre-decided amount.
In order to examine the start of the crack growth, the crack growth criterion is employed. The LEFM typically utilizes SIFs as a fracture criterion. Various techniques of estimating the path of a crack are used, such as the maximum circumferential stress theory, theory of maximum energy release, and theory of minimum energy density. At any stage of crack propagation, a FE model is defined. This model is given in the first step as an input for the modeling. The algorithm output is then generated via the models in the subsequent steps. At each stage, as the crack grows, the geometry elements are deleted and reconstructed using an adaptive technique and updated for the next propagation process. Figure 1 demonstrates the simulation procedure used to model quasi-static crack growth. The main steps of this procedure are explained in detail by [11,14]. Metals 2021, 11, x FOR PEER REVIEW 3 of 17

Displacement Extrapolation Technique (DET)
The DET is based on the nodal displacement around the crack tip. The construction of quarter-point elements around the crack tip is generally needed for this procedure. Generally, the existence of the quarter-point element is essential in order to correctly represent the linear elastic singularity (1/ r ) for stresses and strains at the crack tip. The polynomial isoparametrically representative of the singularity is typically obtained by moving the mid-side nodes adjacent to the crack tip to a quarter-length edge closer to the crack tip. Crack tip elements based on this method were separately suggested by [22,23]. In this study, the natural triangle-quarter-point element was selected as the type of cracktip element and its configuration follows the schematic formation of the rosette around the crack-tip, as seen in Figure 2.

Displacement Extrapolation Technique (DET)
The DET is based on the nodal displacement around the crack tip. The construction of quarter-point elements around the crack tip is generally needed for this procedure. Generally, the existence of the quarter-point element is essential in order to correctly represent the linear elastic singularity (1/ √ r) for stresses and strains at the crack tip. The polynomial isoparametrically representative of the singularity is typically obtained by moving the mid-side nodes adjacent to the crack tip to a quarter-length edge closer to the crack tip. Crack tip elements based on this method were separately suggested by [22,23]. In this study, the natural triangle-quarter-point element was selected as the type of crack-tip element and its configuration follows the schematic formation of the rosette around the crack-tip, as seen in Figure 2. For the calculation of stress intensity factors, the displacement extrapolation method [24] was used as follows: where E is the modulus of elasticity, ν is the Poisson's ratio, κ is the elastic parameter and L is the quarter-point element length. The u ′ and v′ are the displacement components in the x' and y' directions, respectively. The subscriptions represent their position, as seen in Figure 2.

Adaptive Mesh Refinement
To minimize expected errors after a finite element solution has been achieved, an adaptive mesh refinement technique is used. The method of adaptive mesh refinement measures the mesh's adequacy and refines the mesh wherever the estimated error is large. Until user-definable error tolerance is reached, the system iterates the mesh refinement and solution. Because the precision of the solution depends on these tolerance limits, it is important for the use of adaptive mesh generators to provide a good understanding of the FEM in an effective manner. The method is referred to as adaptive since at all times the process relies on previous results. The adaptive remeshing method was carried out on the basis of the posteriori stress error standard scheme to achieve the optimum mesh from [16]. The software adopted a frontal solver that is an effective direct solver used to solve a linear equation system. In h-type adaptive mesh refinement, the major point is to obtain the ratio of element normal stress error to the average normal stress error for the entire domain, which is also known as the relative stress norm error, and a new size can be predicted from this ratio for the refinement method. The mesh size is defined in the procedure of each element as: For the calculation of stress intensity factors, the displacement extrapolation method [24] was used as follows: where E is the modulus of elasticity, ν is the Poisson's ratio, κ is the elastic parameter defined by and L is the quarter-point element length. The u and v are the displacement components in the x' and y' directions, respectively. The subscriptions represent their position, as seen in Figure 2.

Adaptive Mesh Refinement
To minimize expected errors after a finite element solution has been achieved, an adaptive mesh refinement technique is used. The method of adaptive mesh refinement measures the mesh's adequacy and refines the mesh wherever the estimated error is large. Until user-definable error tolerance is reached, the system iterates the mesh refinement and solution. Because the precision of the solution depends on these tolerance limits, it is important for the use of adaptive mesh generators to provide a good understanding of the FEM in an effective manner. The method is referred to as adaptive since at all times the process relies on previous results. The adaptive remeshing method was carried out on the basis of the posteriori stress error standard scheme to achieve the optimum mesh from [16]. The software adopted a frontal solver that is an effective direct solver used to solve a linear equation system. In h-type adaptive mesh refinement, the major point is to obtain the ratio of element normal stress error to the average normal stress error for the entire domain, which is also known as the relative stress norm error, and a new size can be predicted from this ratio for the refinement method. The mesh size is defined in the procedure of each element as: where A e is the area of the triangle element. The norm stress error for each element is defined by whereas the average norm stress error for the whole domain is where m indicates the total number of elements in the whole domain and σ * is the smoothed stress vector. In the finite element treatment the integration with the isoparametric triangular element is converted by the summation of quadratures following the Radau rules [25] as follows: where W P is a weighting factor and is J e is the Jacobian matrix. Similarly, where t e is the element thickness for a plane stress condition and t e = 1 for a plane strain condition. Therefore, the relative stress norm error ζ e for each item is considerably less than some identified value [26]. Thus, And the relative stress error level of the new element is defined as permissible error as This implies that any element with ε e > 1 must be optimized and the new mesh size must be predicted. The asymptotic convergence rate criteria are used, which assumed e e ∝ h p e (11) where p is the approximation of the polynomial order. In the analysis, p = 2 is used for the approximation of finite elements as a quadratic polynomial. The predicted sizes of the new element are stated as follows: where h e is the old element size and p is the order of the interpolation shape function.
Convergence of the mesh is dependent of the size of the new element, which defines how many elements in a model are required to ensure that the results of an analysis are not affected by changing the mesh size. System response (stress, deformation) converge with decreasing element size to a repeatable solution. Further refinement of the mesh does not affect results because the model and its results are now independent of the mesh.
The present mesh is known as the new background mesh and the advancing front method is replicated according to the amount of mesh refinements set by the user.
The mesh optimization is used in the final stage of the mesh generation in order to enhance the shape of the elements. The topological structure of the mesh is fixed in the process of mesh smoothing, i.e., the element's nodal connections are not changed, but the inner nodes are repositioned to create triangles with much improved shapes. The most effective computational smoothing algorithm is the well-known Laplacian smoothing [27], which repositions the inner node created by its neighboring nodes at the center of the polygon. The new position of an internal node i is computed as where N n is the number of nodes linked to node i. The mesh smoothing process consists of several iterations.

Crack Growth Analysis
The direction of the crack path under linear elastic conditions must be computed to facilitate crack propagation simulation. The maximum circumferential stress theory states that for isotropic materials under mixed loading mode the crack grows in a direction normal to a maximum tangential tensile stress. The tangential stress is estimated in polar coordinates as The direction normal to the tangential maximum stress can be obtained by resolving dσ θ /dθ = 0 for θ. The nontrivial solution is determined by which can be solved as The sign of θ must be opposite the sign of K I I to ensure the optimal opening stress associated with the crack direction [28]. Figure 3 illustrated the two possibilities of the crack growth direction.
The sign of θ must be opposite the sign of II K to ensure the optimal opening stress associated with the crack direction [28]. Figure 3 illustrated the two possibilities of the crack growth direction. In the case of fatigue crack growth, the resulting stress intensity range at each crack tip must exceed the stress intensity threshold, specified as where f is a geometrical and loading function and Δ th σ is the stress range limit. According to Equation (17) 3 2 cos ( / 2) 3 cos ( / 2)sin( / 2) In the case of fatigue crack growth, the resulting stress intensity range at each crack tip must exceed the stress intensity threshold, specified as where f is a geometrical and loading function and ∆σ th is the stress range limit. According to Equation (17), the crack is not propagated if ∆σ < ∆σ th . This equation was practically modified by using another parameter known as the equivalent stress intensity factor range, ∆K Ieq . Therefore, if ∆K Ieq > ∆K th , this indicates commencement of fatigue crack growth. This parameter is set to ∆K Ieq = ∆K I cos 3 (θ/2) − 3∆K I I cos 2 (θ/2) sin(θ/2) In the modified equation of the Paris law, Tanaka [20] derived an innovative law known as the power law for determining crack growth in response to fatigue with the equivalent stress intensity factor (∆K eq ) as where a is the length of the crack, N is the number of cycles, C is the Paris constant (mm/cycle), and m is the Paris exponent. The total number of fatigue lifecycles can be calculated using Equation (19) for an increase in crack length as

Two Internal Non-Colinear Cracks
For this geometry, there were two internal, parallel, non-colinear, and non-angled cracks in a rectangular specimen with dimensions (90 mm/180 mm). The initial crack length was a = 10 mm for both cracks. As seen in Figure 4a, this geometry was subjected to acyclic tension (σ max = 160 N/mm, σ min = 0) at the upper end and restricted at the bottom side. The distance between the two tips was 15 mm in the horizontal direction and 5 mm in the vertical direction. The adaptive dense mesh is shown in Figure 4b. The selected material was aluminum, which has the material properties shown in Table 1.
cracks in a rectangular specimen with dimensions (90 mm/180 mm). The initial crack length was a = 10 mm for both cracks. As seen in Figure 4a, this geometry was subjected to acyclic tension ( max min 160 N/mm, 0 σ σ = = ) at the upper end and restricted at the bottom side. The distance between the two tips was 15 mm in the horizontal direction and 5 mm in the vertical direction. The adaptive dense mesh is shown in Figure 4b. The selected material was aluminum, which has the material properties shown in Table 1.    This specimen contained four crack tips, which made it interesting to observe the interaction between cracks and to further explore the performance of the developed software in the simulation of multiple cracks.
The predicted crack growth is shown in Figure 5a, which closely resembled the experimental result of Tu and Cai (1993), as illustrated in Figure 5c. These predicted crack growth trajectories were also in agreement with the numerical results obtained by [29] using the linear smoothed extended finite element method, which was compared to the numerical results reported by [5] using a meshless method with enriched weight functions, as shown in Figure 5b.
The predicted crack growth is shown in Figure 5a, which closely resembled the experimental result of Tu and Cai (1993), as illustrated in Figure 5c. These predicted crack growth trajectories were also in agreement with the numerical results obtained by [29] using the linear smoothed extended finite element method, which was compared to the numerical results reported by [5] using a meshless method with enriched weight functions, as shown in Figure 5b.  Both cracks demonstrated in the beginning a pure mode I of approximately the same SIF values. After that, the mode II of the SIF increased at tip A above that of tip B while the second mode of SIFs became negative at A, thus making the crack path curve towards  the other break. Eventually the second mode of the SIFs at A tended to decrease as crack tip A moved closer when the first mode at B increased continually. Finally, the equivalent mode I of the SIF at B exceeded the fracture toughness and unstable fracture occurred at crack tip B. The fatigue life of the structure was evaluated as 6840 cycles, which was in good agreement with the results obtained by [5] using a meshless method, as shown in Figure 7, as well as with the numerical results obtained by [29].  Both cracks demonstrated in the beginning a pure mode I of approximately the same SIF values. After that, the mode II of the SIF increased at tip A above that of tip B while the second mode of SIFs became negative at A, thus making the crack path curve towards the other break. Eventually the second mode of the SIFs at A tended to decrease as crack tip A moved closer when the first mode at B increased continually. Finally, the equivalent mode I of the SIF at B exceeded the fracture toughness and unstable fracture occurred at crack tip B. The fatigue life of the structure was evaluated as 6840 cycles, which was in good agreement with the results obtained by [5] using a meshless method, as shown in Figure 7, as well as with the numerical results obtained by [29].

PMMA Beam Specimen
The PMMA beam geometry offers a benchmark evaluation based on the numerical and experimental work of [31]. The beams were made of polymethyl methacrylate (PMMA), which is a standard material option for crack path investigations as it is relatively homogeneous and exhibits brittle fracture behavior at room temperature. The specimen was under a cyclic point load and acted on the top mid-span position with a value of 4.448 kN. The properties of the materials were taken as modulus of elasticity, E = 205 GPa, yield stress σ y = 516 MPa, threshold stress intensity factor ∆k th = 80 MPa √ mm, ∆K IC = 730 MPa √ mm, Paris law coefficient, C = 1.2 × 10 −11 , Paris law exponent m = 3, and Poisson's ratio ν = 0.3. The thickness of the specimen was 12.7 mm and there were two different configurations depending on the initial crack length (a) and its position (b), as shown in Table 2. The specimen's geometry and the initial adaptive dens mesh are shown in Figure 8. in Figure 8.

Case I
The simulated crack growth for this specimen moved between the bottom and midhole and reached the mid-hole on the right side. It presented a significant increase in the

Case I
The simulated crack growth for this specimen moved between the bottom and midhole and reached the mid-hole on the right side. It presented a significant increase in the KII component of the shear stress intensity factor across the cracks, which forced the stepsizes of the crack to be shortened. The findings of the crack trajectory during propagation were excellently consistent with the experimental results of the crack trajectory [32], the numerical results obtained by [33] using A polygonal extended finite element method (XFEM) with numerical integration for linear elastic fracture mechanics, the XFEM results using ABAQUS software obtained by obtained by [34], and with the numerical results using the coupled extended meshfree-smoothed meshfree method presented by [35], as shown in Figure 9a-e, respectively. The maximum principal stress distribution is shown in Figure 10.
were excellently consistent with the experimental results of the crack trajectory [32], the numerical results obtained by [33] using A polygonal extended finite element method (XFEM) with numerical integration for linear elastic fracture mechanics, the XFEM results using ABAQUS software obtained by obtained by [34], and with the numerical results using the coupled extended meshfree-smoothed meshfree method presented by [35], as shown in Figure 9a-e, respectively. The maximum principal stress distribution is shown in Figure 10. Figure 9. Final crack growth path for case I: (a) present study, (b) experimental results [32], (c) [33] with permission of Elsevier 2019, (d) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020. The results of this simulation were compared with those from XFEM using the smooth nodal stress technique by Peng et al. 2017, as shown in Figure 11, with good agreement. It was found that as the crack approached the hole, the SIFs appeared to change to a greater amplitude. The predicted fatigue life for this specimen was compared to the analytical results calculated by [36] using Paris and Walker models, as shown in Figure 12, with good agreement.  [32], (c) [33] with permission of Elsevier 2019, (d) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020.
sizes of the crack to be shortened. The findings of the crack trajectory during propagation were excellently consistent with the experimental results of the crack trajectory [32], the numerical results obtained by [33] using A polygonal extended finite element method (XFEM) with numerical integration for linear elastic fracture mechanics, the XFEM results using ABAQUS software obtained by obtained by [34], and with the numerical results using the coupled extended meshfree-smoothed meshfree method presented by [35], as shown in Figure 9a-e, respectively. The maximum principal stress distribution is shown in Figure 10. Figure 9. Final crack growth path for case I: (a) present study, (b) experimental results [32], (c) [33] with permission of Elsevier 2019, (d) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020. The results of this simulation were compared with those from XFEM using the smooth nodal stress technique by Peng et al. 2017, as shown in Figure 11, with good agreement. It was found that as the crack approached the hole, the SIFs appeared to change to a greater amplitude. The predicted fatigue life for this specimen was compared to the analytical results calculated by [36] using Paris and Walker models, as shown in Figure 12, with good agreement. The results of this simulation were compared with those from XFEM using the smooth nodal stress technique by Peng et al. 2017, as shown in Figure 11, with good agreement. It was found that as the crack approached the hole, the SIFs appeared to change to a greater amplitude. The predicted fatigue life for this specimen was compared to the analytical results calculated by [36] using Paris and Walker models, as shown in Figure 12, with good agreement.

Case II
According to Table 2, the differences between this case and the previous case were the initial crack length and its position from the mid-span, which were 38.1 mm and 127 mm, respectively. The crack moved above the lower hole in this specimen and stopped at the central hole from the left. The results of the crack trajectory during propagation were excellently close to the experimental results of the crack trajectory obtained by [32], XFEM results using ABAQUS software obtained by [34], as well as the numerical results obtained by [35] using the coupled extended meshfree-smoothed meshfree method, as shown in Figure 13a-d, respectively. The distribution of the von Mises stress is shown in Figure 14.

Case II
According to Table 2, the differences between this case and the previous case were the initial crack length and its position from the mid-span, which were 38.1 mm and 127 mm, respectively. The crack moved above the lower hole in this specimen and stopped at the central hole from the left. The results of the crack trajectory during propagation were

Case II
According to Table 2, the differences between this case and the previous case were the initial crack length and its position from the mid-span, which were 38.1 mm and 127 mm, respectively. The crack moved above the lower hole in this specimen and stopped at the central hole from the left. The results of the crack trajectory during propagation were excellently close to the experimental results of the crack trajectory obtained by [32], XFEM results using ABAQUS software obtained by [34], as well as the numerical results obtained by [35] using the coupled extended meshfree-smoothed meshfree method, as shown in Figure 13a-d, respectively. The distribution of the von Mises stress is shown in Figure 14. Figure 13. Final crack growth path for case I: (a) present study, (b) experimental results [32], (c) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020. The findings of the study for the dimensionless stress factor were compared with those achieved in XFEM with the smooth nodal stress system [37], as seen in Figure 15, with identical results.  [32], (c) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020. excellently close to the experimental results of the crack trajectory obtained by [32], XFEM results using ABAQUS software obtained by [34], as well as the numerical results obtained by [35] using the coupled extended meshfree-smoothed meshfree method, as shown in Figure 13a-d, respectively. The distribution of the von Mises stress is shown in Figure 14.
(a) (b) (c) (d) Figure 13. Final crack growth path for case I: (a) present study, (b) experimental results [32], (c) [34] with permission of Elsevier 2018, and (e) [35] with permission of Elsevier 2020. The findings of the study for the dimensionless stress factor were compared with those achieved in XFEM with the smooth nodal stress system [37], as seen in Figure 15, with identical results. The findings of the study for the dimensionless stress factor were compared with those achieved in XFEM with the smooth nodal stress system [37], as seen in Figure 15, with identical results.

Conclusions
The results of the developed program simulation were compared with experimental and numerical data for the two internal non-colinear cracks and the three-point bending beam with three holes with two different configurations. The developed program combines the adaptive mesh refinement with increasing mesh density in the required area only in order to reduce the computational time while increasing the solution accuracy.

Conclusions
The results of the developed program simulation were compared with experimental and numerical data for the two internal non-colinear cracks and the three-point bending beam with three holes with two different configurations. The developed program combines the adaptive mesh refinement with increasing mesh density in the required area only in order to reduce the computational time while increasing the solution accuracy. The norm stress error is taken as a posterior estimator for the h-type adaptive refinement. With this series of simulations, the capability of the developed program was demonstrated to accurately predict the crack path trajectory, stress intensity factors, and fatigue life under constant amplitude loading. In these simulation sequences, holes act as a crack stopper and attract a crack trajectory to growth. Such findings support that the algorithm can be used to identify crack-stopping holes used in damage tolerance designs.