A Numerical Modelling of Mixed Mode Crack Initiation and Growth in Functionally Graded Materials

The main objective of this work is to present a numerical modeling of crack propagation path in isotropic functionally graded materials (FGMs) under mixed-mode loadings. The displacement extrapolation technique (DET) and the maximum circumferential stress (MCS) criterion are investigated in the context of crack growth in functionally graded beam subject to three and four bending conditions. Using the Ansys Parametric Design Language (APDL), the variation continues of the material properties are incorporated by specifying the material parameters at the centroid of each finite element (FE) and the crack direction angle is evaluated as a function of stress intensity factors (SIFs) at each increment of crack extension. In this paper, two applications are investigated using an initial crack perpendicular and parallel to material gradient, respectively. The developed approach is validated using available numerical and experimental results reported in the literature.


Introduction
Functionally graded material (FGM) is a material solution and a material concept used for a new class of advanced composites characterized by gradual variation in composition, microstructure and material properties. These graded materials have emerged from the need to enhance material performance. FGMs are ideal candidates for the applications requiring multifunctional performance. For example, the ceramic/metal FGMs can be designed to reduce thermal stresses and take advantage of the heat and corrosion resistances of ceramic and the mechanical strength of metal. Due to the variation of properties, the fracture and crack propagation behavior in FGMs is very complicated.
In contrast to the extensive investigations on the mode-I crack propagation in FGMs; few studies on the mixed-mode crack propagation of FGMs have been reported. Oral and al. 1 performed experimental and numerical investigations on crack initiation in FGMs under mixed-mode loadings. Zhang and Paulino 2 used cohesive zone models to simulate two-dimensional mixed mode dynamic crack propagation in FGMs. Fantuzzi and al. 3 investigated the dynamic behavior of moderately thick FGM plates with geometric discontinuities and arbitrarily curved boundaries, using the Generalized Differential Quadrature Finite Element Method (GDQFEM) as a numerical approach 4 . Benamara and al. 5 employed the strain energy density (SED) criterion and the displacement extrapolation technique, to study the crack growth in FGMs, under mixed mode loading. Dimitri and al. 6 proposed the application of the level set method combined with the numerical extended finite element method (XFEM) to predict the fracture direction of propagation within a specimen, and to compute the stress intensity factor for cracked plates under different loading conditions. Using the maximum energy release rate criterion and interaction integrals, Kim and Paulino 7 simulated the crack propagation in FGMs, under mixed-mode and non-proportional loading. Kirugulige and Tippur 8 studied the influence of material gradient on the crack propagation during a mixed-mode dynamic fracture test. Gu and Asaro 9 analysed the mixed-mode fracture of FGMs with consideration of a crack perpendicular to the material gradient. Doan and al. 10 investigated a numerical solution of dynamic crack propagation in FGMs with different configurations, using a hybrid phase field model, which is particularly suitable for dynamic crack propagation. Steigemann and al. 11 presented the numerical modeling of quasi-static crack growth in FGMs using the compact tension shear (CTS) specimen. Benamara and al. 12 examined the effect of inclusions and cavities on the crack trajectories in FGMs, using the minimum strain energy density (MSED) and maximum circumferential stress (MCS) criteria. Rousseau and Tippur 13 performed experimental and numerical investigations on crack growth in an epoxy/ glass FGM beam subjected to four-point bending. Tilbrook and al. 14 performed experimental and theoretical studies of fracture and fatigue behaviour of FGMs and discussed the effects of material gradation on stress intensity factors and crack extension direction. Cheng and al. 15 22 studied a numerical crack propagation using the scaled boundary finite element method (SBFEM) for isotropic and orthotropic FGMs, the crack propagation direction is determined from the maximum circumferential stress (MCS) criterion. In this approach, the computational domain is discretized with a combination of quadtree cells and polygon meshes. 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. Calculation by these methods can be carried out far from the crack tip, so as not to involve the singularities. The main disadvantage of these methods is that the stress intensity factor components, /and /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.
This paper focuses on two-dimensional mixed-mode crack propagation in FGMs using the finite element method (FEM). Using the APDL code 23 , the displacement extrapolation technique (DET) and the MCS approach are used to determine numerically the SIFs and the crack direction in functionally graded beam subjected to three and four bending conditions, respectively.
The present paper is organized as follows: Section 2 presents the numerical evaluation of SIFs and the crack growth direction. Section 3 clarifies the methodology of crack extension modeling in FGMs. To validate all the previous development, Section 4 presents many numerical results which, whenever possible, are compared to available experimental and numerical solutions. Finally, the major conclusions are summarized in Section 5.

Numerical evaluation of SIFs
To calculate the SIFs for material gradients, general load cases, and complex crack geometries, numerical methods are habitually employed: the modified crack-closure integral MCC (Rybicki and Kanninen 24 , Kim and Paulino 25 ), the interaction integral (Topal and Dag 26 ), the displacement correlation technique DCT (Sabuncuoglu and al. 27 ). In this paper, the displacement extrapolation technique DET 28,29,30 are used to calculate the SIFs K I and K II for homogeneous materials are modified to calculate the SIFs FGMs. The advantage of this method is that, it is well suited for multiple cracks and it can be performed faster. The SIFs Factors can be determined from: Where: E tip and ν tip are the Young's modulus and the Poisson's ratio at the crack-tip location, respectively.
K tip =3-4 ν tip for plane strain, K tip =3-ν tip /1+ ν tip for plane stress. L is the length of the element side connected to the crack-tip. u i and v i (i=b, c, d and e) are the nodal displacements at nodes b, c, d and e in the x and y directions, as shown in Fig. 1.
In this paper, the special quarter point finite elements proposed by Barsoum 31 are used to obtain a better approximation of the field near crack-tip. The mid-side node of the element in the crack-tip is moved to 1/4 of the length of the element.

Crack direction criteria
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 energy release rate (MEER) criterion based on Griffith's theory 32 and the minimum strain energy density (MSED) criterion introduced by Sih 33 .
In the present paper, the maximum circumferential stress (MCS) criterion proposed by Erdogan and Sih 34 is investigated to determine the crack direction angle at each increment of crack extension. They postulated that the crack extension angle is begins along the maximum value of the circumferential tensile stress, and can be expressed as follows (3) In order to ensure that the opening stress associated with the crack direction of the crack extension is maximum, the sign of θ 0 should be opposite to the sign of K II . The two possibilities are illustrated in Fig. 2.
The MCS theory is simple, easily implementable and gives good results. This approach requires several remeshing at each increment of crack propagation. It also offers the advantage of being able to be used for elastoplastic materials in the case of confined plasticity (Bocca et al. 35 ). The displacement extrapolation technique used in the calculation is explained. The advantage of this method is that, it is well suited for multiple cracks and it can be performed faster.

Methodology of Crack Modeling
In this paper, the APDL code has been employed for creating the program to simulate the mixed-mode crack propagation in functionally graded beams. The displacement correlation technique and the maximum circumferential stress theory are used, to determine the SIFs and the crack direction, respectively.
The automatic mixed mode crack growth is characterized by successive propagation steps performed without user interaction. Each step consists of: 1. Setting the geometrical with initial crack and input material properties data of the problem.

Discretization of the model by plane183 elements:
For the mesh generation of the cracked plate, the element type 'PLANE183' of ANSYS code is used, as shown in Fig. 3(a). The main characteristic of this element is the fact that the middle node of the element is translated on one quarter of the element side. This transformation changes the shape function to model the singularity at the crack tip ( Fig. 3(b)).
Better stress field modeling will give more accurate results about the SIFs. 3. Mesh generation; refining around the crack-tip: The finite 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. 4. FE analysis. 5. Calculation of stress intensity factors: Using finite element implementations the SIF coefficients may be calculated using displacement at points near  crack tip. At each increment of crack extension Δa, the expressions for SIFs K I and K II are calculated using Eqs.1 and 2. 6. Calculation of crack propagation direction: The crack direction angle θ is evaluated as a function of SIFs at each increment of crack extension Δa, using Eq. 3, 7. Crack arrests? If yes, go to step 8. If no, go to step 9. 8. Stop. 9. Automatic deletes of the mesh and the crack segment, as shown in Fig. 3(a). and 4(b). 10. Automatic Insertion of the new crack segment for new crack tip position (Fig. 4(c)). 11. Return to step 3: The algorithm is repeated until ultimate failure of material or by using another criterion for termination of the simulation process. Fig. 4 illustrates the mixed mode crack propagation mechanism proposed in this study. Fig. 5 shows the crack propagation procedure at each step, based on the combination of the FE analysis and the MCS approach.

Numerical Analysis and Validation
The main objective of this study is to present a numerical modeling of crack propagation path in functionally graded materials (FGMs) under mixed-mode loadings. In order to confirm the validity of present approach, two examples are presented: • Four-point bending specimen with crack perpendicular to material gradation. • Three-point bending specimen with crack parallel to material gradation.

Four-point bending specimen with crack perpendicular to material gradation
In this problem, consider a notched functionally graded beam under four-point bending. The beam contains solid A-glass spheres of mean diameter 42 µm dispersed within an epoxy matrix. A crack of length a=5.5 mm is assumed to exist perpendicular to the material gradation. Fig. 6(a) shows the specimen geometry, applied forces and boundary conditions. Three crack positions: A, B and C, located respectively at ξ = 0.17, 0.58 and 1 are considered in this study. In Fig. 6(a), u x and u y indicate the displacements according to (Ox) and (Oy), respectively.
It is noted that this geometry was used by Rousseau and Tippur 13 to carry out experimental tests of fracture and crack propagation in FGMs and then, has been investigated in several numerical works 17,21,22,37 .
The structure considered is meshed by quadratic elements with 8 nodes and particularly, special elements were used to characterize the singularity around the crack. This mesh will be used to calculate the SIFs K I and K II using DET in order to then determine the crack direction at each increment of crack extension. The determination of the SIFs K I and K II for three crack positions is performed under plane stress conditions. Fig. 6(a) shows Geometrical model of functionally graded beam. Fig. 6(b) illustrates the FE mesh of bending specimen with different crack positions and, Fig. 6(c) shows the mesh detail around the crack tip. Fig. 7 illustrates the variations of Young's modulus E and Poisson's ratio ν. The numerical values of material properties E, ν and fracture toughness K IC in the graded region are illustrated in Table 1.     In order to validate the finite element simulations obtained from present approach, we adopt the experimental data reported by Rousseau and Tippur 13 for crack initiation angles θ 0 and crack trajectories of initial crack located at ξ=0.17, ξ=0.58, and ξ=1.00, respectively see figures 8(a)~(c) Using present method, Fig.9(a)~(c) shows the numerical results obtained for crack trajectories and crack initiation angles θ 0 .Crack propagation simulation was carried out assuming plane stress conditions and a crack propagation length Δa =1 mm constant. The FE results are performed for three initial crack positions, using MCS criterion. There is reasonable agreement in crack initiation angles and crack trajectories between numerical and experimental data reported by Rousseau and Tippur 13 In Fig.10, we performed a FE simulations of crack propagation initiated at position A (ξ=0.17). The crack path is predicted, using the MCS criterion and the displacement extrapolation technique DET. The predicted crack path obtained from displacement extrapolation technique DET is compared with those obtained by displacement correlation technique DCT. Excellent agreement is observed between both approaches. Fig. 11 compares the predicted crack path of the developed method with reported by Ean and al. 21 , using To better show the influence of the crack position located in the functionally graded beam on the final prediction of the crack path, we have shown in Fig.13 the three trajectories obtained for three positions of initial crack: A, B and C, respectively.
The results proved that the crack trajectory was dependent on the initial crack location in graded region. The estimated propagation path for crack initiated in position A (ξ=0.17) has a maximum value of initial crack angle/, the more the initial crack position ξ moves towards the position (ξ=1), the more the initial crack angle decreases. The prediction of the final cracking path is related to the initial position of the crack in the FGM material. Table 2 compares the prediction initiation angle / at initial load step estimated from present model with those obtained in the references 13,17,21,37 , using experimental data and numerical results based on the finite element method (FEM), extended finite element method (XFEM) and scaled boundary finite element method (SBFEM). The results in Table 2 clearly indicate that the results are very close to each other.

Three-point bending specimen with crack parallel to material gradation
In order to validate the FE simulations by present approach, we adopt the experimental data of mixed-mode crack propagation in ZrO 2 /NiCr FGMs, reported by Xin and al. 38 . Fig. 14 shows a functionally graded beam subjected to three-point bending conditions. In this problem, the direction of the material gradient is parallel to the crack direction.
In this study two types of material are analyzed: FGM-A corresponds to the specimen with a crack situated parallel to material gradation and propagating along the direction of decreasing elastic modulus.
Non-FGM: with the material composition of 50% NiCr and 50% ZrO 2 is also tested under the same experimental conditions. Fig. 15 illustrates the variation of the elastic modulus along the width of the FGM-A specimen. The elastic modulus decreases monotonously from 201 to 119 GPa with increasing the NiCr content.
The experimental results reported by Xin and al. 38 for mixed mode crack propagation paths in FGM-A and non-FGM specimens are shown in Fig. 16(a)-(b), respectively. In Fig. 16(a), the notch-tip of the FGM-A specimen exists in the 100% ZrO 2 layer and the crack propagates through the graded region including the layers with 100%, 90%, 80%, 70%, 60% and 50% ZrO2. Fig. 15b corresponds to the crack propagation path in non-FGM specimen. scaled boundary polygon elements. Excellent agreement was observed between two configurations. Fig. 12 compares the predicted crack path obtained from present method for ξ= 0.17; 0.58 and 1, with obtained from experimental data reported by Rousseau and Tippur 13 . The comparison shows that the finite element simulations are very acceptable for predicting the final propagation path in functionally graded beam. In Fig. 12, the coordinates (0, 5) corresponds to the initial crack tip.    The three-point bending specimen is meshed with quadratic elements with 8 nodes and particularly, elements were used to characterize the singularity in the vicinity of the crack-tip. The FE model consists of approximately 2686 eight-node isoparametric elements. The crack propagation simulation was carried out assuming plane strain conditions with crack increment length Δa=1 mm. The value of Poisson's ratio is taken as 0.3 in the entire model. Fig. 17 shows a closeup view of final crack growth path obtained from present approach, using FGM-A. Fig. 18 shows the numerical results obtained for crack trajectories in FGM-A and non-FGM. These results are compared with the XFEM simulations and from experimental data reported in reference 38 . The comparison shows that the FE simulations obtained by developed method are very acceptable for predicting the crack propagation path in the case of functionally graded and non-graded materials and further confirms the validity of the present approach.

Conclusion
In this study, the quarter-point singular elements proposed by Barsoum are used to consider the singularity of stress and deformations fields at crack-tip in FGMs.
The displacement extrapolation technique was used to obtain the SIFs at crack-tip and to predict the final crack trajectory in isotropic FGMs by evaluation, for each propagation step, the kinked angle using maximum circumferential stress theory.
The methodology of mixed mode crack propagation modeling proposed in this study has been used successfully to predict the final path of initial crack located perpendicular and parallel to material gradation. The crack trajectories obtained from present study and those reported in the literature agree very well with each other and further confirms the validity of the present approach.
Based on the results, it was recommended to add further development the APDL code to simulate crack propagation in orthotropic FGMs, under thermal and thermo-mechanical loadings.