Finite Element Model of Crack Growth under Mixed Mode Loading

In this paper, in order to predict the crack growth trajectory and to evaluate the SIF under mixed modes (I & II), one proposes a new finite element program for crack growth using the source code written in FORTRAN. The fin ite element mesh is generated using an advancing front method, where the generation of the background mesh and the construction of singular elements are also added to this developed programme to facilitate the crack process and the fracture analysis. Displacement Extrapolation Technique (DET) was employed to evaluate the SIFs under mixed mode loading conditions. Therefore, the accuracy of both SIF`s values and the crack path predictions results are compared and validated with other relevant published research work. However, the assessment indicated that this developed fin ite element programme is reliable and robust to evaluate the SIFs and predicts the crack trajectories successfully based on the applied loading conditions.


Introduction
Structural failu re can be generally associated with one or more fracture of the materials making that structure. In fact, fractu re mechanics deals with the irrevers ible process of rupture due to nucleation or sudden crack and crack gro wth. The influence of cracking on the structural strength has been wid ely ap p reciated s in ce th e end o f th e 19th centu ry. Ho wever, so me aspects o f its natu re and influen ce st ill remain un kno wn [1]. Th e use o f crack p ropagat ion laws based on stress intensity factor range is the most successful engineering applicat ion of fracture mechanics. In the elastic fracture analysis, the stress intensity factors sufficient ly define the stress field close to the crack tip and prov ide fundamental in fo rmat io n o f ho w th e crack is go ing to p rop ag ate. Bas ically , th e estimatio n meth ods can be categ o rized in to two g ro u p s , th os e bas ed o n field extrapolation near the crack t ip and those which make use of the energy release when the crack propagates. The latter group includes the J-contour integration, the virtual crack extension and the strain energy release rate method. The main d isadv antag e o f these methods is th at th e st ress intens ity factor co mponents, K I and K II in mi xed mode prob lems are either imposs ib le o r very d ifficu lt t o be

Finite Element Model
Most established research was on computer base by applying numerical methods which were mostly uses commercial software such as; [2] an enriched partit ion of unity finite element method (PUFEM ), which is known as one of the meshless method to calculate the SIF in LEFM under plane stress and plane strain condition. Another, numerical examp les were performed to evaluate the generated stress intensity factors directly fro m the scaled boundary finite element solution for the singular stress field [3]. More recent work was developed variat ional mesh-free method, to evaluate the stress intensity factors of mixed mode crack problems, using the element free Galerkin method [4]. Accordingly, another new method was proposed to determine the SIFs for the indentation problem based on the conservation integral [5]. In contrast, Numerical experiments were conducted to evaluate the effectiveness of two proposed techniques on near crack tip singularity [6]. In fact, an observation was stated that a general theoretical model for automatically evaluating the increments of crack growth during a loading process does not exist yet [7]. Experiments work on a central crack specimen with holes, to predict the crack path and to evaluate the mixed mode stress intensity factors (K I & K II ) [8].
Somehow, this limit nowadays still characterizes the computer codes for crack growth simulat ion in elastic and elasto-plastic materials. Generally, the software should address three different aspects of the problem, specifically are; (i) SIFs evaluation, (ii) the crack propagation direction, and (iii) mesh modification to accommodate the crack advancement. The first and second aspects can be solved with no particular d ifficult ies and many commercial FE codes have built-in capabilit ies to evaluate the SIF. The last aspect is much mo re co mplex and is rarely solved by FE codes, e.g. FRANC2D, FRANC3D and ZENCRA CK [9]. Eventually, the FE model is modified in order to accommodate the obtained evaluation on short and straight crack propagation directions. The whole process is repeated many times and both is time-consuming and a source of errors if it is performed manually. This is part icularly t rue for the mesh modification. However, software on the other hand, is capable in handling the solution process in an automat ic way.
However, in the current study, we developed a new fin ite element programme for crack growth using the source code written in FORTRAN. In order to predict the crack growth trajectory and to evaluate the SIFs under mixed modes (I & II) loading in the frame of LEFM. The displacement extrapolation techniques with the adaptive refinement mesh method are used, to determine the stress intensity factors on two different geometries. Specifically are; four points bending plate and, centre double cracked plate with two-holes. Additionally, the obtained results of the current study are evaluated and validated with other relevant experimental and numerical results selected from the literature.

Generation of Singular Element
An unstructured triangular mesh is automatically generated by employing the advancing front method. The singular elements have to be constructed correctly to get a proper field of singularity around the crack tip as shown in Fig. 1. The number o f elements depends on the distributed nodes around the crack tip, wh ich can be set by the user as shown in Fig. 2.
In fact, the natural triangular quarter point elements are used [10] instead of the collapsed quadrilateral element, which is proposed by Barsoum[11]. The success of the adaptivity in general depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator. The importance of these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specificat ion of the boundary and mesh size d istribution over the domain. Moreover, several comprehensive numerical tools have been developed to enhance the accuracy at the crack tip. In this work, the unstructured triangle mesh is automatically generated by employing the advancing front method. Many researchers applied the FEM with remeshing to study the fracture propagation and its SIFs analysis ( [12]; and [13]), though it is not an easy task. In contrast, it is almost impossible to automatically remesh finite elements of an arbitrarily growing crack [14]. One of the advantages of the advancing front method is that the new triangle element formation is coinstantaneous with new node generation, and this advantage makes it possible to control the shape and size of the element through the adjustment of the location of the node [15]. However, a lot o f intersection checking between the new generated triangular elements and the existing elements must be computed in order to ensure that the triangular elements are valid. Accordingly, an observation was stated that the algorithm for the advancing front method has been shown to be robust in two-dimensional mesh generation of triangles and could be extended easily to generate quadrilaterals [16].

Mesh Generati on and Adapti ve Refinement Mesh
In general, the s maller mesh size g ives more accurate fin ite element appro ximate solution. Ho wever, reduction in the mesh size leads to greater computational effort. The adaptive mesh refinement is emp loyed as the optimization scheme. This scheme bases on a posteriori error estimator which is obtained from the solution fro m the previous mesh. Basically, the success of the adaptivity in overall is depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator.
The importance of these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specification of the boundary and mesh size distribution over the do main. An examp le of geo metry shows the whole process of generating the mesh is illustrated in Fig.3. The geometry of a p late with six holes and two notches is illustrated Fig. 3a. Six connector lines as shown in Fig. 3b, are fo rcing the internal boundaries to be the continuous part of the external boundary. Fig. 3c shows the cutting out of the rosette templates around each crack tip. The background mesh for this domain is then set up automatically using dichotomy technique as shown in Fig.3d and Fig. 3e shows the conventional mesh being generated by the advancing front method. The first generation produces mesh with initial size set by user. Later, during adaptive refinement, this first generated mesh will be taken as the background mesh. In Fig.  3f, for each rosette template, quarter-point elements are then constructed. Fig.3g, shows the enlargement of the quarter-point element at each crack t ip. In general, the smaller mesh size gives more accurate fin ite element approximate solution. However, reduction in the mesh size leads to greater computational effort.
The strategy used to refine the mesh during analysis process as follo ws [17]: (i) Determine the error norm for each element Where σ is the stress field obtained fro m the fin ite element calculation and σ* is the smoothed stress field inclusion of quarter point elements. Where h e is the old element size and p is the order of the interpolation shape function.

The Dis placement Extrapolation Techni que
One of the simplest and most frequently used methods is displacement ext rapolation technique [18]. It consists typically in the effective SIF concept by which, the fracture evaluation can be easily carried out [19]. The crack length was evaluated by linear-elastic analysis fro m the compliance of single-edge-notched specimen in three-points bending test [20]. A new fin ite element of quasi-static brittle fracture was developed, based on computational framework for three-dimensional cracks propagation bodies. They uses consistent thermodynamical framework for crack propagation in elastic solids [3].
The asymptotic expression for the displacement normal to crack p lane, v under mode I loading [21]: Where K I is the stress intensity factor for mode I, E is the modulus of elasticity, v is the Poisson's ratio, K an elastic parameter defined in equation (6) The near tip nodal displacements at nodes b, c, d and e shown below in Fig. 3, are of interest. The displacements are extrapolated by evaluating Equation (5) along the crack faces = ± . Particu larizing for nodes b and c on the singular element at the upper face of the crack yields: Where L, is the length of the element side which is connected to the crack tip. Equations (7) and (8) can be solved to obtain the value of K I as: The nodal displacements at the other two nodes can be evaluated by the similar way. If one considers the all four nodal displacements, the stress intensity factor expression becomes: Under pure mode II loading, the displacement tangential to the crack plane, u is given by: (1 ) Similarly, the stress intensity factor for mode II, using the nodal displacements of the four nodes can be estimated by: 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 d irection of crack trajectory such as the maximu m circu mferential stress theory, the maximu m energy release rate theory and the minimu m strain energy density theory. The maximu m circu mferential stress theory asserts that, for isotropic materials under mixed-mode loading, the crack will propagate in a direction normal to maximu m tangential tensile stress. In polar coordinates, the tangential stress is given by 2 1 3 cos cos sin 2 2 2 2 The direct ion normal to the maximu m tangential stress can be obtained by solving Which can be solved as:  In order to ensure that the opening stress associated with the crack direction of the crack extension is maximu m, the sign of should be opposite to the sign of [19]. The two possibilities are illustrated in Fig. 4. The criterion for crack to propagate fro m crack t ip is based on the material toughness, C K . If the calculated stress intensity factor, I C K K ≥ then the crack will propagate to the direction 0 θ expressed by Equation (11). The crack increment length a ∆ is taken 10%-20% of the initial crack length a , inversely proportional to the ratio of II I K K . The ratio represents the mixed mode p roportionality, therefo re shorter increment length should be taken to carefully justify the crack path curvature when is relat ively large compare to [13]. Thus the crack length increment is approximated by the Lagrange interpolat ion as:

Numerical Analysis and Validation
A developed finite element model of crack growth for brittle materials has been employed on four points bending geometry and, double centre cracked plate with two holes under mixed mode loading conditions.

Four Poi nts Bending Geometry
having the dimensions of B = 2 mm, W = 5 mm and length > 43 mm, made of Al 2 O 3 -Ceramics. An enlargement of mesh around crack t ip represented the elements around the crack tip is illustrated in Fig.3. Four points bending geometry with final adaptive mesh are shown in Fig. 5.
The stresses in this loading case were also computed. The geometric functions of F I and F II as follo ws [22]: In this study, a comparison and analysis were applied between the present results with those empirical results selected from the literature ( [3] and [22]). The stress intensity  Table 1 and  Table 2 with d/W = 0.5. The dimensionless form of the estimated stress intensity factor was obtained by using Eq. (17), presented for the range of 0.1 < a/ W < 0.4. In order to show the trends of the dimensionless SIFs with the ratio value of initial crack per width under Mode I conditions. Fig. 6 shows, the comparison between the current study results of the dimensionless SIFs and the emp irical results [22]. The relat ion showed that the values of SIFs and a/W increased from 0.31 with a/W = 0.1 to the maximu m dimensionless SIF value of 2.1 accordingly with a/W = 0.8, as shown in Fig. 6.   This section presents comparison of the current results with these results of Gürse and Miehe [3] in terms of crack propagation trajectories. Fig. 8 shows, the comparison results of predicted crack gro wth of the current study with those results of [3]. In this case, the crack direction was dominated by Mode I stress intensity factor, which became larger by the crack growth non-linear to the top right of the geometry. Due the shear stress which is generated by of Mode II loading condition. Further, the effectiveness of the shear stress was significant on crack trajectory with different direct ions, which is mostly reflected on the material properties as shown in Fig. 9. However, the representation of the deformation was controlled by the user, to show the deformat ion step by step with direct ions through the initial crack propagation until the final crack as illustrated in fig. 9. Obviously, the current results of this developed fin ite element programme has good agreement as shown in fig. 8b.

Center Double cracked Pl ate with Two-Holes under Mi xed Mode Loading
A comp licated mixed mode fracture problem was studied to demonstrate the performance o f the developed programme on double centre cracked plate with two holes. The geometry and its final adaptive mesh are shown in Fig. 10. Meanwh ile, the enlargement of the rosette around both crack tips is shown in Fig. 11.  The geometry was made of high carbon steel and had the following parameters: Young's modulus E = 2.1 ×10 5 MPa and Poisson's ratio ν = 0.3. The trend showed the relationship between the stress intensity factors of I K and II K as illustrated in Fig. 12.
In the beginning, the crack exh ibited a pure Mode I state, which increased the K I values, whereas in Mode II state, values became lower so that the crack curved closely to the left side of the hole and then curved downwards where and tended to decrease. Therefore, when the crack direction declined fro m both holes with the curvature, both and decreased. However, when the crack growth linearly tended to increased again whereas, decreased. Furthermore, based on this phenomenon, it could be understood that Mode I significant effect than Mode II, and this was due to the brittleness of this material. A comparison study between the current study and the experimental results [8] was performed, to validate the accuracy and reliability of this developed FE model. Notably, all the crack lengths were measured along the cracked path. The crack trajectory of the nu merical results for this study is shown in Fig. 13. It appeared that, in the v icin ity of the holes, the crack direct ion was curved as a consequence of the mixed mode (I and II) loading. The crack was propagated non-linearly towards the hole by the attraction of the holes, which are generated high stresses from both sides of the crack t ip.  [8] revealed the same behaviour of the crack growth trajectory as illustrated in Fig. 14. Furthermore, Fig. 15 shows the symmetric results of the maximu m principal stress with the crack propagation trajectory. Figure 14. Current work results of an enlargement of a crack trajectory for the CCT specimen with two holes Figure 15. The triangle points correspond to experimentally obtained crack path, where the full line holds for one parameter and the dash line for two parameter fracture mechanics [8] The predicted crack trajectory fro m the numerical results of the current results was obviously shows, a similarity and a good agreement to the experimental results [8]. In fact, the crack direction is influenced whenever there is hole in the geometry and this is due to the huge of stresses were generated around the holes.

Conclusions
In this paper, a co mprehensive Finite Element model was developed using the source code written in FORTRAN© with an advancing front method for crack growth analysis. Displacement Extrapolat ion Technique (DET) was emp loyed to evaluate the SIFs under mixed mode loading conditions. The automatic crack propagation was characterised by the successive propagation steps performed without user interaction. The developed FE programme allo wed the user to control the process by changing the size of the crack gro wth increment, maximu m and minimu m element size in the mesh and init ial mesh size for the background mesh. Overall, all these advantages were successfully revealed the reliability and the capability of this new developed FE model in dealing with plane crack behaviours. In fact, in term of flexib ility the user could also control the problem type, which is either a plane strain or a plane stress to clarify the deformation. The comparisons and analysis shown that, this developed program was truly evaluated and validated with relevant research works under different methods selected from the literature.