Fatigue Fracture of Functionally Graded Materials Under Elastic-Plastic Loading Conditions Using Extended Finite Element Method

Abstract

In this chapter, extended finite element method (XFEM) has been used to simulate the fatigue crack growth problems in functionally graded material (FGM) in the presence of hole, inclusion and minor crack under elastic and plastic conditions. The fatigue crack growth analysis of alloy/ceramic FGMs, alloy and equivalent composite is done by XFEM in the presence of multiple discontinuities under mode-I mechanical load. The validity of linear elastic fracture mechanics (LEFM) theory is limited to the brittle materials. Therefore, the elastic plastic fracture mechanics (EPFM) theory needs to be utilized to characterize the plastic behavior of the material. A generalized Ramberg-Osgood material model has been used to model the stress-strain behavior of the material. Plasticity has been checked by Von Mises Yield criteria. J-integral has been used to calculate the SIF. Crack growth direction is determined by maximum principal stress criteria.

1. Introduction

Development of novel materials improves performance and efficiency of the structures, and also leads to development of advanced and sophisticated structures. This complex process of materials, structures and technology has led to the development of composite materials. Strength and stiffness plays a key role in evaluating the worth of the material. These characteristics provide strength to the structure to retain its desired shape and size under loading or any other external action.

Cracks/flaws are inevitable in all engineering materials. Loading under severe environmental conditions may either initiate new cracks or may cause the propagation of pre-existing cracks in the structures. Theoretically, fracture can be defined as the breaking or rupturing of a material resulting into its separation into two or more pieces.

Composite materials manifested in the middle of the twentieth century. Composites are naturally occurring or engineered materials made from two or more constituents with different chemical or physical properties distinct boundary among constituents. Lightweight composite materials with high strength to weight and stiffness to weight ratios have been used successfully in aircraft industry and other engineering applications. Under high temperature conditions the strength of the metal is deteriorated whereas, ceramics have excellent resistance to heat.

FGMs can be referred as multiphase composite materials in which the composition or microstructure or both are spatially varied which lead to a certain gradation in the local material properties. FGMs can be defined as multi-phase composites. FGMs are synthesized such that they own continuous variations in volume fractions of their components in space to return a pre-established composition. FGMs possess continuously varying properties in one or more than one direction and the form non-homogeneous macrostructure due to these variations. By gradually varying the volume fraction of the constituents, FGMs exhibit a smooth and continuous change from one surface to another, thus reducing interface problems, and minimizing thermal stress concentrations. The ceramic phase of FGMs provides a good resistance to heat, while the metal phase provides a strong mechanical performance and hence reduces the possibility of catastrophic failure.

The major advantages of FGM over conventional materials are firstly, FGM satisfies the working conditions for which it is specifically developed. Secondly, it is economical as it reduces material costs for particular engineering applications. Thirdly, it can reduce the magnitude of residual and thermal stresses generated under working conditions. Finally, FGMs exhibit better fracture toughness and bond strength. This is normally achieved by using a ceramic layer connected with a metallic layer. FGMs have wide area of engineering applications like in the computer circuit and aerospace industries. FGMs have typical applications is in aircraft and automotive industries as thermal barrier coatings (TBCs).

In general, all structural components are subjected to thermo-mechanical cyclic load. The fatigue life of these components is generally predicted without considering the effect of defects/discontinuities present in component. However, FGMs are commonly made by sintering process, which are porous in nature. These discontinuities at the vicinity of a major crack tip lead to increase the effective SIF at the major crack tip due to which the life of the components get depreciated. Hence, the analysis of FGMs in the vicinity of discontinuities becomes very important from the design point of view. To widen the spectrum of applications of FGMs, the fatigue/fracture behavior should be properly evaluated.

Over the years, greater understanding of fracture mechanics has undoubtedly prevented a significant number of structural failures. Fracture mechanics approach for the design of structures includes flaw size as one of the key variables. Fracture toughness replaces strength of material as a relevant material attribute, and its evaluation is mainly done in composites using the $J$-integral approach [1]. Failure of FGM has always been a trending domain of research for scientists and engineers due to the wide spectra of their engineering applications.

Fatigue and quasi-static fracture are two forms of crack growth phenomenon. Fatigue fracture refers to the slow propagation of cracks under cyclic loading conditions where the stress intensity factors are below the fracture toughness of the material. Quasi-static fracture is observed near the end of the fatigue life when the increased crack length leads to stress intensity factors which are above the fracture toughness [2, 3].

In many cases, multiple cracks may exist in the components. There interaction resulting in the variation of stress intensity factor, stress distribution and propagation direction of the major crack. In the past, the failure of structures was analyzed in the presence of multiple cracks [4]. Some efforts have been made using analytical, experimental and simulation techniques to analyze the effect of interaction among multiple cracks [5, 6, 7].

Although, many analytical [8, 9] and experimental methods [10, 11] have been explored for the calculation of fracture parameters even then the drawbacks associated with experimental investigation and scarcity of analytical solution have impelled the analysts towards alternative techniques. Numerical methods hold the promise in this regard.

Many numerical methods are available to simulate the problems of fatigue failure in materials. These include finite element method (FEM), boundary element method, hybrid boundary node method [12, 13], meshfree methods [14, 15, 16, 17, 18] and extended finite element method [19, 20]. Out of these methods, FEM has been widely used for solving a wide variety of engineering and industrial problems [21, 22, 23, 24, 25, 26]. It has achieved a remarkable success in solving various linear and non-linear problems [27, 28, 29, 30, 31, 32, 33, 34]. Despite its numerous advantages and unparalleled success, it is not well-suited for solving the problems involving crack propagation. In crack growth problems, element edges provide natural lines along which a crack can grow. This is advantageous if the crack path is known a priori, but in most of the fracture phenomenon, the crack path is unknown. Thus, FEM requires a conformal mesh and re-meshing to ensure that the element boundaries coincide with the moving discontinuities (crack). Moreover, crack tip singularity cannot be accurately modeled by standard finite element approximation. Therefore, the modeling of crack growth becomes quite tedious and time consuming due to the modification in mesh topology at each stage of crack propagation. To overcome this difficulty, a new method known as extended finite element method (XFEM) has been developed to model arbitrary discontinuities without a need of conformal mesh or re-meshing. Level set method (LSM) is used in conjunction with XFEM for defining as well as tracking the geometry of cracks and other discontinuities like holes and inclusions. To cope up with these problems, XFEM has been adopted as a tool for the analysis of fatigue crack propagation in FGM.

2. Calculation of SIF for FGM

A domain based interaction integral approach can be used for calculating the stress intensity factors for homogeneous, bi-layer and functionally graded materials under thermal as well as mechanical loading. In this chapter, interaction integral approach will be extended to calculate the SIFs for FGM and bi-layered FGM under mechanical loads. The interaction integral is calculated based on J-integral. The J-integral for an elastic body subjected to thermo-mechanical load is given as,

$$J=\underset{{\mathrm{\Gamma }}_o}{\mathrm{\oint }}\left(\tilde{W}{\delta }_{1j}-{\sigma }_{\mathit{ij}}\frac{\partial u_i}{\partial x_1}\right)n_{j}d\mathrm{\Gamma }$$E1

For the interaction integral calculation of an elastic body, consider two equilibrium states i.e. state 1, the actual state with given boundary conditions and state 2, an auxiliary state of the cracked body. The parameters for auxiliary state are represented with superscript $a$. The final expression for the interaction integral takes the form [35]

$$\begin{array}{l}{M}_{12}=\underset{A_o}{\int }\left({\sigma }_{\mathit{ij}}\frac{\partial u_i^a}{\partial x_1}+{\sigma }_{\mathit{ij}}^a\frac{\partial u_i}{\partial x_1}-{\sigma }_{\mathit{ik}}^a{\epsilon }_{\mathit{ik}}^m{\delta }_{1j}\right)\frac{\partial q}{\partial x_j}\mathit{dA}+\\ \underset{A_o}{\int }\left({\sigma }_{\mathit{ij}}\left(S_{\mathit{\text{ijkl}}}^{\mathit{tip}}-S_{\mathit{\text{ijkl}}}\left(\mathbf{x}\right)\right)\frac{\partial {\sigma }_{\mathit{kl}}^a}{\partial x_1}\right)q\mathit{dA}\end{array}$$E2

where, the auxiliary field for the FGM may be taken from [36] as.

$${\sigma }_{\mathit{ij}}^a=C_{\mathit{\text{ijkl}}}^{\mathit{tip}}\frac{1}{2}\left(\frac{\partial u_k^a}{\partial x_l}+\frac{\partial u_l^a}{\partial x_k}\right),{\epsilon }_{\mathit{ij}}^a=S_{\mathit{\text{ijkl}}}\left(\mathbf{x}\right){\sigma }_{\mathit{kl}}^a\text{and}{\epsilon }_{\mathit{ij}}^a\ne \frac{1}{2}\left(\frac{\partial u_i^a}{\partial x_j}+\frac{\partial u_j^a}{\partial x_i}\right)$$E3

The SIFs are calculated from the interaction integral as [36]:

Mode-I SIF is given as,

$$K_I=\frac{M_{12}{E}^{\ast }{\text{cosh}}^2\left({\mathit{\pi \epsilon }}^{\mathit{tip}}\right)}{2}\text{with}{K}_I^a=1\text{and}{K}_{\mathit{II}}^a=0$$E4a

Mode-II SIF is given as,

$$K_{\mathit{II}}=\frac{M_{12}{E}^{\ast }{\text{cosh}}^2\left({\mathit{\pi \epsilon }}^{\mathit{tip}}\right)}{2}\text{with}{K}_I^a=0\text{and}{K}_{\mathit{II}}^a=1$$E4b

where, $E^{\ast }=\frac{2{\overline{E}}_1{\overline{E}}_2}{{\overline{E}}_1+{\overline{E}}_2}$ with ${\overline{E}}_i=\left\{\begin{array}{l}{E}_i^{\mathit{tip}}\text{for plane stress}\\ E_i^{\mathit{tip}}/\left(1-{\left({\nu }_i^{\mathit{tip}}\right)}^2\right)\text{for plane strain}\end{array}\right.$ with $i=1,2$

3. Fatigue crack growth

Here we use Paris law for stable crack propagation, the generalized Paris’s law is given as:

$$\frac{\mathit{da}}{\mathit{dN}}=C{\left(\Delta K_{\mathit{Ieq}}\right)}^m$$E5

Where, $a$ is the crack length and $N$ is the number of loading cycles. $C$ and $m$ are material properties to find the rate of crack growth. At each crack tip, the local direction of crack growth ${\theta }_c$ can be calculated by the maximum principal stress theory [37]. Crack is assumed to grow in a direction perpendicular to the maximum principal stress. Thus, by enforcing the condition that the local shear stress is zero for $\theta ={\theta }_c$,

$$K_I\text{sin}{\theta }_c+K_{\mathit{II}}\left(3\text{cos}{\theta }_c-1\right)=0$$E6

The solution of Eq. (6) gives

$${\theta }_c=2{\text{tan}}^{-1}\left(\frac{K_I-\sqrt{K_I^2+8K_{\mathit{II}}^2}}{4K_{\mathit{II}}}\right)$$E7

According to this criterion, the equivalent mode-I SIF is obtained as

$$K_{\mathit{Ieq}}=K_I{\text{cos}}^3\left(\frac{{\theta }_c}{2}\right)-3K_{\mathit{II}}{\text{cos}}^2\left(\frac{{\theta }_c}{2}\right)\text{sin}\left(\frac{{\theta }_c}{2}\right)$$E8

For stable crack propagation, the generalized Paris’ law for FGM is given as

$$\frac{\mathit{da}}{\mathit{dN}}=C\left(\mathbf{x}\right){\left(\mathrm{\Delta }{K}_{\mathit{Ieq}}\right)}^{m\left(\mathbf{x}\right)}$$E9

where, $C\left(\mathbf{x}\right)$ and $m\left(\mathbf{x}\right)$ are the functions of the location.

In the numerical example, the crack growth value Δa is assumed and the corresponding number of cycles ΔN is computed from Eq. (9). When multiple crack tips are present, the crack growth value Δa is assumed for the most dominant crack tip, corresponding ΔN is computed and then at the other crack tips the crack growth is computed corresponding to ΔN. Eventually, when the maximum value of $K_{\mathit{Ieq}}$ for any crack tip becomes more than the fracture toughness $K_{\mathit{IC}}$ at corresponding location then the simulation is terminated. At this point, the total number of cycles elapsed is the fatigue life of the FGM.

4. Modeling of the properties of FGM

In this chapter, the results have been presented for a FGM plate as shown in Figure 1. The FGM plate is manufactured by reinforcing an alloy with ceramic. The volume fraction of ceramic is varied in the $x$-direction to get a material property variation in the $x$-direction. It is assumed that at $x=0$ the FGM have the properties of the alloy and at $x=L$ properties of ceramic. The major crack is always taken at the center of the FGM plate in the $x$-direction. The interface, when present is also in the same direction. The material properties of the aluminum alloy and alumina used in FGM are tabulated in Table 1 [38, 39].

Figure 1.

Geometry of the FGM plate along with its dimensions.

Material properties	Aluminum alloy	Alumina
Elastic modulus $E\left(\mathrm{GPa}\right)$	70	300
Poisson’s ratio, $\nu$	0.33	0.21
Coefficient of thermal expansion $\gamma$ (/^°C)	25 × 10^−6	8.2 × 10^−6
Fracture toughness $K_{\mathit{IC}}\left(\mathrm{MPa}\sqrt{m}\right)$	29	3.5
Paris law parameter $C$ in $m/\mathit{\text{cycle}}{\left(\mathrm{MPa}\sqrt{m}\right)}^{-m}$	10^−12	2.8 × 10^−10
Paris law parameter, $m\left(\mathbf{x}\right)$	3	10

Table 1.

Material properties of aluminum alloy and alumina.

The variation of the elastic modulus for FGM is modeled as

$$E\left(x\right)=E_{\mathit{\text{alloy}}}{e}^{\alpha x}\text{where}\alpha \text{is\hspace{0.17em}given}\mathrm{as}\alpha =\frac{1}{L}\text{ln}\left(\frac{E_{\mathit{\text{ceramic}}}}{E_{\mathit{\text{alloy}}}}\right)$$E10

A plot of $E\left(\mathbf{x}\right)$ for L = 100 mm is shown in Figure 2. The fatigue life of FGM has been compared with the same of the aluminum alloy and an equivalent composite of aluminum alloy/alumina. The equivalent composite considered in this example has the same overall volume fractions of aluminum alloy and ceramic as the FGM. The volume fractions of ceramic and aluminum alloy in the FGM are obtained as

$$V_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)=\frac{E\left(x\right)-E_{\mathit{\text{alloy}}}}{E_{\mathit{\text{ceramic}}}-E_{\mathit{\text{alloy}}}}=\frac{E_{\mathit{\text{alloy}}}{e}^{\alpha x}-E_{\mathit{\text{alloy}}}}{E_{\mathit{\text{ceramic}}}-E_{\mathit{\text{alloy}}}}$$E11a

$$V_{\mathit{\text{alloy}}}^{\mathit{FGM}}\left(x\right)=1-V_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)$$E11b

Figure 2.

Variation of modulus of elasticity along the length of the plate.

In this example, the equivalent composite is assumed to have the same amount of metal and ceramic. The volume fraction of alumina in the equivalent composite is calculated as

$$V_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}=\frac{1}{L}\underset{0}{\overset{L}{\int }}{V}_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)\mathit{dx}$$E11c

where, $L$ is the length of the plate. For $L=100\mathrm{mm}$, $V_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}=38.28$% and $V_{\mathit{\text{alloy}}}^{\mathit{\text{composite}}}=61.72$%. The variation in volume fraction of ceramic (alumina) in the FGM is shown in Figure 3. The volume fraction for the equivalent composite has also been indicated. Now, using the rule of mixtures for the equivalent composite

$$E_{\mathit{\text{composite}}}=E_{\mathit{\text{alloy}}}{V}_{\mathit{\text{alloy}}}^{\mathit{\text{composite}}}+E_{\mathit{\text{ceramic}}}{V}_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}$$E12

we get $E_{\mathit{\text{composite}}}=158.04\mathrm{GPa}$. The Poisson’s ratio for the equivalent composite as well as for the FGM may be calculated as [40]

$$\nu \left(x\right)=\frac{{\nu }_{\mathit{\text{alloy}}}{V}_{\mathit{\text{alloy}}}^{\mathit{FGM}}\left(x\right)E_{\mathit{\text{ceramic}}}+{\nu }_{\mathit{\text{ceramic}}}{V}_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)E_{\mathit{\text{alloy}}}}{V_{\mathit{\text{alloy}}}^{\mathit{FGM}}\left(x\right)E_{\mathit{\text{ceramic}}}+V_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)E_{\mathit{\text{alloy}}}}$$E13a

$$\text{and}{\nu }_{\mathit{\text{composite}}}=\frac{{\nu }_{\mathit{\text{alloy}}}{V}_{\mathit{\text{alloy}}}^{\mathit{\text{composite}}}{E}_{\mathit{\text{ceramic}}}+{\nu }_{\mathit{\text{ceramic}}}{V}_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}{E}_{\mathit{\text{alloy}}}}{V_{\mathit{\text{alloy}}}^{\mathit{\text{composite}}}{E}_{\mathit{\text{ceramic}}}+V_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}{E}_{\mathit{\text{alloy}}}}$$E13b

Figure 3.

Variation of volume fraction of ceramic along the length of the plate.

The Poisson’s ratio is shown in Figure 4.

Figure 4.

Variation of Poisson’s ratio along the length of the plate.

The coefficient of thermal expansion for the FGM by the rule of mixtures is calculated as

$$\gamma \left(x\right)={\gamma }_{\mathit{\text{alloy}}}{V}_{\mathit{\text{alloy}}}^{\mathit{FGM}}\left(x\right)+{\gamma }_{\mathit{\text{ceramic}}}{V}_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)$$E14a

The coefficient of thermal expansion for the equivalent composite is given by

$${\gamma }_{\mathit{\text{composite}}}={\gamma }_{\mathit{\text{alloy}}}{V}_{\mathit{\text{alloy}}}^{\mathit{\text{composite}}}+{\gamma }_{\mathit{\text{ceramic}}}{V}_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}$$E14b

The value of the coefficient of thermal expansion for the equivalent composite is calculated using Eq. (14b) and is found to be ${\gamma }_{\mathit{\text{composite}}}$ = 18.57 × 10⁻⁶/^°C. A variation of coefficient of thermal expansion for the FGM is shown in Figure 5.

Figure 5.

Variation of coefficient of thermal expansion along the length of the plate.

The fracture toughness of the FGM as well as the equivalent composite may be expressed as a function of the volume fraction of the ceramic by the following formula given by [41]

$$\begin{array}{l}{K}_{\mathit{IC}}\left(x\right)=\frac{K_{\mathit{IC}}^{\mathit{\text{alloy}}}+K_{\mathit{IC}}^{\mathit{\text{ceramic}}}}{2}\\ +\frac{K_{\mathit{IC}}^{\mathit{\text{alloy}}}-K_{\mathit{IC}}^{\mathit{\text{ceramic}}}}{2}\left(\sqrt{1-V_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)}-\sqrt{V_{\mathit{\text{ceramic}}}^{\mathit{FGM}}\left(x\right)}\right)\end{array}$$E15a

$$\begin{array}{l}{K}_{\mathit{IC}}^{\mathit{\text{composite}}}=\frac{K_{\mathit{IC}}^{\mathit{\text{alloy}}}+K_{\mathit{IC}}^{\mathit{\text{ceramic}}}}{2}\\ +\frac{K_{\mathit{IC}}^{\mathit{\text{alloy}}}-K_{\mathit{IC}}^{\mathit{\text{ceramic}}}}{2}\left(\sqrt{1-V_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}}-\sqrt{V_{\mathit{\text{ceramic}}}^{\mathit{\text{composite}}}}\right)\end{array}$$E15b

The variation of fracture toughness is shown in Figure 6.

Figure 6.

Variation of fracture toughness along the length of the plate.

The Paris law parameters are assumed to have exponential variation in a manner similar to the elastic modulus. Thus, the variation in the parameters of Paris equation is taken as

$$C\left(x\right)=C_{\mathit{\text{alloy}}}{e}^{\vartheta x},\text{where},\vartheta =\frac{1}{L}\text{ln}\left(\frac{C_{\mathit{\text{ceramic}}}}{C_{\mathit{\text{alloy}}}}\right)$$E16

$$m\left(x\right)=m_{\mathit{\text{alloy}}}{e}^{\varsigma x},\text{where},\varsigma =\frac{1}{L}\text{ln}\left(\frac{m_{\mathit{\text{ceramic}}}}{m_{\mathit{\text{alloy}}}}\right)$$E17

For the equivalent composite, we find the location at which the volume fraction of ceramic in the FGM is same as that of the equivalent composite. This location $\overline{x}$ may be found by either Figure 5 or by using the formula $\overline{x}=\frac{1}{\alpha }\text{ln}\left(\frac{E_{\mathit{\text{composite}}}}{E_{\mathit{\text{alloy}}}}\right)$, where $\alpha$ is defined in Eq. (10). For the present example $\overline{x}=56\mathrm{mm}$. The Paris law parameters of the equivalent composite is assumed to be same as that of the FGM at $x=\overline{x}$. Thus,

$$C_{\mathit{\text{composite}}}=C_{\mathit{\text{alloy}}}{e}^{\vartheta \overline{x}}$$E18a

$$m_{\mathit{\text{composite}}}=m_{\mathit{\text{alloy}}}{e}^{\varsigma \overline{x}}$$E18b

The values of $C$ and $m$ for the equivalent composite comes out to be $C_{\mathit{\text{composite}}}=2.34\times {10}^{-11}\mathrm{m}/\text{cycle}{\left(\mathrm{MPa}\sqrt{\mathrm{m}}\right)}^{-\mathrm{m}}$ and $m_{\mathit{\text{composite}}}=5.88$. The variation of C and $m$ are shown in Figures 7 and 8 respectively.

Figure 7.

Variation of $C\left(\mathbf{x}\right)$ along the length of the plate.

Figure 8.

Variation of $m\left(\mathbf{x}\right)$ along the length of the plate.

Plastic behaviour for FGM can be modeled using Ramberg Osgood equation [42]

$$\epsilon =\frac{\sigma }{E}+{\left(\frac{\sigma }{H}\right)}^{1/n}$$E19

Here, $H$ is the strength coefficient and $n$ is the strain hardening exponent. The value of $n=0.0946$ is used for the present example. The values of the parameters of Paris equation are taken as $C=3\times {10}^{-11}$ and $m=3$. In actual case the path of crack growth is curved but in this study the linear crack growth path is taken. Linear crack extension length $\mathit{\Delta a}$ for an edge crack is kept constant. For a center crack maximum crack extension length $\mathrm{\Delta }{a}_{\mathit{max}}$ is kept on principal crack tip. The principle crack tip is the crack tip where ${\mathrm{\Delta }K}_{\mathit{Ieq}}$ maximum. Crack increment at the other crack tip is given by:

$$\mathrm{\Delta }a=\mathrm{\Delta }{a}_{\mathit{max}}{\left(\frac{{\mathit{\Delta K}}_{\mathit{Ieq}}}{{\mathit{\Delta K}}_{\mathit{Ieq}\mathit{max}}}\right)}^m$$E20

The crack tip extension at the principal crack tip is $\mathrm{\Delta }{a}_{\mathit{max}}$ and at the other crack tip extension is smaller. The crack extension takes place $K_{\mathit{Ieq}\mathit{max}}<K_{\mathit{IC}}$ . Crack becomes unstable when $K_{\mathit{Ieq}\mathit{max}}>K_{\mathit{IC}}$. Simulation continues until this condition is met. Here, $K_{\mathit{Ieq}\mathit{max}}$ is the equivalent SIF for mode-I at principal crack tip and $K_{\mathit{IC}}$ is the material property called fracture toughness or critical SIF. $K_{\mathit{IC}}$ for FGM is given by [35]

$$K_{\mathit{IC}}\left(x\right)=K_{\mathit{IC}}^{\mathit{\text{ceramic}}}{\left[\frac{E\left(x\right)}{1-{\nu }_{\mathit{FGM}}^2}\left\{V_m\left(x\right)\frac{1-{\nu }_{\mathit{\text{alloy}}}^2}{E_{\mathit{\text{alloy}}}}{\left(\frac{K_{\mathit{IC}}^{\mathit{\text{alloy}}}}{K_{\mathit{IC}}^{\mathit{\text{ceram}}\mathit{ic}}}\right)}^2+\left(1-V_m\left(x\right)\right)\frac{1-{\nu }_{\mathit{\text{ceramic}}}^2}{E_{\mathit{\text{ceramic}}}}\right\}\right]}^{1/2}$$E21

where $K_{\mathit{IC}}\left(x\right)$ is the fracture toughness of the FGM at point x. $K_{\mathit{IC}}^{\mathit{\text{alloy}}}$ and $K_{\mathit{IC}}^{\mathit{\text{ceramic}}}$ are the fracture toughness of the alloy and ceramic, while νalloy and νceramic are Poisson’s ratios for the alloy and ceramic respectively. $V_m\left(x\right)$ denotes the volume fraction for the alloy at point $x$.

The constitutive relation for the elastic-plastic material is given as

$$\sigma \left(u\right)=D_{\mathit{ep}}\left(\mathrm{x}\right)\epsilon \left(u\right)$$E22

where x is the vector of x and y-coordinates, $D_{\mathit{ep}}\left(\mathrm{x}\right)$ is elastic-plastic constitutive matrix varying in x-direction. The elastic constitutive matrix can be written for plane stress condition as

$$D_e\left(\mathbf{x}\right)=\frac{E\left(\mathbf{x}\right)}{\left\{1-\nu {\left(\mathbf{x}\right)}^2\right\}}\left[\begin{array}{lll}1& \nu \left(\mathbf{x}\right)& 0\\ \nu \left(\mathbf{x}\right)& 1& 0\\ 0& 0& \frac{1-\nu \left(\mathbf{x}\right)}{2}\end{array}\right]=D$$E23a

and for plane strain condition as

$$D_e\left(\mathbf{x}\right)=\frac{E\left(\mathbf{x}\right)}{\left\{1-2\nu \left(\mathbf{x}\right)\right\}\left\{1+\nu \left(\mathbf{x}\right)\right\}}\left[\begin{array}{lll}1-\nu \left(\mathbf{x}\right)& \nu \left(\mathbf{x}\right)& 0\\ \nu \left(\mathbf{x}\right)& 1-\nu \left(\mathbf{x}\right)& 0\\ 0& 0& \frac{1-2\nu \left(\mathbf{x}\right)}{2}\end{array}\right]=D$$E23b

The incremental theory of plasticity [43] has been used to model the elastic-plastic constitutive relation for a material. An incremental stress vector $\mathit{d\sigma }$ and incremental strain vector $\mathit{d\epsilon }$ are such that $\mathit{d\sigma }=D_{\mathit{ep}}.\mathit{d\epsilon }$. Where $D_{\mathit{ep}}$ is the elastic-plastic constitutive matrix, which is determined as discussed under:

Total strain increment is the sum of elastic and plastic strains

$$\mathit{d\epsilon }=d{\epsilon }_e+{\mathit{d\epsilon }}_p$$E24

Elastic incremental strain and stress is determined

$$\mathit{d\sigma }=D_{e}d{\epsilon }_e$$E25

$$F\left(\sigma \right)=f\left(\overline{\sigma }\right)$$E26

where $\sigma$ is the stress tensor and $\overline{\sigma }$ is the equivalent stress, $F$ and $f$ are two different failure functions. By the flow rule the incremental strain is related to the gradient of a function known plastic potential. If the plastic potential function and the failure function is same, then the following relation is obtained,

$${\mathit{d\epsilon }}_p=\nabla F.d\lambda$$E27

Plastic modulus ${\rm H}$ is given as

$${\rm H}=\frac{d\overline{\sigma }}{d{\overline{\epsilon }}_p}$$E28

For a given strain energy $\mathit{\delta w}$, and according to the definitionof ${\mathit{d\epsilon }}_p$ we must have,

$$\mathit{\delta w}=\overline{\sigma }.d{\overline{\epsilon }}_p$$E29

According to the Von Mises criteria, $F=J_2$, where $J_2$ is the second invariant of deviatoric stress tensor. So, we must have $F\left(\overline{\sigma }\right)=\frac{{\overline{\sigma }}^2}{3}$, Thus Eqs. (26) and (27) result in

$$\mathit{d\sigma }=D_e\left\{\mathit{d\epsilon }-d{\epsilon }_p\right\}$$E30

After taking the derivatives from both sides of failure criteria equation

$$\left(\frac{\partial F}{\partial \sigma }\mathit{d\sigma }\right)=\left(\frac{\partial f}{\partial \overline{\sigma }}.\frac{\partial \overline{\sigma }}{\partial {\overline{\epsilon }}_p}.\frac{\partial {\overline{\epsilon }}_p}{\partial w}.\frac{\partial w}{\partial {\overline{\epsilon }}_p}.\partial {\epsilon }_p\right)$$E31

For simplicity we take $\frac{\partial F}{\partial \sigma }=a$, $\frac{\partial f}{\partial \overline{\sigma }}=\overline{a}$

$$a.\mathit{d\sigma }=\overline{a}.{\rm H}.\left(\frac{1}{\sigma }\right)\sigma .{\mathit{d\epsilon }}_p$$E32

$\mathit{d\lambda }$ is calculated by omitting $\mathit{d\sigma }$ between Eqs. (30) and (31) and substituting ${\mathit{d\epsilon }}_p$ from Eq. (27). By substituting $\mathit{d\lambda }$ in Eq. (27), the final form of material matrix is obtained as [43]

$$D_{\mathit{ep}}=D_e-D_p$$E33

$$\text{where}{D}_p=\frac{\mathit{Da} a^{T}D}{\frac{\overline{a}}{\overline{\sigma }}{\rm H}{\sigma }^T+a^T\mathit{Da}}$$E34

5. XFEM: Introduction and formulation for cracks and discontinuities in FGM

XFEM or the extended finite element method is a numerical technique which allows crack modeling irrespective of the mesh, and eliminates the cumbersome process of remeshing in problems involving change in the crack geometry like crack growth. XFEM models a crack by enriching the standard finite element approximation with some functions, which are obtained from the theoretical background of the problem. Moving discontinuities are tracked by the level set method. XFEM is a numerical method, based on the finite element method (FEM) that is especially designed for treating discontinuities. The formulation is done as discussed in [35, 44]. The solution of FGM differs from homogeneous materials only in the spatial gradation in the material properties. After calculating the values of stress and strain, the SIF is determined.

6. Numerical examples and discussion

The FGM plate considered in all the numerical simulations has 100% aluminum alloy on one side and 100% alumina on the other side. The volume fraction of alumina changes from 0% on one side to 100% on the other side so as to produce an FGM. The equivalent composite is equivalent to the FGM in the sense that both the FGM and the composite plate contain the same amount of aluminum alloy and alumina. The fatigue crack growth analysis of alloy/ceramic FGMs, aluminum alloy and equivalent composite is done by XFEM in the presence of multiple cracks, holes and inclusions under mode-I mechanical load and their fatigue life are compared. The constituents of the FGM plate are aluminum alloy and alumina. A major crack of large initial length is assumed to exist at the edge of the plate. The major crack is assumed to be in the direction of material gradation. The fatigue crack growth analyses of the FGM, the equivalent composite and the aluminum alloy plates have been carried out in the presence of minor cracks, holes and inclusions till the final failure of the plate under mode-I mechanical load. The effect of these small defects on the fatigue life as well as on the crack path has been investigated in detail.

6.1. Example 1

A rectangular FGM plate of length ($L$) 100 mm. and height ($D$) 200 mm. with 100% aluminum alloy on left side and 100% ceramic (alumina) on right side is considered. Property variation is taken in x-direction, where x = 0 to x = 100 mm. The plate with a major edge crack of length a = 20 mm is analyzed under plane strain condition in the presence of multiple discontinuities. In all simulations, the plate dimensions, initial crack length and material properties are taken to be same. The properties of FGM, composites and aluminum alloy are already described in Table 1. The material properties of the inclusions are taken as ${\rm E}=20\mathrm{GPa}$ and $\nu =0.2$. The plate domain is discretized using uniformly distributed 117 nodes in x-direction and 235 nodes in y-direction. The fatigue crack growth analysis is performed by taking a crack increment of $\mathrm{\Delta }a=\frac{a}{10}$= 2 mm. A cyclic tensile load varying from ${\sigma }_{\max }=70\mathrm{MPa}$ to ${\sigma }_{\min }=0\mathrm{MPa}$ is applied in all the simulations. The geometric discontinuities like holes, inclusions and minor cracks are added in the plate in addition to the major edge or center crack to analyze their effect on the fatigue life of the material. The fatigue life of the FGM, equivalent composite and aluminum alloy are obtained under mode-I loading, and are compared with each other.

6.2. Plate with a major edge crack under linear elastic condition

Figures 9 and 10 show a plate with a major edge crack of length a = 20 mm at the left and right edge respectively. These plates have been analyzed under plane strain condition using a uniform mesh of 117 by 235 nodes. The plots of the fatigue life for different materials are shown in Figure 11. From these figures, it is seen that the equivalent composite withstands 7885 cycles before it fails while the FGM with crack on alloy side undergoes 15,561 cycles and pure aluminum alloy undergoes 19,145 cycles before failure. It is also observed that when a major crack initiates from the ceramic (alumina) rich side then it fails much earlier (4872 cycles) as compared to when the crack initiates from the aluminum alloy side.

Figure 9.

Plate with an edge crack on the alloy rich side under mode-I loading.

Figure 10.

Plate with an edge crack on the ceramic rich side under mode-I loading.

Figure 11.

A plot of crack extension with number of cycles.

These plots show that when a crack is present on the ceramic rich side, the life diminishes by a considerable extent as compared to when a crack is present on the alloy rich side. The equivalent composite shows the minimum life except in case when a crack is present on the ceramic side. It is also observed that the crack follows nearly a straight path in all the materials.

6.3. Plate with a major edge crack, minor cracks, holes and inclusions under linear elastic condition

In this case, a major crack of length a = 20 mm. is taken at the left and the right edge of the plate ($100\times 200\mathrm{mm}$) as shown in Figures 12 and 13 respectively. In addition to the major edge crack, 36 minor cracks, 15 holes and 15 inclusions are randomly distributed in the plate. The length of the minor cracks varies from 3.5 to 4.5 mm, and orientation varies from 0 to 60^° randomly. The holes and inclusions have variations in their radii from 3 to 4.5 mm randomly. A cyclic mode-I mechanical load is applied at the top edge of the plate. The plots for crack extension with number of cycles are shown in Figure 14.

Figure 12.

Plate with an edge crack on the alloy rich side under mode-I loading.

Figure 13.

Plate with an edge crack on the ceramic rich side under mode-I loading.

Figure 14.

A plot of crack extension with number of cycles.

It is also observed that the crack deflects in all the materials. Moreover, it is seen that the number of cycles to failure in case of aluminum alloy is about 18,111 cycles whereas in case of FGM with crack on the alloy and ceramic rich sides is 14,622 cycles and 3111 cycles respectively. The fatigue life of the composite plate is found to 6956 cycles. Thus, it can be stated that due to the presence of minor cracks, holes and inclusions, the life of the aluminum alloy is reduced by about 5.42%, whereas the fatigue life of the FGM with crack on the alloy and ceramic rich sides goes down by 6.03 and 36.15% respectively. The fatigue life of the equivalent composite is reduced by 11.78%.

6.4. Example 2

A rectangular FGM plate of length ($\mathit{Lt}$) 100 mm. and height ($\mathit{Ht}$) 200 mm. with 100% copper nickel alloy on left side and 100% ceramic (alumina) on right side is considered. Property variation is taken in x-direction, where x varies from x = 0 to x = 100 mm. A uniform traction of 100 MPa is applied on the top edge of the rectangular domain along y direction. Cyclic loading is applied at top edge of the plate with a maximum value of ${\sigma }_{\mathit{max}}=100\mathrm{MPa}$ and minimum value of ${\sigma }_{\mathit{min}}=0\mathrm{MPa}$. A uniform mesh of size 117 × 235 nodes is used for the analysis in each case. The values of SIFs are computed at the tip of the major crack. The variation of SIF with crack length is plotted in each case. The material properties are taken from Table 2 [43].

Material properties	Values
Elastic modulus of copper nickel alloy $E_{\mathit{\text{alloy}}}\left(\mathrm{GPa}\right)$	160
Elastic modulus of alumina (ceramic) $E_{\mathit{\text{ceramic}}}\left(\mathrm{GPa}\right)$	386
Elastic modulus of soft inclusion $E_{\mathit{\text{inclusion}}}\left(\mathrm{GPa}\right)$	100
Elastic modulus of Hard inclusion $E_{\mathit{\text{inclusion}}}\left(\mathrm{GPa}\right)$	400
Poisson’s ratio of copper nickel alloy ${\nu }_{\mathit{\text{alloy}}}$	0.35
Poisson’s ratio of alumina (ceramic) ${\nu }_{\mathit{\text{ceramic}}}$	0.21
Poisson’s ratio of inclusion ${\nu }_{\mathit{\text{inclusion}}}$	0.3
Poisson’s ratio of inclusion ${\nu }_{\mathit{FGM}}$	0.23
Fracture toughness of copper nickel alloy $K_{\mathit{IC}}^{\mathit{\text{alloy}}}\left(\mathrm{MPa}\sqrt{m}\right)$	79
Fracture toughness of alumina (ceramic) $K_{\mathit{IC}}^{\mathit{\text{ceramic}}}\left(\mathrm{MPa}\sqrt{m}\right)$	5
Paris constant $C$ in $m/\mathit{\text{cycle}}{\left(\mathrm{MPa}\sqrt{m}\right)}^{-m}$	3 × 10^−11
Paris exponent $m$	3

Table 2.

Material property table.

6.5. A major crack in FGM plate under elastic: Plastic loading condition

In this case, a major crack of length $a=20\mathrm{mm}$ is taken at the edge of the domain (100 × 200 mm) as shown in Figure 15. Cyclic loading is applied at the top edge of the FGM plate, and a crack propagates due to this loading. The plots of SIF with crack length for an crack configuration is shown in Figure 16. The failure crack length obtained for edge crack is 0.0402 m.

Figure 15.

FGM plate with an edge crack.

Figure 16.

Plot for variation of SIF with crack length.

6.6. A major edge crack in FGM plate with holes, inclusions and minor cracks under elastic-Plastic loading condition

In this case, a major crack of length $a=20$ mm is taken at the edge of the domain (100 × 200 mm) is taken as shown in Figure 17. Minor cracks, holes and inclusions are randomly distributed in the plate. All 36 minor cracks have varying length randomly from 3.5 to 4.5 mm, with varying orientation from 0 to 60°. In addition to these 15 inclusions are also distributed in the domain randomly. The holes and inclusions have variation in their radii from 3 to 4.5 mm. A cyclic mode-I loading is applied due to which the major crack propagates. The plots for SIF variation with crack length of edge crack is shown in Figures 18 and 19 for soft and hard inclusions respectively. The failure crack length for edge crack is obtained 0.0384 and 0.0392 m. for soft and hard inclusions respectively.

Figure 17.

FGM plate with an edge crack, 15 inclusions, 15 holes and 36 minor cracks.

Figure 18.

Plot for variation of SIF with crack length for soft inclusions.

Figure 19.

Plot for variation of SIF with crack length for hard inclusions.

7. Conclusions

In this chapter we have discussed the simulation of cracks in a FGM plate has been carried out in the presence of multiple inhomogeneities by XFEM using both linear elastic as well as elastic-plastic formulations. SIF has been calculated at the tip of the major crack using interaction integral approach. The variation in the SIF at the tip of the major crack has been studied when multiple inhomogeneities are present in the domain. From this study it is observed that minor cracks have least effect in the FGM plate’s failure crack length, whereas soft inclusions have moderate effect and holes have the most severe effect. It is found that the FGM plate’s life increases in each case when soft inclusions are replaced by hard inclusions. Hence the presence of the hard inclusions in the plate increases the failure crack length of the plate i.e. plate survives more.

Nomenclatures

W˜:	Strain energy
σ:	Stress
ε:	Strain
u:	Deformation
C:	Compliance matrix
E:	Modulus of elasticity
γ:	Coefficient of thermal expansion
KIC:	Critical stress intensity factor (Fracture toughness)