Development of Efficient Finite Element Software of Crack Propagation Simulation using Adaptive Mesh Strategy

The purpose of this study is on the determination of 2D crack paths and surfaces as well as on the evaluation of the stress intensity factors as a part of the damage tolerant assessment. Problem statement: The evaluation of SIFs and crack tip singular stresses for arbitrary fracture structure are a challenging problem, involving the calculation of the crack path and the crack propagation rates at each step especially under mixed mode loading. Approach: This study was provided a finite element code which produces results comparable to the current available commercial software. Throughout the simulation of crack propagation an automatic adaptive mesh was carried out in the vicinity of the crack front nodes and in the elements which represent the higher stresses distribution. The finite element mesh was generated using the advancing front method. The adaptive remising process carried out based on the posteriori stress error norm scheme to obtain an optimal mesh. The onset criterion of crack propagation was based on the stress intensity factors which provide as the most important parameter that must be accurately estimated. Facilitated by the singular elements, the displacement extrapolation technique is employed to calculate the stress intensity factor. Crack direction is predicted using the maximum circumferential stress theory. The fracture was modeled by the splitting node approach and the trajectory follows the successive linear extensions of each crack increment. The propagation process is driven by Linear Elastic Fracture Mechanics (LEFM) approach with minimum user interaction. Results: In evaluating the accuracy of the estimated stress intensity factors and the crack path predictions, the results were compared with sets of experimental data, benchmark analytical solutions as well as numerical results of other researchers. Conclusion/Recommendations: The assessment indicated that the program was highly reliable to evaluate the stress intensity factors and successfully predicts the cracks trajectories. Based on the results, it was recommended to add further development in the software to simulate crack propagation in elastoplastic materials.


INTRODUCTION
For many structures, crack propagation is an important failure mechanism requiring accurate numerical models to implement simulations essential for failure prediction. To perform numerical simulation, the computational methods must be applied to determine fracture response and reliability of cracked structure. Accurate popular method is the Finite Element Method (FEM), which has been extensively used for fracture analysis of cracks [1] .
In general, numerical methods such as the Boundary Element Method (BEM) and the Finite Element Method (FEM) are used in the fracture analysis of structures, because of the complex.shape and continuously changing path of the growing crack. Nowadays, crack propagation problems are usually solved by means of the FEM [2] . The propagation process is divided into several short steps. In each one, SIFs are evaluated and the direction of propagation is calculated using the propagation criterion that is best suited to the particular problem. Finally, the FE model is modified in order to accommodate a short, straight propagation of the crack in the direction evaluated. The whole process is repeated many times and both is timeconsuming and a source of errors if performed manually. This is particularly true for the modification of the mesh. Software capable to handling the solution process in an automatic way is welcome.
Nodal relaxation is frequently used to release nodes, in order to enable the crack tip to propagate through the mesh. In adaptive mesh refinement, most analysts favor either the Delaunay technique or the advancing front method over other techniques when generating meshes due to the quality of unstructured meshes generated [3] .
The present study code has been developed to enable the user to determine 2D-cracks under mixed mode loading with the aid of automatic adaptive mesh finite element. This program is written in FORTRAN language.

The simulation model:
The developed code is a simulation software for the evaluation of 2D crack growth processes under general mixed-mode loading conditions. Based on the Finite Element (FE) method and under consideration of fracture mechanical boundary conditions this software predicts the crack growth rate in 2D components. Throughout the simulation of crack propagation an automatic adaptive mesh is carried out in the vicinity of the crack front nodes and in the elements which represent the higher stresses distribution. The finite element mesh is generated using the advancing front method. In order to suit the requirements of the fracture analysis, the generation of the background mesh and the construction of singular elements have been added to the developed program. The adaptive mesh refinement is employed as the optimization scheme. This scheme bases on a posteriori error estimator which is obtained from the solution from the previous mesh. Here stress error norm is taken as the error estimator. The strategy used to refine the mesh during analysis process is adopted from Alshoaibi et al. [4,5] .

Stress intensity factor and crack growth analysis:
In this study, the displacement extrapolation technique [6] is used to calculate the stress intensity factors as follows: where, E is the modulus of elasticity, ν the Poisson's ratio, κ an elastic parameter equal to 3-4ν for plane strain and (3-ν)(1+ν) for plane stress and L is the length of the element side connected to the crack tip. The near tip nodal displacements at nodes b, c, d and e, shown in Fig. 1 are of a great interest. The displacement tangential and normal to crack plane is denoted as u ' and v' respectively.
In order to simulate crack propagation under linear elastic condition, the crack path direction must be determined. There are several methods use to predict the direction of crack trajectory such as the maximum circumferential stress theory, the maximum energy release rate theory and the minimum strain energy density theory.
The maximum circumferential stress theory asserts that, for isotropic materials under mixed-mode loading, the crack will propagate in a direction normal to maximum tangential tensile stress. In polar coordinates, the tangential stress is given by: The direction normal to the maximum tangential stress can be obtained by solving dσ θ /dθ = 0 for θ. The nontrivial solution is given by: which can be solved as: The automatic crack propagation is characterized by successive propagation steps, performed without user interaction. Each step consists of: • Finite element analysis • Computation of stress intensity factors • Determination of direction of crack extension and new crack tip location • Crack geometry updating • Automatic remeshing Node splitting and relaxation: The opening of the new crack segment, which up to this point was closed, is done by splitting and releasing the nodes along its faces. At this stage, the unbalance is only due to the nonzero nodal forces acting on both sides of the new crack segment. These nodal forces were selfequilibrating as long as the crack segment was closed, but become external forces upon nodal splitting and are no longer in equilibrium with the new boundary conditions. When the condition of crack growth is satisfied at the particular crack tip then the node at the crack tip has to be split into two nodes to simulate the crack opening. Before that the displacement must be updated to the boundary node coordinates if it needs to show the deformation. The length of the splitting is set by the user but must be very small value. The direction of the splitting is taken upwards and downwards at the half of the angle between the segment that firstly contains the current crack tip and the segment that connecting the current crack tip to the approximated next crack tip. This direction is taken in order to avoid interception between the new boundary lines and the updated existence boundary lines.
Numerical results and validation: Double edge cracked plate with two holes: The geometry and the final adaptive mesh of this specimen are shown in Fig. 2. The plate was simply fixed at the bottom edge and subjected to a concentrated load at the upper edge.
The material properties are Young modulus, E = 98 GPa and Poisson's ratio ν = 0.3. Each crack was found to grow towards the nearest hole since it creates a stress drop as seen in Fig. 3a. Once the crack tip has moved beyond the hole, the crack reorients horizontally in mode I loading (Fig. 3b) since they modified the stress distribution at each other's tip as shown in Fig. 3c. Eventually, the cracks were attracted again by the opposite hole and curved towards the holes (Fig. 3d). The results illustrated that the cracks are propagated with the same length since they are symmetric. This proofed [7] statement which said that, two cracks can propagate with the same length if they are symmetric. As can be shown in Fig. 3d, a slight convergence of the crack paths could be detected at the areas close to the holes. The present research simulated the crack propagation for the double edge crack plate simultaneously. Each crack tip had its own refinement with its own concentric mesh that did not depend on each other. The influence of each crack on the other is presented in terms of Von Mises equivalent stress field, as shown in Fig. 4a, which also compared with the Von Mises stress for the same geometry obtained by Bouchard et al. [7] using FORGE2 finite element code as shown in Fig. 4b. Comparison between the present results and those obtained by Bouchard et al. [7] shows that there are a good agreement for crack trajectory and Von Mises stress distribution. Single edge cracked plate with three holes under mixed mode loading: A more complicated mixed mode fracture problem has been studied to demonstrate the performance of the developed program to model the propagation of cracks from initial notches in a single edge cracked plate with three holes and accurately predict the values of stress intensity factors. The geometry of the single edge cracked plate with three holes and its final adaptive mesh are shown in Fig. 5  and 6. The material is polymethylmethacrylate (PMMA beam), which has the Young modulus of E = 300 ksi and Poisson's ratio of ν = 0.3. It is clear that PMMA is not a metal material; however, the results obtained are applicable to any kind of material as long as LEFM is valid. The plate is simply supported near the lower corners and subjected to a concentrated load at the center of the upper edge. The experiments were conducted by [8] with two cases of the initial crack length, a and its location, b, as shown by the table in Fig. 5. This configuration was chosen because two different crack growth trajectories had been predicted by the finite element modeling of Experiment [8] Adaptive FEM [6] (a) (b) Fig. 7: (a): Present results on adaptive finite element meshes and the crack growth trajectory for the single edge cracked plate with three holes and (b): Experimental [8] and FEM [6] trajectory (case 1) the holed specimens, depending on the initial notch geometry. For the first case, the initial crack length, a and its location, b, are 1.5 and 5.0 units, respectively.
The results of the adaptive finite element meshes and the crack growth trajectory are shown in Fig. 7a. Figure 7a shows that the crack growth trajectory passes just above the lower hole and ended at the middle hole. As can be shown in Fig. 7 the crack trajectory still moved towards the hole although the initial crack was far from the hole compared to that in Fig. 7.
For the second case with the initial crack length, a = 1.0 and location b = 4.0 units, respectively, the crack propagates toward the middle hole as shown in Fig. 8a. Experiment [8] Adaptive FEM [6] (a)   [8] and FEM [6] trajectory (case 2) In both cases, the predicted crack growth trajectories very closely resemble the experimental results conducted by [8] and also with the numerical results simulation using adaptive finite element method obtained by [6] as shown in Fig. 7b and 8b.
As expected, there was significant difference in crack trajectories between Fig. 7 and 8. The results proofed that, the crack trajectory was dependent on the initial crack location. An initial location close to the hole could cause the crack path to go towards the hole, while the crack path is far from initial locations remote from the hole.
For case 2 the results of the present study for the crack path prediction is also compared qualitatively with the numerical results using enriched element free Galerkin method (EFGM) [9] and the extended finite element method results (X-FEM) [10] as shown in Fig. 9a and 9b respectively with a good agreement.
The predicted values of the stress intensity factors for mode I and mode II during crack propagation steps are compared with those results observed by [8] for the same geometry and boundary conditions as shown in Fig. 10 and 11. Almost there are good agreements for stress intensity factors values in both modes for case 1.
(a) (b) Fig. 9: Crack propagation path for case 2 (a): EFGM [9] and (b): X-FEM [10] Fig. 10: Comparison among mode I SIF histories for case 1 between present study results and Bittencourt et al. [8] results Fig. 11: Comparison among mode II SIF histories for case 1 between present study results and Bittencourt et al. [8] results However, the values of SIFs differed somewhat when the crack was in the neighborhood of a hole which were shown in Fig. 10 and 11 by the sharp notches decreases in the curves of SIFs.

CONCLUSION
The presented developed research is a powerful tool for the simulation of 2D crack propagation processes. The methodology developed for simulating automatic crack propagation in 2D structural components has been proven very practical and robust. The developed code is assessed with several test specimens and the outcomes are compared to the similar works either by numerically or experimentally; some of which were found in literatures and recent publications. Moreover, the developed software demonstrated that effective and economical prediction of crack propagation paths as well as the calculation of the mixed mode stress intensity factors for arbitrary 2D structural components.