1. Introduction

Structural failure can be generally associated with one or more fracture of the materials making that structure. In fact, fracture mechanics deals with the irreversible process of rupture due to nucleation or sudden crack and crack growth. The influence of cracking on the structural strength has been widely appreciated since the end of the 19th century. However, some aspects of its nature and influence still remain unknown[1]. The use of crack propagation laws based on stress intensity factor range is the most successful engineering application of fracture mechanics. In the elastic fracture analysis, the stress intensity factors sufficiently define the stress field close to the crack tip and provide fundamental information of how the crack is going to propagate. Basically, the estimation methods can be categorized into two groups, those based on field extrapolation near the crack tip 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 disadvantage of these methods is that the stress intensity factor components, K_I and K_II in mixed mode problems are either impossible or very difficult to be separated. Nevertheless, the first groups which are based on near-tip field fitting procedures require finer meshes to produce a good numerical representation of crack-tip fields generated to facilitate the calculation.

2. 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 partition 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 examples were performed to evaluate the generated stress intensity factors directly from the scaled boundary finite element solution for the singular stress field[3]. More recent work was developed variational 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 simulation 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 difficulties and many commercial FE codes have built-in capabilities to evaluate the SIF. The last aspect is much more complex and is rarely solved by FE codes, e.g. FRANC2D, FRANC3D and ZENCRACK[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 particularly true for the mesh modification. However, software on the other hand, is capable in handling the solution process in an automatic way.

However, in the current study, we developed a new finite 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.

2.1. 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.

Figure 1. The cut and patch procedure of generating singular elements around a crack tip

The number of elements depends on the distributed nodes around the crack tip, which 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 specification of the boundary and mesh size distribution over the domain. Moreover, several comprehensive numerical tools have been developed to enhance the accuracy at the crack tip.

Figure 2. A typical arrangement of the natural triangle quarter-point elements around a 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 of 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].

2.2. Mesh Generation and Adaptive Refinement Mesh

In general, the smaller mesh size gives more accurate finite element approximate solution. However, reduction in the mesh size leads to greater computational effort. 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. 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 domain. An example of geometry shows the whole process of generating the mesh is illustrated in Fig.3.

Figure 3. The mesh generation stages

The geometry of a plate with six holes and two notches is illustrated Fig. 3a. Six connector lines as shown in Fig. 3b, are forcing 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 tip. In general, the smaller mesh size gives more accurate finite 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 follows[17]:

(i) Determine the error norm for each element

(1)

Where σ is the stress field obtained from the finite element calculation and σ* is the smoothed stress field inclusion of quarter point elements.

(ii) Determine the average error norm over the whole domain

(2)

Where m is the total number of elements in the whole do ain. (iii)Determine a variable, for each element as

(3)

Where η is a percentage that measures the permissible error for each element. If the size of the element is reduced and vice versa. (iv) The new element size is determined as

(4)

Where h_e is the old element size and p is the order of the interpolation shape function.

2.3. The Displacement Extrapolation Technique

One of the simplest and most frequently used methods is displacement extrapolation 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 from the compliance of single-edge-notched specimen in three-points bending test[20]. A new finite 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 plane, v under mode I loading[21]:

(5)

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), A_i are parameters depending on the geometry and load on the specimen, and are the polar coordinates defined at the crack tip.

(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 . Particularizing for nodes b and c on the singular element at the upper face of the crack yields:

(7)

(8)

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:

(9)

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:

(10)

Under pure mode II loading, the displacement tangential to the crack plane, u is given by:

(11)

Similarly, the stress intensity factor for mode II, using the nodal displacements of the four nodes can be estimated by:

(12)

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

(13)

The direction normal to the maximum tangential stress can be obtained by solving . The nontrivial solution is given by

(14)

Which can be solved as:

(15)

Figure 4. Sign of the propagation angle

In order to ensure that the opening stress associated with the crack direction of the crack extension is maximum, 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 from crack tip is based on the material toughness, . If the calculated stress intensity factor, then the crack will propagate to the direction expressed by Equation (11). The crack increment length is taken 10%-20% of the initial crack length , inversely proportional to the ratio of . The ratio represents the mixed mode proportionality, therefore shorter increment length should be taken to carefully justify the crack path curvature when is relatively large compare to [13]. Thus the crack length increment is approximated by the Lagrange interpolation as:

(16)

3. 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.

3.1. Four Points Bending Geometry

having the dimensions of B = 2 mm, W = 5 mm and length > 43 mm, made of Al₂O₃-Ceramics. An enlargement of mesh around crack tip 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 follows[22]:

(17)

Figure 5. (a) Four-point bending geometry and (b) the final adaptive mesh

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 factors evaluation (), are tabulated in 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.

Table 1. Dimensionless stress intensity factor for the central cracked plate for Mode I

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 empirical results[22].

Table 2. Dimensionless stress intensity factor for the central cracked plate for Mode II

The relation showed that the values of SIFs and a/W increased from 0.31 with a/W = 0.1 to the maximum dimensionless SIF value of 2.1 accordingly with a/W = 0.8, as shown in Fig. 6.

Figure 6. Dimensionless SIF vs. (a/W) for mode I

Therefore, under Mode II a condition, the values of the dimensionless SIFs was decreased gradually, where a/W was increased respectively the maximum value of 0.8 as shown in Fig. 7. From these results, it indicated that under Mode I loading condition, high stresses were generated when a/W = 0.8, which increased the dimensionless SIFs and also exceeded the maximum stiffness of this material. Whereas, under Mode II loading condition shows a low effect which influenced the values of the dimensionless SIFs. In fact, this was due to the brittleness of this material.

Figure 7. Dimensionless SIF vs. (a/W) for mode II

At the beginning, the stress was higher under Mode I and after that, a stress exceeded the criterion value it’s decreased gradually where produced a low effect under Mode II.

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 growth 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.

Figure 8. Comparison results of the crack trajectory; (a) current study (b) Gürse & Miehe[3]

Figure 9. Three steps of crack propagation direction of four-point bending geometry

Further, the effectiveness of the shear stress was significant on crack trajectory with different directions, 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 deformation step by step with directions through the initial crack propagation until the final crack as illustrated in fig. 9. Obviously, the current results of this developed finite element programme has good agreement as shown in fig. 8b.

3.2. Center Double cracked Plate with Two-Holes under Mixed Mode Loading

A complicated mixed mode fracture problem was studied to demonstrate the performance of the developed programme on double centre cracked plate with two holes. The geometry and its final adaptive mesh are shown in Fig. 10. Meanwhile, the enlargement of the rosette around both crack tips is shown in Fig. 11.

Figure 10. The geometry and its final adaptive mesh

Figure 11. The enlargement of the rosette around both crack tips

The geometry was made of high carbon steel and had the following parameters: Young’s modulus E = 2.1 ×10⁵ MPa and Poisson’s ratio ν = 0.3. The trend showed the relationship between the stress intensity factors of and as illustrated in Fig. 12.

In the beginning, the crack exhibited 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 from 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.

Figure 12. The relationship of KI & KII with crack extension

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 numerical results for this study is shown in Fig. 13. It appeared that, in the vicinity of the holes, the crack direction 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 tip.

Figure 13. Numerical result of the crack path for the CCT specimen with two holes

These generated stresses obviously influenced the crack propagation direction as well as the values of SIFs. The experiment results[8] revealed the same behaviour of the crack growth trajectory as illustrated in Fig. 14. Furthermore, Fig. 15 shows the symmetric results of the maximum 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 from 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.

4. Conclusions

In this paper, a comprehensive Finite Element model was developed using the source code written in FORTRAN© with an advancing front method for crack growth analysis. Displacement Extrapolation Technique (DET) was employed 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 allowed the user to control the process by changing the size of the crack growth increment, maximum and minimum element size in the mesh and initial 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 flexibility 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.