# Numerical analysis of quasi-static fracture in functionally graded materials

###### Abstract

This work investigates the existing capabilities and limitations in numerical modeling of fracture problems in functionally graded materials (FGMs) by means of the well-known finite element code ABAQUS. Quasi-static crack initiation and growth in planar FGMs is evaluated. Computational results of fracture parameters are compared to experimental results and good agreement is obtained. The importance of the numerical fit of the elastic properties in the FE model is analyzed in depth by means of a sensitivity study and a novel method is presented. Several key computational issues derived from the continuous change of the material properties are also addressed and the source code of a user subroutine USDFLD is provided in the Appendix for an effective implementation of the property variation. The crack propagation path is calculated through the extended finite element method (X-FEM) and subsequently compared to available experimental data. Suitability of local fracture criteria to simulate crack trajectories in FGMs is discussed and a new crack propagation criterion is suggested.

###### Keywords:

Functionally graded material (FGM) Finite element method (FEM) Fracture mechanics Crack propagation Extended finite element method (X-FEM)^{†}

^{†}journal:

∎

## 1 Introduction

Functionally Graded Materials (FGMs) are those whose composition and hence their properties vary gradually as a function of the position. Since their introduction by Kawasaki and Watanabe (1987) in high temperature metal/ceramic aerospace components, FGMs have found a wide range of commercial applications including cutting tools, biomedical devices, optical fibers and wear resistant coatings (Uemura, 2003). In many of these applications, FGMs provide an attractive way for the designer to tailor the microstructure to specific operating conditions, while minimizing the difficulties associated with discrete material interfaces. Very often, however, fracture resistance constitutes the primary design criterion, and this fact has led to the development of a special branch of fracture mechanics devoted to the failure of this class of materials.

Until now, the fracture of an FGM under quasi-static loading, which is one of the predominant modes of material failure, has been investigated extensively (Eischen, 1987; Jin and Noda, 1994; Erdogan, 1995). The primary conclusion of these investigations is that the classical inverse square root singular nature of the stress field is preserved in FGMs, but the stress intensity factor (SIF) is influenced by the non-homogeneity of the material. Therefore, in a linear-elastic cracked FGM, SIFs play a significant role since they characterize the crack-tip stress and strain fields. The non-singular -stress, which represents the stress parallel to crack faces, is another factor affecting the crack growth behavior (Becker Jr et al, 2001).

The finite element method (FEM) has been widely used for fracture analyses of FGMs. Eischen (1987) has evaluated mixed-mode SIFs by means of the path-independent integral. Bao and Wang (1995) have investigated periodic cracking in graded ceramic/metal coatings. Gu et al (1997) have evaluated SIFs using the standard -integral. Bao and Cai (1997) have studied delamination cracking in graded ceramic/metal substrate under mechanical and thermal loads. Anlas et al (2000) have calculated SIFs by means of the path-independent integral. Marur and Tippur (2000) have investigated a crack normal to the material gradient by means of both the FEM and experiments. Dolbow and Gosz (2002) have calculated the SIFs through the extended finite element method (X-FEM). Kim and Paulino (2002a) have evaluated mixed-mode SIFs by means of the path-independent integral, the modified crack closure (MCC) and the displacement correlation technique (DCT). The -stress has also been computed by means of the FEM. Becker Jr et al (2001) studied -stress and finite crack kinking in FGMs. Kim and Paulino (2003) used a unified approach of the interaction integral method for evaluating SIFs and -stress in FGMs. Chapa-Cabrera and Reimanis (2002b, a) have also used the FEM to investigate crack kinking in graded composites.

On the experimental side, the difficulty and cost of manufacturing large-size fracture specimens amenable to testing has led most investigators to develop model FGMs. Of particular interest in this research is a model FGM based on polyethylene 1% carbon monoxide co-polymer (ECO), manufactured by selective exposure to ultraviolet (UV) irradiation (Lambros et al, 1999). These specimens have material characteristics mimicking ceramic-metal FGMs, i.e., stiffer and more brittle at one end, becoming gradually less stiff and more ductile at the other. Cracks in these FGM ECO specimens have been analyzed by Abanto-Bueno and Lambros (2006) in the experimental work that has been chosen to validate the present numerical analysis. The difficulty in performing this type of experiments has led many analysts to adopt numerical schemes and solve FGM-related fracture problems. Although boundary integral formulations have been used in some cases (Zhang et al, 2003; Riveiro and Gallego, 2013), the FEM is by far the approach most commonly adopted.

This work evaluates the performance of numerical tools in the computational assessment of cracks in FGMs by means of the well-known ABAQUS finite element (FE) code. Computational results of fracture parameters (SIFs and -stress) are compared with available experimental results and good agreement is obtained. The importance of the numerical fit of the elastic properties in the finite element model is analyzed by means of a sensitivity study and a new method is presented and evaluated. FEM capabilities in various key issues from the numerical point of view, such as the implementation of the property variation at the element level or the effect of the material gradation in the computation of fracture parameters, are examined in depth and, in order to overcome the existing limitations in commercial FE packages, the source code of a user subroutine USDFLD is provided and several improvements are suggested.

The crack propagation path is simulated through the X-FEM and a good agreement with the experimental results of Abanto-Bueno and Lambros (2006) is obtained. This is of particular interest since work previously reported in the literature on this subject is limited. Suitability of local crack propagation criteria to simulate crack trajectory in FGMs is discussed and a novel crack propagation criterion is proposed.

## 2 Model formulation

### 2.1 Specimen geometry and parameters

The experimental results reported in this study are taken from those obtained by Abanto-Bueno and Lambros (2006). Details of the experimental procedure are described in that work and in others related (Lambros et al, 1999; Li et al, 2000), so only the ones relevant to the current research will be described here.

Abanto-Bueno and Lambros (2006) manufactured polymeric model FGMs based on selective ultraviolet (UV) irradiation of polyethylene cocarbon monoxide (ECO). ECO is a very ductile semicrystalline copolymer that undergoes accelerated mechanical degradation when exposed to UV light, so that by gradually irradiating a sheet of the material from one end to the other, a sample with continuous in-plane property gradation from stiff and brittle to more compliant and more ductile can be obtained. A very thin sheet of in-plane dimensions 300150 mm was irradiated for times varying from 5 h to 300 h. Once irradiated, the sheet was divided parallel to the irradiation direction, and two samples of 150150 mm were obtained. One of these was then cut perpendicularly to the irradiation direction into 15 strips of 10 mm width, which were used in uniaxial tension tests to measure the Young’s modulus , failure stress , and failure strain as a function of distance along the ECO sheet. The remaining 150150 mm sample from the original sheet was used to generate SENT fracture specimens. Therefore, although the material property variation was measured independently of the fracture experiments, both originate from the same manufacturing process.

Abanto-Bueno and Lambros (2006) monitored the near-tip field using the optical technique of digital image correlation (DIC). The DIC measured displacement field was then used to extract the fracture parameters by performing a least square minimization of the asymptotic expression of the displacement field in the vicinity of the crack tip; as in the case of an homogeneous material, but with material properties evaluated at the crack tip position (Eischen, 1987).

The testing protocol of Abanto-Bueno and Lambros (2006) included mixed mode fracture experiments on the base homogeneous material and various graded FGM samples. Mixed-mode fracture is inherent to FGMs since for a crack inclined to the property gradation direction, the stress state near the crack tip is mixed-mode irrespective of the far field loading. In order to validate and develop a complete numerical investigation of the fracture process of FGMs, Abanto-Bueno and Lambros experimental work is especially interesting because evaluates the three characteristic geometries of mixed-mode fracture in FGMs. Thereby, near-tip mixity can be attained either by asymmetric external loading, as in the homogeneous case, or by placing the notch at an angle to the direction of mechanical property variation, or by a combination of both. The effect of each of these cases was investigated using three specimens labeled here FGM I, II, and III. The geometry, dimensions and measured variation of local material properties of the three specimens are shown in Table 1 and Fig. 1.

H | W | h | a | ||
---|---|---|---|---|---|

(mm) | (mm) | (mm) | (mm) | (rad) | |

FGMI | 75 | 70 | 37.5 | 30 | |

FGMII | 90 | 70 | 32 | 26 | |

FGMIII | 90 | 70 | 32 | 25 |

### 2.2 Finite element model

The simulations were developed in the latest version of ABAQUS (2013). The specimens dimensions correspond to those reported by Abanto-Bueno and Lambros (2006). The sample thickness, common to all homogeneous and FGM cases, was 0.406 mm and thereby plane stress conditions are assumed. Mimicking the experimental procedure, loading is applied as a fixed vertical displacement along the upper edge of the specimen, the vertical displacement is constrained in the lower edge and, in order to remove rigid body motion, the horizontal displacement is also set to zero at the lower right hand corner. Given the fact that in both the homogeneous and graded cases the material used fails by crazing while showing very little shear yielding, linear-elastic behavior is assumed in this work as it was also assumed in the data analysis of Abanto-Bueno and Lambros (2006).

#### 2.2.1 Application of material gradient

The assignment of material properties must reflect the property distribution in the FGM specimen being simulated. However, almost all of the FE approaches mainly concentrate on homogeneous materials or piecewise homogeneous materials; specific FE formulations relating to nonhomogeneous materials with continuously varying properties are scarce. Nevertheless the inclusion of continual spatial variation of properties in the FE formulation does not entail a computational problem, as the stiffness matrix may be determined by averaging across each element. Material properties can vary between elements or between integration points. The former leads to a discontinuous step-type variation in properties. Assigning element properties individually, or dividing a structure into numerous areas and then assigning properties to areas (Bao and Wang, 1995) may be inappropriate in failure analysis or crack path predictions, where local stress values may be of critical importance.

Santare and Lambros (2000) developed a formulation for graded elements, which automatically interpolate material properties within the element. These can substantially improve the solution quality based on the same mesh density, especially for higher-order graded elements. Kim and Paulino (2002b) have also investigated elements with an internal property gradient and reached similar conclusions. Their work differs in that the former samples the material properties directly at the Gauss integration points of the element, while the latter adopts a generalized isoparametric formulation.

Rousseau and Tippur (2000) developed a technique to assign spatially varying properties at integration points by defining properties as a function of temperature and providing the model with an initial temperature distribution that matches the elastic modulus variation desired. The assignment of a zero thermal expansion coefficient then eliminates unwanted thermal strains. This technique enjoys great popularity since it can be used in most of the commercial FE packages. However, it is not suitable for thermomechanical analyses and does not allow for differences between the gradient profiles of the Young’s modulus and the Poisson’s ratio. Furthermore, contrary to what is often assumed, Rousseau and Tippur’s technique is not able to define a non-linear continuous variation of the elastic properties in most of the FE codes since, in order to obtain a consistent variation between mechanical and thermal strains, nodal temperature values are interpolated within the element through shape functions one order lower than those used for the displacements. In the case of ABAQUS, an average value of the temperature in the nodes is passed to the integration points when using linear elements and an approximate linear variation is assumed in quadratic elements. The former produces a step-type variation in the elastic properties (i.e., homogeneous elements) and the latter translates into a piecewise linear variation, regardless of the order of the function describing the elastic gradient. Therefore, one should be careful when assigning non-linear material property variations by means of Rousseau and Tippur’s technique since, for coarse meshes and some gradient profiles, it can bring inaccuracies in the results. In ABAQUS this can be overcome by defining the gradation of properties through a user subroutine UMAT or USDFLD, since both are called at integration points. However, if a UMAT subroutine is used, the mechanical constitutive behavior of the material must also be programmed and hence, it is not possible to use the material models already implemented in ABAQUS. Consequently, the material gradient is implemented in this work through a USDFLD user subroutine. Material elastic properties are defined as a function of a field variable and its variation throughout the specimen is programmed in the subroutine. In addition, when computing the SIFs, the elastic properties in the crack tip must be defined and therefore a UFIELD subroutine is also embedded in the FORTRAN code in order to take into consideration as well the elastic properties variation at the nodes.

The source code of the subroutine is provided in the Appendix in order to allow other engineers to implement an effective continuous variation of the material elastic properties without requiring programming efforts. Another option could be to use the research codes FGM-FRANC2D (Kim, 2003) or WARP3D (Healy et al, 2012) since both include the gradation effect at the element level, based on the nodal-values approach (Kim and Paulino, 2002a). Both are freely distributed, open-source finite element codes with extended capabilities for fracture in FGMs, though the former is not yet available to the public and the latter does not have the capability to model plane stress conditions (Walters et al, 2006).

#### 2.2.2 Numerical fit of the material elastic properties

The variation of composition in FGMs depends on the production technique (Lambros et al, 1999; Butcher et al, 1999; Parameswaran and Shukla, 2000) and generally, the property variation tends not to mirror that of composition. If the spatial composition profile is known, property variation may be predicted by means of theoretical mixing laws. Their use is frequent in composites (Hashin, 1983) and has also been extended to FGMs (Reiter et al, 1997; Gasik, 1998). In these cases the assignment of the material property variation in the model is done straightforwardly, fitting the variation of the elastic properties through a function with the shape of the theoretical prediction, following the procedure mentioned in the previous paragraph. However, in such models predicted property variation is largely based in the assumed composite structure and, therefore, is usually limited in applicability and accuracy due to the geometric and micromechanical assumptions upon which the theoretical mixing laws are based.

Hence, material property variation is usually determined directly from experiment, being characterized by a sequence of experimental data as a function of the position, regardless of the form in which these data were obtained. Either by producing and testing individual homogeneous specimens with a range of compositions (Carrillo-Heian et al, 2001; Jedamzik et al, 2000), or by testing the graded material directly by means of indentation or ultrasonic techniques (Krumova et al, 2001) or by cutting and testing small, effectively homogeneous, specimens from a larger graded sample (Lambros et al, 1999; Rousseau and Tippur, 2000; Butcher et al, 1999), as in the case of the experimental work (Abanto-Bueno and Lambros, 2006) that serves as basis for the validation of the numerical model presented in this paper. To the authors’ knowledge, the numerical fit of this experimental data, its implementation in the numerical model and its effect on the computational calculations have not received the attention of the scientific community. Mostly, the numerical fit is based on a general approximation of all the experimental data, either by assuming a linear change in the material properties (Rousseau and Tippur, 2000) or by means of a polynomial function through a least squares fit (Oral et al, 2008a, b). Following the criterion of the authors of the experimental study (Oral et al, 2008a), a fourth-order polynomial function was chosen to approximate the data (Fig. 2). But, as it can be seen in Fig. 2, it is impossible to completely remove the differences between the measured elastic properties and the polynomial curve fitting. And, even though no systematic study of the problem has been published yet, it is reasonable to expect that, in fracture analyses of FGMs, an accurate fit of the elastic properties near the crack tip could be overriding, due to the dependence of the fracture parameters’ magnitude on the crack direction, property profile, crack-tip position and specimen geometry.

In order to rate the effect of the numerical fit of the experimental data in the calculation of the fracture parameters a complete sensitivity study is developed. First, the effect of an accurate fit of the local values of the elastic properties at the crack tip is analyzed. In order to do that, the value of the experimental point that characterizes the local elastic properties in the crack tip is modified a small amount (5% - 10%). For each of the FGM specimens evaluated, the new data points derived from these changes are shown in Table 2 and the corresponding polynomial curve fits can be seen in Fig. 3. SIFs and -stress for each curve fit are calculated and compared.

FGM I | FGM II | FGM III | |
---|---|---|---|

-10% | 242.72856 | 334.18539 | 347.14386 |

-5% | 256.21348 | 352.751245 | 366.42963 |

Expt. data | 269.6984 | 371.3171 | 385.7154 |

+5% | 283.18332 | 389.882955 | 405.00117 |

+10% | 296.66824 | 408.4481 | 424.28694 |

Next, with the aim of evaluating the influence of an accurate fit of the material gradient profile in the computation of the SIFs and the -stress, several polynomial curve fits of different orders are considered, containing all of them the experimental data characterizing the local property values in the crack tip, as shown in Fig. 4. Fracture parameters for each curve fit are calculated and compared. Based on the conclusions of the sensitivity analysis a new method is developed and evaluated.

#### 2.2.3 Calculation of fracture parameters

By far the most common concern pertaining to linear elastic fracture mechanics analysis is the accurate prediction of SIFs in arbitrarily shaped cracked bodies. There are usually several ways to calculate fracture parameters once the stress and displacement fields have been obtained. In the displacement-based stress intensity factor computation techniques, the SIFs are obtained by extrapolating from the displacement ahead of the crack tip in the asymptotic expression. These methods have the advantage that almost no additional calculation is necessary, but they require a high degree of mesh refinement and often suffer from instability as the crack tip is approached (Anderson, 2005). Also, for the FGM case, choosing the appropriate correlation points can be a difficult task (Tilbrook et al, 2005b). A more often-used procedure, the domain integral method, which is an energy approach based on the -integral (Rice, 1968) that has been proved to be efficient for homogeneous materials, is used in this work.

From a numerical and computational perspective, one of the challenges concerns the need for examining the limiting case of a vanishing contour for the proper evaluation of the -integral for crack tips in FGMs. This need stems from the fact that for some non-homogeneous materials and crack tip orientations, the integrand in the -integral is not divergence free (Chen et al, 2000). As a result, an evaluation of the integral on open contours will exhibit path dependence. Thus, the standard -integral along an integral path is defined as:

(1) |

where is the strain energy density and is the outward normal to the path . For a closed boundary as shown in Fig. 5, the -integral is formulated such that:

(2) |

Applying Gauss’s divergence theorem gives:

(3) |

where

(4) |

Equations of equilibrium in the absence of body forces take the form , and therefore:

(5) |

Substituting Eqs. (4) and (5) in Eq. (3) renders:

(6) |

Along crack sides, and the traction is also zero. Consequently and thus:

(7) |

In an homogeneous material, since , and , the -integral is path independent. For FGMs, generally , therefore and , whereby the -integral is related to the integral path. When the material properties only vary along the axis, and in consequence, for this case, the -integral is still path independent.

Considering a smooth function which has he value of unity on and zero on , the -integral given in (1) can be written in terms of a closed boundary integral:

(8) |

Applying the divergence theorem and assuming a Young’s modulus variation along the axis, the corresponding expression for the domain integral is obtained:

(9) |

Being the derivative of under the second integral with respect to the coordinate in . So, comparing with the homogeneous case, the second integral is an additional term which represents the effect of non-homogeneity. But if the domain integral is evaluated in a region sufficiently small around the crack tip, the value of the second term in Eq. (9) involving the derivative of is very small, essentially negligible (Gu et al, 1997). The same conclusions can be extended to the Interaction Integral (Shih and Asaro, 1988) as it is calculated using the same procedure. Under mixed-mode conditions many finite element packages, including ABAQUS, provide an interaction integral method to compute the SIFs since it is not straightforward to do it from a known -integral. The interaction integral is also used for the computation of the non-singular term of the -stress.

Therefore, a sufficiently refined mesh near the crack tip was used in all the specimens analyzed in order to allow for comparison between all cases. Although the standard domain integral combined with a fine mesh allows to reach accurate results, this greatly increases the numerical costs. In order to exploit the advantage of the contour integral, which does not require a very fine mesh, one should use the above mentioned research codes, FGM-FRANC2D or WARP3D, which take into account the additional term of the Domain Integral necessary to maintain the path independence. Therefore, the computation of the derivative of the strain energy density with respect to the crack front (normal) direction is a pending task for future developments of commercial FE software.

In this work, the entire specimens were modeled using eight-noded quadrilateral plane-stress elements with reduced integration (CPS8R). An example of the mesh used in the simulations is shown in Fig. 6. The near-tip discretization consisted of a focused mesh centered at the crack tip and the inverse square root singularity of the strain field at the crack tip is obtained by collapsing quadrilateral elements into triangular elements and placing the mid-side nodes at quarter-points from the crack tip. Approximately 30000 elements were used in total.

#### 2.2.4 Simulation of crack propagation

Crack propagation was simulated by means of the X-FEM. The X-FEM is a numerical technique for modeling discontinuities by local enrichment functions in the area of interest. This method allows to follow crack paths independently of the finite element mesh; being this feature especially important for FGMs, since the gradation of the mechanical properties may lead to complex propagation paths also in simple symmetric tests (Rousseau and Tippur, 2000). Cracks are modeled by means of the X-FEM using the cohesive segments method (Remmers et al, 2008) implemented in ABAQUS software. This method is based on the principle of phantom nodes (Song et al, 2006) and a traction-separation cohesive behavior.

In quasi-brittle materials a crack is usually assumed to grow when the tensile strength is attained at the current crack tip; propagation then occurs perpendicularly to the direction of the maximum in-plane principal stress. Therefore, crack initiation and propagation direction was predicted using the local criterion of the maximum principal stress (MPS). Based on the concept of local homogenization near the crack tip (Gu and Asaro, 1997), fracture criteria originally developed for homogeneous materials can be extended to non-homogeneous materials such as FGMs. The values of the initial failure stress and the subsequent fracture energy characterizing damage evolution are assigned taking into consideration the local failure properties recorded by Abanto-Bueno and Lambros (2006) (see Fig.1). Fracture energy variation was provided following the same procedure as with the elastic gradient implementation, although, based on the local propagation criteria established, the changes of fracture resistance as a function of the position only influence the amount of applied load necessary to propagate the crack. Notice that there is no effect of the load magnitude on the crack trajectory within the framework of linear elastic analysis. All the specimens were discretized with the same mesh density by means of approximately 450 four-noded quadrilateral plane-stress elements with reduced integration (CPS4R). An example of the mesh used in the simulations is shown in Fig. 7.

## 3 Results and discussion

### 3.1 Crack initiation

#### 3.1.1 Homogeneous Materials

The base homogeneous material studied is an ECO sheet irradiated uniformly for 50 h under UV light. Its elastic properties were measured by Abanto-Bueno and Lambros (2006) as Young’s modulus E=280 MPa and Poisson’s ratio . Geometry and dimensions are shown in Fig. 8.

By means of the procedure outlined in the previous section, values of , and corresponding to the experimentally recorded instant of crack initiation were calculated. Table 3 shows a comparison of the results obtained with the corresponding experimental values extracted by Abanto-Bueno and Lambros (2006) and the results obtained numerically by the same authors in collaboration with other researchers in a related work (Oral et al, 2008b). Good agreement is observed, especially for the non-singular term of the -stress, which is usually more difficult to obtain.

Results | ||||
---|---|---|---|---|

Expt. | Abanto-Bueno and Lambros (2006) | 0.903 | 0.245 | -0.784 |

Num. | Oral et al (2008b) | 0.793 | 0.212 | -0.992 |

Present | 0.812 | 0.291 | -0.647 |

Results | Case | ||||
---|---|---|---|---|---|

Expt. | Abanto-Bueno and Lambros (2006) | 0.554 | 0.039 | -4.272 | |

Num. | Oral et al (2008b) | I | 0.551 | -0.022 | -2.149 |

Present | 0.599 | -0.018 | -1.526 | ||

Expt. | Abanto-Bueno and Lambros (2006) | 0.755 | 0.179 | -0.069 | |

Num. | Oral et al (2008b) | II | 0.722 | 0.204 | -0.673 |

Present | 0.734 | 0.257 | -0.539 | ||

Expt. | Abanto-Bueno and Lambros (2006) | 0.969 | 0.224 | -0.930 | |

Num. | Oral et al (2008b) | III | 0.878 | 0.230 | -0.870 |

Present | 0.908 | 0.304 | -0.763 |

#### 3.1.2 Functionally Graded Materials

Young’s modulus distribution for the three graded specimens is shown in Fig. 2. In all cases, a constant value of 0.45 is considered for the Poisson’s ratio, in accordance with the assumption of the experimental work Abanto-Bueno and Lambros (2006). As described in the previous section, material property variation is implemented via a user subroutine USDFLD although in the present isothermal analysis, where the Poisson’s ratio is constant and a very fine mesh is used, results agree with those obtained using Rousseau and Tippur’s technique and therefore verify the subroutine implementation. Computed values of fracture parameters are compared with available experimental and numerical results in Table 4.

As in the homogeneous case, results agree reasonably well with the available experimental data, although there are some differences that need to be analyzed. The greatest discrepancies arise in the -stress that, being a second-order term, is more difficult to extract, both numerically and experimentally. Besides, analyzed geometries bring great changes in its value, which is also affected by the direction of the material property variation. All the results, and especially the -stress, are affected by the method chosen to extract them. As mentioned before, Abanto-Bueno and Lambros, in both their experimental (Abanto-Bueno and Lambros, 2006) and their related numerical work (Oral et al, 2008b) fit the experimentally measured or numerically computed displacements in the asymptotic displacement equation to extract the value of the fracture parameters, while in this work these are computed by means of the Domain Integral without the need of any postprocessing. Also, due to the limitations of the DIC method, experimentally measured displacements exclude a rectangular region around the crack tip. Whereas in the numerical model data as close to the crack tip as possible are considered, in accordance with the local homogenization criterion.

In the first graded specimen, a negative value of is obtained, in contraposition with the positive value extracted from the experiments. The value is relatively small and whence the difference is minimal, but the sign affects the direction of crack kinking. The crack growth analysis conducted by Abanto-Bueno and Lambros (2006) shows a positive kink direction and therefore, the negative sign of obtained in the numerical simulation correctly predicts the positive sign of subsequent crack kinking.

Given the unavoidable experimental error and taking into consideration that differences found were duly justified, results show that finite element calculations allow to obtain fracture parameters for the precise instant of crack initiation with reasonable accuracy.

### 3.2 Sensitivity analysis

#### 3.2.1 Effect of an accurate fit of the local material properties in the crack-tip

Among other factors, fracture parameters in non-homogeneous materials depend on the local properties at the crack tip. In order to quantify the influence of an accurate fit of the measured elastic properties in the vicinity of the crack, the value of the experimental point that characterizes the local value of the Young’s modulus in the crack tip is slightly modified (see Table 2 and Fig. 3).

(%) | (%) | (%) | |
---|---|---|---|

-10% | 0.0145% | 0.0219% | 0.0131% |

-5% | 0.0037% | 0.0055% | 0.0033% |

+5% | 0.0037% | 0.0055% | 0.0033% |

+10% | 0.0149% | 0.0219% | 0.0131% |

Test points | (%) | (%) | (%) |
---|---|---|---|

-10% | 0.591 (1.34%) | -0.017 (5.55%) | -1.509 (1.11%) |

-5% | 0.595 (0.67%) | -0.017 (5.55%) | -1.517 (0.59%) |

Expt. data | 0.599 | -0.018 | -1.526 |

+5% | 0.603 (0.67%) | -0.017 (5.55%) | -1.535 (0.59%) |

+10% | 0.608 (1.50%) | -0.017 (5.55%) | -1.543 (1.11%) |

Expt. Results | 0.554 | 0.039 | -4.272 |

As expected and it can be seen in Fig. 3, a small change in the value of one of the nearly 20 experimental points fitted does not bring many differences between the different least square fits derived from the changes proposed. In order to quantify these deviations, the average difference along the specimen width between the numerical fits emerged from the proposed changes and the numerical fit of the original experimental data is calculated in a percentage scale by means of the norm. Results obtained from this relationship are shown in Table 5.

Fracture parameters computed for each curve derived from these changes are presented in Tables 6 to 8. Percentage differences between the calculated values and the values of the fracture parameters corresponding to the numerical fit of the original experimental data are also shown for comparison purposes.

Test points | (%) | (%) | (%) |
---|---|---|---|

-10% | 0.720 (1.91%) | 0.252 (1.95%) | -0.527 (2.23%) |

-5% | 0.727 (0.95%) | 0.255 (0.78%) | -0.533 (1.11%) |

Expt. data | 0.734 | 0.257 | -0.539 |

+5% | 0.741 (0.68%) | 0.260 (1.17%) | -0.544 (0.93%) |

+10% | 0.747 (1.77%) | 0.263 (2.33%) | -0.550 (2.04%) |

Expt. Results | 0.755 | 0.179 | -0.069 |

Test points | (%) | (%) | (%) |
---|---|---|---|

-10% | 0.895 (1.43%) | 0.298 (1.97%) | -0.759 (0.52%) |

-5% | 0.902 (0.66%) | 0.301 (0.99%) | -0.761 (0.26%) |

Expt. data | 0.908 | 0.304 | -0.763 |

+5% | 0.914 (0.66%) | 0.306 (0.66%) | -0.764 (0.13%) |

+10% | 0.920 (1.32%) | 0.308 (1.32%) | -0.766 (0.39%) |

Expt. Results | 0.969 | 0.224 | -0.930 |

Even though the differences between the numerical fits that emerged from the proposed changes are extremely small in the three graded specimens, very significant discrepancies, between 30 and 1000 times the differences between the numerical curve fits, can be found among the values of the fracture parameters. As a result more attention should be paid to an accurate fit of the experimental elastic properties in the vicinity of the crack.

#### 3.2.2 Influence of the shape of the numerical curve fit

Besides the local properties in the crack tip, the shape of the material gradient could also influence the results. In order to evaluate the effect of an accurate fit of the elastic gradient profile, several polynomial curve fits of different orders are considered, containing all of them the experimental data characterizing the local property values in the crack tip, as seen in Fig. 4. Computed values for each curve fit of the SIFs and the -stress for each graded specimen are presented in Tables 9 to 11. Experimental results of Abanto-Bueno and Lambros (2006) are also shown for comparison.

Curve fit | |||
---|---|---|---|

Pol. order 3 | 0.584 | -0.016 | -1.499 |

Pol. order 4 | 0.599 | -0.018 | -1.526 |

Pol. order 9 | 0.601 | -0.018 | -1.523 |

Expt. Results | 0.554 | 0.039 | -4.272 |

Curve fit | |||
---|---|---|---|

Pol. order 3 | 0.745 | 0.261 | -0.539 |

Pol. order 5 | 0.742 | 0.259 | -0.525 |

Pol. order 9 | 0.743 | 0.260 | -0.533 |

Expt. Results | 0.755 | 0.179 | -0.069 |

Curve fit | |||
---|---|---|---|

Pol. order 3 | 0.975 | 0.329 | -0.768 |

Pol. order 5 | 0.972 | 0.327 | -0.774 |

Pol. order 9 | 0.965 | 0.322 | -0.772 |

Expt. Results | 0.969 | 0.224 | -0.930 |

Tables 9 to 11 show sensitivity in all the fracture parameters to the different curve fit considered in the three specimens evaluated. Although changes in the order of the polynomial function of the least squares fit do not affect the SIFs and the -stress as much as the variation of the local elastic properties at the crack tip, results prove the influence of the variation of the elastic properties by itself in the calculation of fracture parameters. Therefore, significant differences between the curve fit and the experimentally measured elastic properties can bring inaccuracies in the calculations. Consequently, the numerical fit must be as accurate as possible in the vicinity of the crack, but without neglecting the relevance of a rigorous fit of the rest of experimental data. In order to solve this problem a new method is proposed: a point to point lineal fit.

#### 3.2.3 Point to point lineal fit

With the aim of removing the differences between the data fit and the experimental data arises the possibility to implement directly in the numerical model the values of the Young’s Modulus extracted from the experiments. In order to do so, elastic properties are designated as a function of a field variable (or temperature) and each one of the experimental registered values of the Young’s Modulus in the specimen is defined with an assigned value of the field variable that depends on its position in the sequence of experimental data. Afterwards, the specimen is provided with a lineal distribution of the field variable that allocates the points depending on its position in the specimen. The finite element software assigns by default the value of the elastic properties in places where they are not defined considering a linear variation between the nearest two points registered. This results in a linear change between the data points as is indicated by the red line in Fig. 9 for the three specimens evaluated. The polynomial curve fit used until now is also presented with the aim of establishing a comparison.

In order to be accurate, this method is limited by the fact that the experimental data needs to be separated by a distance as similar as possible. However, this is usually not an obstacle since when property variation is determined directly from experiment, as in most cases, measurements are performed on points of the material containing a constant separation from each other with the aim of correctly record the variation of the material properties. Through the same procedure that has been used so far, fracture parameters are calculated and shown in Table 12. Experimental data from Abanto-Bueno and Lambros (2006) and the results from the common used polynomial fit are also shown.

Procedure | Case | |||
---|---|---|---|---|

Experimental | 0.554 | 0.039 | -4.272 | |

Polynomial fit | I | 0.599 | -0.018 | -1.526 |

PtP lineal fit | 0.602 | -0.018 | -1.529 | |

Experimental | 0.755 | 0.179 | -0.069 | |

Polynomial fit | II | 0.734 | 0.257 | -0.539 |

PtP lineal fit | 0.718 | 0.251 | -0.517 | |

Experimental | 0.969 | 0.224 | -0.930 | |

Polynomial fit | III | 0.908 | 0.304 | -0.763 |

PtP lineal fit | 0.935 | 0.296 | -0.942 |

As seen in Table 12, even though the experimental error is unknown, results improve when using the point to point lineal fit. Especially in the third specimen, where the polynomial function loses accuracy when fitting the experimental data in the vicinity of the crack (Fig. 9).

### 3.3 Crack growth

Crack propagation trajectories of the three specimens were calculated based on the X-FEM and the local fracture criterion of the MPS. Results are compared with available experimental data from Abanto-Bueno and Lambros (2006) in Fig. 10 and the performance of local fracture criterion in the prediction of the crack propagation path is analyzed.

As seen in Fig. 10 numerical prediction approaches reasonably well the experimental data. In the case of the first graded specimen, in both the numerical (Fig. 10 (b)) and the experimental results (Fig. 10 (a)) it can be appreciated that the crack propagation path deviates from the natural trajectory of the first mode of fracture. This is due to the fact that the crack is initially orientated perpendicular to the material gradient, creating an asymmetric stress field around the crack tip that induces mixed-mode loading causing crack deflection. As expected, the crack propagated towards the more compliant part. This is understood to occur as it results in greater release of elastic potential energy.

Crack propagation path for FGMII is shown in Fig. 10 (c) and (d). Under the same loading conditions in an homogeneous specimen, the crack propagation path would register a small initial deviation, due to the asymmetry of the loading within the initial crack orientation, to grow then according to the trajectory distinctive of the first mode of fracture. So that the initial asymmetry in the exact moment of the onset of crack propagation due to the position of the crack within the direction of mechanical property variation and the external loading translates over the trajectory exclusively in a small initial deflection, that comes across very well in the numerical simulation and which is also referred in Abanto-Bueno and Lambros experimental work (Abanto-Bueno and Lambros, 2006). After this initial deviation, the crack grows following the propagation path distinctive to the first mode of fracture since is oriented parallel to the material property variation

In FGMIII at the precise instant of crack initiation there is a mixed-mode fracture induced by the asymmetry of the loading within the crack position. However, in an homogeneous specimen with the same loading conditions, after an initial deviation the crack should grow according to the idiosyncratic trajectory of the first mode of fracture. As seen in Fig. 10 (e) and (f), this is not the case, as once the crack trajectory tries to get horizontal, forms an angle with the gradient direction and therefore, the asymmetric stress field around the crack tip induces mixed-mode loading causing crack deflection.

Although numerical predictions show reasonably good agreement with the experimental data, a close observation of the results of the first and third cases shows that the deviation of the trajectory from the crack propagation path corresponding to the first mode of fracture is slightly higher in the experimental case. This may be due to many factors and it is not surprising that there are some discrepancies between predicted and observed deflection trajectories, given the stochastic nature of crack initiation from a notch and the possible influences of the microstructure. In fact, assumed local homogenization criteria is expected to work well only for ideally brittle non-homogeneous materials where the toughness is taken to be independent of direction at a fixed point (Gu and Asaro, 1997). However, in a FGM specimen, due to its change in material composition, this is often not the case, and therefore the material toughness variation as a function of the position cannot be neglected. So that, in order to develop a more accurate crack propagation criterion for non-homogeneous materials, this should be stated as:

(10) |

Where is the energy release rate and is the critical value of the energy release rate according to which the crack starts to propagate, the toughness of the FGM. Although in the specimens evaluated in this work the spatial variation of the fracture resistance is unknown, based on the variation of the failure properties registered by Abanto-Bueno and Lambros (2006) and shown in Fig. 1, one could expect that it would increase the deviation predicted based on the effect of the elastic gradient of the material, as it varies in the same direction in the all the specimens evaluated.

Nevertheless, very little work has been published on crack propagation in graded regions. Crack trajectory becomes more complex after deviating from a straight path and therefore, appropriate models for the continued propagation of cracks forming an angle with the material gradient will require a fracture criterion which incorporates the effects of variations in local mechanical properties, crack-tip toughness, bridging and residual stresses, as well as changing crack shape (Tilbrook et al, 2005a). To the authors’ knowledge, this fracture criterion has not been developed yet.

## 4 Concluding remarks

In this work, performance of numerical tools to study the structural integrity of advanced materials is evaluated. Crack initiation and growth in planar FGMs is investigated by means of the well-known finite element package ABAQUS.

With the aim of overcoming the limitations of Rousseau and Tippur’s technique (Rousseau and Tippur, 2000), material gradient was implemented by means of a user subroutine USDFLD and its template is provided in the Appendix in order to facilitate the work of other other researchers and practitioners. When computing fracture parameters in a precise instant of the fracture process, commercially available finite element software has proven to obtain accurate results, although the use of the before mentioned research codes, FGM-FRANC2D or WARP3D, is recommended when the crack is not perpendicular to the direction of the material property variation in order to reduce CPU time.

However, when predicting crack propagation paths, available finite element codes can only get accurate results by means of local crack propagation criteria if the elastic gradient is the dominant influence on crack propagation. To develop a general purpose formulation and implementation for predicting crack trajectories in FGMs a new fracture criterion that incorporates the influence of variations in crack-tip toughness as well as other mechanical effects, is needed.

In addition, this paper emphasizes that more attention needs to be paid to the numerical fit of the experimental data characterizing the change in the material properties, and a new method is proposed to improve the accuracy of the results. The commonly used criterion based on a general approximation of the overall experimental data could introduce inaccuracies in the calculations within a fracture analysis, as small changes in the local elastic properties at the crack-tip have been shown to have significant influence in the values of the fracture parameters. Furthermore, an accurate fit of the shape of the material gradient is also relevant, as even by itself it influences the fracture parameters. A point to point lineal fit has proven to reach more accurate results.

###### Acknowledgements.

The authors gratefully acknowledge the financial support from the Ministry of Science and Innovation of Spain through the grant DPI2010.21590.CO2.01.## Appendix A. User subroutine USDFLD FGMII

C User subroutine for the implementation of

C a continuous variation of the material

C elastic properties between integration

C points.

C Code layout is provided the simplest way

C possible and comments are included

C conveniently with the aim that it can serve

C as a template. Fortran statements are

C placed at the beginning of each row based

C on text formatting, make sure they are

C placed in the appropriate column before

C compiling the code.

C The user subroutine USDFLD will be called

C at all the integration points of the model

C for which the material definition includes

C user-defined field variables. As a matter

C of example, the elastic gradient profile of

C the second specimen evaluated is

C implemented (FGMII).

SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,

1 T,CELENT,TIME,DTIME,CMNAME,ORNAME,NFIELD,

2 NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,NDI,

3 NSHR,COORD,JMAC,JMATYP,MATLAYO,LACCFLA)

C

INCLUDE ’ABA_PARAM.INC’

C

CHARACTER*80 CMNAME,ORNAME

CHARACTER*3 FLGRAY(15)

DIMENSION FIELD(NFIELD),STATEV(NSTATV),

1 DIRECT(3,3),T(3,3),TIME(2)

DIMENSION ARRAY(15),JARRAY(15),JMAC(*),

1 JMATYP(*),COORD(*)

C In accordance with the procedure adopted

C by ABAQUS, variables are defined implicitly

C and therefore, its type is determined by

C its first letter.

C X and Y are variables of the type real

C that store the coordinates of the

C corresponding integration point relative to

C the horizontal and vertical axis,

C respectively.

X=COORD(1)

Y=COORD(2)

C As seen in Fig. 2, the material elastic

C properties are assumed to vary according

C to a fourth order polynomial function. Note

C that the origin of the coordinate axes has

C been placed following the notation adopted

C in Fig. 2.

C It is also important to note that, in the

C material definition, a direct relation

C between the Young’s Modulus (MPa) and the

C FIELD1 variable is designated so that both

C variables adopt the same values and

C therefore vary equally.

FIELD(1)= - 916659.61963*X**4 +

1 410881.971401*X**3 - 49875.3879242*X**2 -

2 185.636233064*X + 437.202232184

RETURN

END

C In this case, a UFIELD subroutine is used

C in conjunction with the USDFLD subroutine

C defined above in order to implement the

C material elastic properties variation not

C only in the integration points of the

C model, but also in the nodes. This is only

C necessary when computing the stress

C intensity factors in a fracture mechanics

C analysis, since the elastic properties

C in the crack tip must be defined. The

C terminology used is the same as in the

C USDFLD subroutine.

SUBROUTINE UFIELD(FIELD,KFIELD,NSECPT,KSTEP,

1 KINC,TIME,NODE,COORDS,TEMP,DTEMP,NFIELD)

C

INCLUDE ’ABA_PARAM.INC’

C

DIMENSION FIELD(NSECPT,NFIELD), TIME(2),

1 COORDS(3), TEMP(NSECPT), DTEMP(NSECPT)

C

X=COORDS(1)

Y=COORDS(2)

FIELD(1,1)= -916659.61963*X**4 +

1 410881.971401*X**3 - 49875.3879242*X**2 -

2 185.636233064*X + 437.202232184

RETURN

END

## References

- Abanto-Bueno and Lambros (2006) Abanto-Bueno J, Lambros J (2006) An experimental study of mixed mode crack initiation and growth in functionally graded materials. Exp Mech 46(2):179–196
- ABAQUS (2013) ABAQUS (2013) ABAQUS version 6.13 documentation. Providence (RI, USA): Dassault Systèmes Simulia Corp.
- Anderson (2005) Anderson TL (2005) Fracture Mechanics. Fundamentals and Applications. 3rd Ed., Taylor and Francis CRC Press
- Anlas et al (2000) Anlas G, Santare MH, Lambros J (2000) Numerical calculation of stress intensity factors in functionally graded materials. Int J Frac 104(2):131–143
- Bao and Cai (1997) Bao G, Cai H (1997) Delamination cracking in functionally graded coating/metal substrate systems. Acta Mech 45(3):1055–1066
- Bao and Wang (1995) Bao G, Wang L (1995) Multiple cracking in functionally graded ceramic/metal coatings. Int J Solids Struct 32(19):2853–2871
- Becker Jr et al (2001) Becker Jr TL, Cannon RM, Ritchie RO (2001) Finite crack kinking and T-stresses in functionally graded materials. Int J Solids Struct 38(32-33):5545–5563
- Butcher et al (1999) Butcher RJ, Rousseau CE, Tippur HV (1999) A functionally graded particulate composite: preparation, measurements and failure analysis. Acta Mater 47(1):259–68
- Carrillo-Heian et al (2001) Carrillo-Heian EM, Carpenter RD, Paulino GC, Gibeling JG (2001) Dense layered molybdenum disilicide-silicon carbide functionally graded composites formed by field activated synthesis. J Am Ceram Soc 84(5):962–8
- Chapa-Cabrera and Reimanis (2002a) Chapa-Cabrera J, Reimanis IE (2002a) Crack deflection in compositionally graded Cu-W composites. Philos Mag A 82(17-18):3393–3403
- Chapa-Cabrera and Reimanis (2002b) Chapa-Cabrera J, Reimanis IE (2002b) Eï¬ects of residual stress and geometry on crack kink angles in graded composites. Eng Fract Mech 69(14-16):1667–1678
- Chen et al (2000) Chen J, Wu L, Du S (2000) A modified J integral for functionally graded materials. Mechanics Research Communications 27(3):301–306
- Dolbow and Gosz (2002) Dolbow JE, Gosz M (2002) On the computation of mixed-mode stress intensity factors in functionally graded materials. Int J Solids Struct 39(3):2557–2574
- Eischen (1987) Eischen JW (1987) Fracture of nonhomogeneous materials. Int J Frac 34:3–22
- Erdogan (1995) Erdogan F (1995) Fracture mechanics of functionally graded materials. Compos Eng 5(7):753–770
- Gasik (1998) Gasik MM (1998) Micromechanical modelling of functionally graded materials. Compos Mater Sci 13:42–55
- Gu and Asaro (1997) Gu P, Asaro RJ (1997) Crack deflection in functionally graded materials. Int J Solids Struct 34(24):3085–3098
- Gu et al (1997) Gu P, Dao M, Asaro RJ (1997) A simplified method for calculating the crack-tip field of functionally graded materials using the domain integral. J Appl Mech - T ASME 34(1):1–17
- Hashin (1983) Hashin Z (1983) Analysis of composite materials - a survey. J Appl Mech - T ASME 50:481–505
- Healy et al (2012) Healy B, Gullerud A, Koppenhoefer K, Roy A, RoyChowdhury S, et al MW (2012) WARP3D Release 17.3.2 manual. Univ. of Illinois at Urbana-Champaign, Urbana, Ill
- Jedamzik et al (2000) Jedamzik R, Neubrand A, RÃ¶del J (2000) Characterisation of electrochemically processed graded tungsten/copper composites. Mater Sci Forum 308-311:782–7
- Jin and Noda (1994) Jin ZH, Noda N (1994) Crack-tip singular fields in nonhomogeneous materials. J Appl Mech - T ASME 61:738–740
- Kawasaki and Watanabe (1987) Kawasaki A, Watanabe R (1987) Finite element analysis of thermal stress of the metal/ceramic multi-layer composites with compositional gradients. J Jpn Inst Met 51:525–529
- Kim (2003) Kim JH (2003) Mixed-mode crack propagation in functionally graded materials. PhD thesis, University of Illinois at Urbana-Champaign
- Kim and Paulino (2002a) Kim JH, Paulino GH (2002a) Finite element evaluation of mixed-mode stress intensity factors in functionally graded materials. Int J Numer Methods Engrg 53(8):1903–1935
- Kim and Paulino (2002b) Kim JH, Paulino GH (2002b) Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials. J Appl Mech - T ASME 69(4):502–514
- Kim and Paulino (2003) Kim JH, Paulino GH (2003) T-stress, mixed-mode stress intensity factors, and crack initiation angles in functionally graded materials: A unified approach using the interaction integral method. Comput Method Appl M 192(11-12):1463–1494
- Krumova et al (2001) Krumova M, Klingshirn C, Haupert F, Friedrich K (2001) Microhardness studies on functionally graded polymer composites. Compos Sci Technol 61(4):557–63
- Lambros et al (1999) Lambros J, Santare MH, Li H, Sapna III GH (1999) A novel technique for the fabrication of laboratory scale model functionally graded materials. Exp Mech 39(3):184–190
- Li et al (2000) Li H, Lambros J, Cheeseman BA, Santare MH (2000) Experimental investigation of the quasi-static fracture of functionally graded materials. Int J Solids Struct 37(27):3715–3732
- Marur and Tippur (2000) Marur PR, Tippur HV (2000) Numerical analysis of crack-tip fields in functionally graded materials with a crack normal to the elastic gradient. Int J Solids Struct 37(38):5353–5370
- Oral et al (2008a) Oral A, Abanto-Bueno JL, Lambros J, Anlas G (2008a) Crack Initiation Angles in Functionally Graded Materials under Mixed Mode Loading. In: AIP Conference Proceedings, vol 973, pp 248–253
- Oral et al (2008b) Oral A, Lambros J, Anlas G (2008b) Crack Initiation in Functionally Graded Materials Under Mixed Mode Loading: Experiments and Simulations. J Appl Mech - T ASME 75(5):0511,101 (8 pages)
- Parameswaran and Shukla (2000) Parameswaran V, Shukla V (2000) Processing and characterisation of a model functionally graded material. J Mater Sci 35:21–29
- Reiter et al (1997) Reiter T, Dvorak GJ, Tvergaard V (1997) Micromechanical models for graded composite materials. J Mech Phys Solids 45(8):1281–302
- Remmers et al (2008) Remmers JJC, de Borst R, Needleman A (2008) The simulation of dynamic crack propagation using the cohesive segments method. J Mech Phys Solids 56:70–92
- Rice (1968) Rice JR (1968) A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech - T ASME 35:379–86
- Riveiro and Gallego (2013) Riveiro MA, Gallego R (2013) Boundary elements and the analog equation method for the solution of elastic problems in 3-D non-homogeneous bodies. Comput Meth Appl Mech Eng 263:12–19
- Rousseau and Tippur (2000) Rousseau CE, Tippur HV (2000) Compositionally graded materials with cracks normal to the elastic gradient. Acta Mater 48:4021–4033
- Santare and Lambros (2000) Santare MH, Lambros J (2000) Use of graded finite elements to model the behaviour of nonhomogeneous materials. J Appl Mech - T ASME 67:819–822
- Shih and Asaro (1988) Shih CF, Asaro RJ (1988) Elastic-plastic analysis of cracks on bimaterial interfaces: part I â small scale yielding. J Appl Mech - T ASME 58(2):299–316
- Song et al (2006) Song JH, Areias PMA, Belytschko T (2006) A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer Methods Engrg 67:868–893
- Tilbrook et al (2005a) Tilbrook MT, Moon RJ, Hoffman M (2005a) Crack propagation in graded composites. Compos Sci Technol 65:201–220
- Tilbrook et al (2005b) Tilbrook MT, Moon RJ, Hoffman M (2005b) Finite element simulations of crack propagation in functionally graded materials under flexural loading. Engng Fract Mech 72(16):2444–67
- Uemura (2003) Uemura S (2003) The activities of FGM on new application. Mater Sci Forum 423-425:1–10
- Walters et al (2006) Walters MC, Paulino GH, Dodds Jr RH (2006) Computation of mixed-mode stress intensity factors for cracks in three-dimensional functionally graded solids. Journal of Engineering Mechanics 132(1):1–15
- Zhang et al (2003) Zhang C, Savadis A, Savadis G, Zhu H (2003) Transient dynamic analysis of a cracked functionally graded material by a BIEM. Comp Mater Sci 26:167–174