Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 315176, 11 pages

http://dx.doi.org/10.1155/2013/315176

Research Article

Hygrothermal Fracture Analysis of Orthotropic Functionally Graded Materials Using -Integral-Based Methods

Serra Topal¹ and Serkan Dag²

¹Department of Mechanical Engineering, Gazi University, 06570 Ankara, Turkey

²Department of Mechanical Engineering, Middle East Technical University, 06800 Ankara, Turkey

Received 21 April 2013; Accepted 1 August 2013

Academic Editor: Abdelouahed Tounsi

Copyright © 2013 Serra Topal and Serkan Dag. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper puts forward two different -integral-based methods, which can be used to perform mixed-mode fracture analysis of orthotropic functionally graded materials subjected to hygrothermal stresses. The first method requires the evaluation of both components of -integral, whereas the second method employs the first component and the asymptotic crack tip displacement fields. Plane orthotropic hygrothermoelasticity is the basic theory behind the -integral formulation, which is carried out by assuming that all material properties are functions of the spatial coordinates. Developed procedures are implemented by means of the finite element method and integrated into a general purpose finite element analysis software. Temperature and specific moisture concentration fields needed in the fracture analyses are also computed through finite element analysis. Each of the developed methods is utilized in conjunction with the superposition technique to calculate the hygrothermal fracture parameters. An inclined crack located in a hygrothermally loaded orthotropic functionally graded layer is examined in parametric analyses. Comparisons of the results generated by the proposed methods do indicate that both methods lead to numerical results of high accuracy and that the developed form of the -integral is domain independent. Further results are presented so as to illustrate the influences of crack inclination angle, crack length, and crack location upon the modes I and II stress intensity factors.

1. Introduction

Functionally graded materials (FGMs) are heterogeneous composite materials, which were originally proposed to be employed in high temperature applications as protective coatings [1]. However, since then, the concept of introducing gradations in certain physical properties has been realized in a number of other technological applications such as biomedical materials [2], high-performance cutting tools [3], surfaces possessing improved contact-damage resistance [4], and solid oxide fuel cells [5]. The heterogeneity of FGMs stems from the fact that volume fractions of the constituents vary along a particular direction in a predetermined manner. Certain kinds of functionally graded materials are known to be orthotropic in addition to being heterogeneous. For example, FGMs generated by means of plasma spray forming have a lamellar structure and thus contain weak cleavage planes parallel to the bounding planes [6]. Graded materials produced through the use of electron beam physical vapor deposition technique on the other hand have a columnar structure with weak cleavage planes perpendicular to the boundaries [7]. Moreover, fiber-reinforced composites possessing a variable fiber volume fraction can be considered as orthotropic functionally graded materials [8].

Considerable emphasis has been placed on fracture mechanics in research studies pertaining to orthotropic functionally graded materials due to the low fracture toughness of commonly used constituents such as ceramics and plastics. Both analytical and computational methods have been set forth for the purpose of evaluating fracture parameters for orthotropic FGMs. Analytical methods presented are almost invariably based on the approach of singular integral equations [9, 10]. Among the computational methods applied for fracture analysis of orthotropic FGMs, we can mention displacement correlation technique [11], modified crack closure method [12], interaction integral method [13], continuum shape sensitivity technique [14], and the -integral approach [15]. In the studies on fracture analysis of orthotropic FGMs in the literature, the graded medium is assumed to be under either mechanical or thermal loading. However, for certain types of orthotropic functionally graded materials, such as polymer-matrix FGMs with variable fiber spacing, in addition to the mechanical and thermal effects, hygroscopic stresses are also critical. Hygroscopic stresses are induced due to changes in moisture concentration within the polymeric matrix; and when these stresses are sufficiently large they may lead to fracture mechanics-related problems such as cracking, delamination, or fatigue. In many instances, thermal and hygroscopic effects act simultaneously on the structure. In these cases, the combined loading due to changes in temperature and moisture concentration is generally referred to as hygrothermal loading.

This paper presents two different -integral-based computational techniques, which can be used to conduct fracture mechanics analysis of orthotropic functionally graded materials subjected to hygrothermal stresses. The first method developed makes use of both components and of the -integral vector, whereas in the second method the first component is utilized in conjunction with the asymptotic crack tip displacement fields. The formulation of the -integral is carried out by considering the constitutive relations of plane orthotropic hygrothermoelasticity. All material properties are assumed to be functions of the spatial coordinates in the derivations. The formulation yields a domain independent form for the -integral, which comprises area and line integrals. Developed procedures are implemented by means of the finite element method and integrated into the general purpose finite element analysis software ANSYS [16]. Parametric analyses are performed by considering a hygrothermally loaded orthotropic functionally graded layer that contains an inclined edge crack. Comparisons of the hygrothermal fracture parameters computed by the proposed methods indicate that both techniques are capable of producing numerical results of high accuracy. Additional numerical results presented illustrate the influences of factors such as crack length, crack location, and inclination angle on modes I and II stress intensity factors (SIFs).

The organization of the paper is as follows: In Section 2, we outline the formulation of -integral; Section 3 provides the details of the two different fracture analysis methods; finite element analysis techniques used in the implementation are elucidated in Section 4; and numerical results are presented in Section 5. Finally, the paper concludes with Section 6, which contains our final remarks.

2. Formulation of the -Integral for Orthotropic FGMs Subjected to Hygrothermal Loading

Figure 1 depicts an inclined edge crack located in an orthotropic functionally graded material that is subjected to hygrothermal stresses. The medium is assumed to be in a state of either plane stress or strain. and axes are aligned along the principal directions of orthotropy; and the and axes lie parallel and perpendicular to the crack plane, respectively. Thus, in the figure represents the angle of inclination of the edge crack with respect to the principal direction of orthotropy . is an arbitrary area enclosing the crack tip. All material properties are assumed to be functions of the spatial coordinates. Thus, the medium possesses the property of being both orthotropic and inhomogeneous.

Figure 1: An inclined edge crack in an orthotropic functionally graded medium.

Constitutive relations of plane orthotropic hygrothermoelasticity are expressed as follows in the principal coordinate system:

In this equation, and , respectively, stand for strain and stress; and , where is the temperature, is the specific moisture concentration, is the reference temperature, and designates the reference specific moisture concentration. For plane stress, entries of the matrices are given by

and for plane strain

In (2a), (2b), (3a), and (3b), , , and are material parameters of plane orthotropic elasticity, and and are thermal and moisture expansion coefficients, respectively. In general, all material parameters are functions of the spatial coordinates in a graded medium. Moreover, the following relations are valid for an orthotropic material:

Constitutive relations in the crack tip coordinate system comprising the axes and are derived by considering coordinate transformation rules valid for the strain and stress tensors. These constitutive relations are obtained as

The entries of the compliance, thermal expansion, and moisture expansion matrices are provided in the appendix. Note that all of the entries are derived in terms of the inclination angle shown in Figure 1.

The structure of the constitutive relation (5) indicates that the principle of superposition can be used to generate the mixed-mode stress intensity factors corresponding to hygrothermal loading. The results due to thermal loading can be found by assuming , and the results corresponding to hygroscopic loading can be calculated by taking . Utilizing the principle of superposition, these two separate sets of results can be combined to evaluate the results for hygrothermal loading. Thus, it suffices to formulate the problem in terms of and thermal expansion matrix entries , , and . When the results for hygroscopic loading are to be calculated, these quantities need to be replaced by , , , and . The respective sums of the modes I and II stress intensity factors are the final results valid for hygrothermal loading. In what follows, we provide the formulation for thermal loading, which can be used to generate the numerical results for both thermal and hygroscopic loading cases.

For the crack depicted in Figure 1, the -integral is derived to be in the following form:

This equation is valid in the local coordinate system -. here stands for the displacement vector; designates the mechanical strain energy density function; is Kronecker delta; is a piecewise smooth function that is equal to unity at the crack tip and zero on the circumference of area ; denotes the explicit derivative of ; and is the arc length. is the straight line that is initiated at the point where intersects crack faces and terminates at the origin. The outcome of the analysis is independent of the size and shape of area used around the crack tip. In the present study, is specified as a circular area with its center lying at the crack tip. The -function utilized is in the form

where is the radius of area . The integration domains and are depicted in Figure 2.

Figure 2: Integration domains used in the evaluation of the -integral.

Although (6) is to be evaluated in the local coordinate system -, the scalar quantities and can be computed in any coordinate system since they are independent of coordinate transformation. In this study, these two quantities are computed in the principal coordinate system for the sake of simplicity. For plane stress, the expression of in the principal coordinate system is derived as

and for plane strain reads

The explicit derivatives of the mechanical strain energy density function are derived by considering (8), (9a), (9b), and (9c). For the cases of plane stress and strain, these derivatives are given by

Derivatives of with respect to the material parameters and the temperature difference in (10a) and (10b) are found in closed form. The remaining partial derivatives in these equations, that is, the derivatives of the material properties and the temperature difference with respect to spatial coordinates, are calculated numerically during the finite element computations.

3. -Integral-Based Methods

In this section, we elucidate two different -integral-based methods, which can be employed to compute mixed-mode stress intensity factors for orthotropic functionally graded materials under hygrothermal stresses. For both thermal and hygroscopic loading cases, the formulation given in the previous section is applicable provided that the appropriate loading function and material properties are considered. Once mixed-mode stress intensity factors are computed for thermal or hygroscopic loading cases, the resultant SIFs can be obtained by combining the separate stress intensity factors evaluated for these loads. The methods presented in this section allow the evaluation of the mixed-mode SIFs for both thermal and hygroscopic types of loading. The first method makes use of both components of the -integral, and , whereas in the second method we utilize the first component and the asymptotic displacement fields.

3.1. Method I

In this method, we use the relations between the components of the -integral and the mixed-mode stress intensity factors and . Equation (6) indicates that the first component of the -integral comprises solely area integrals, while the second component involves both area and line integrals. The expression of is further simplified by separating the line integral into two parts, one evaluated over a region away from the crack tip and the other evaluated near the crack tip. The integral evaluated over a domain near the crack tip is determined in closed form. Then, the components of the -integral are written as follows:

where is the length of the domain over which the line integral is evaluated analytically. This length is also depicted in Figure 2. The term in (11b) is derived in the form

where is the constant term in the asymptotic expansion of the stress component and is the value of at the crack tip. is given by for plane stress and by for plane strain in the appendix. and are the real and imaginary parts of the two roots of the characteristic equation

whose imaginary parts are positive. in this equation are the crack-tip values of . The expressions of are provided in the appendix for both plane stress and strain.

given by (11a) consists of only area integrals and can be evaluated numerically in a straightforward manner by using the Gauss-quadrature methods. However, the expression of given by (11b) contains , , and , which are unknowns. In order to be able to evaluate , we define a new variable as follows:

For a given loading condition, is calculated for two different values of . Denoting these values by , , the corresponding values by ,   and considering (11b), is expressed as

Once is calculated using (11a) and using (15), mixed-mode stress intensity factors can be evaluated by the following equalities that relate the components of the -integral to and [17]:

where

The nonlinear equation set (16a) and (16b) is solved by employing the Newton-Raphson method. Although in general the outcome of the Newton-Raphson method is sensitive to the initial guesses, it is found that convergence is assured by initially specifying as an arbitrary positive number and as zero.

3.2. Method II

The second method we developed is based on the use of the first component of the -integral and the asymptotic crack tip displacement fields. Referring to Figure 2, asymptotic relative displacements of the crack faces can be written as follows [18]:

where

and and are the roots of the characteristic equation (13) with positive imaginary parts. We define a ratio involving the relative displacements of the crack faces at a radial location in the following form:

From (18a), (18b), and (20) it then follows that

Substituting this result into (16a), is derived as

In this second method, once is determined through (11a), and can be computed by the use of (22) and (21), respectively.

4. Numerical Implementation

The methods described in Section 3 are implemented by means of the finite element method. The proposed procedures are integrated into the general purpose finite element analysis software ANSYS [16]. In thermal fracture analysis, the first step is the determination of the temperature field. In the case of hygroscopic loading, specific moisture concentration distribution needs to be determined at the outset. These fields are computed by solving the governing partial differential equations through the use of finite element method. The governing partial differential equation for the temperature distribution is the heat equation, which is given by

where denotes the temperature; and are the principal coordinates of orthotropy shown in Figure 1; and and are the principal thermal conductivities. The governing partial differential equation for the specific moisture concentration is

in this equation is specific moisture concentration and and are the principal mass diffusivities.

The temperature and specific moisture concentration distributions computed through the solutions of (23) and (24) are used to calculate the components of the -integral and the modes I and II stress intensity factors and for each type of loading, that is, for thermal and hygroscopic loadings. The SIFs generated for these separate loads are then superposed to determine the results valid for hygrothermal loading. Smooth spatial variations of the hygrothermomechanical properties of orthotropic FGMs are taken into account in the finite element analyses by specifying the material properties of each finite element at its centroid. This approach is generally referred to as homogeneous finite element approach and leads to computational results of high accuracy provided that there is an appropriate degree of mesh refinement in the finite element model [11, 19]. The finite element meshes employed in the analyses are constructed by utilizing 6-node triangular elements.

The area and line integrals required to be evaluated in -integral computations are calculated by using the Gauss quadrature in conjunction with the isoparametric finite element concept. Area integrals are computed over circular domains centered at the tip of the crack. The accuracy of the outcome of the numerical procedures developed is influenced by certain parameters needed in the solution. For instance, and values used in (15) and value used in (20) have to be set sufficiently small to evaluate the hygrothermal fracture parameters within a high degree of accuracy. It is found that for both thermal and hygroscopic loading cases, highly accurate results can be obtained by setting , , and , respectively, as , , and , where is the length of the inclined crack.

5. Numerical Results

Hygrothermal fracture analysis methods developed in this study are used to compute the fracture parameters for the problem depicted in Figure 3. The figure illustrates both the geometry and hygrothermal boundary conditions. A graded orthotropic layer of length contains an inclined edge crack whose length is denoted by . The axes and are the principal axes of orthotropy, and the angle of inclination between the crack plane and -axis is symbolized by . We suppose that the medium is in a state of plane strain. and values used in hygrothermal loading are also the reference values for the temperature and the specific moisture concentration. In all analyses and are, respectively, set as 20°C and 0.005. As also mentioned in the previous sections, the modes I and II stress intensity factors for the considered hygrothermal loading are determined by superposing the results calculated for thermal and hygroscopic loads.

Figure 3: Problem geometry and hygrothermal boundary conditions.

The layer is assumed to be a fiber-reinforced composite with a fiber volume fraction decreasing as increases from 0 to . As a result, all of the material properties of the layer become functions of the thickness coordinate . Each of the material property required in the analysis is represented by a power function in the following form:

where and are the values of the material property at and , respectively, and is the exponent of the power function characterizing the nature of the property distribution across the thickness. The properties of the orthotropic FGM layer at and and the symbols designating the exponents are provided in Table 1.

Table 1: Material properties at and and the exponents.

Exponent

(GPa)
80
80

(GPa)
80
40

(GPa)
32
26.7

0.25
0.25

0.25
0.25

0.25
0.3

0.25
0.3

(1/°C)
3(10)⁻⁵
3(10)⁻⁵

(1/°C)
3(10)⁻⁵
6(10)⁻⁵

(1/°C)
3(10)⁻⁵
4.5(10)⁻⁵

0.33
0.33

0.33
0.66

0.33
0.495

(W/m·°C)
0.5
0.5

(W/m·°C)
0.5
1

(m²/s)
3(10) −15
3(10)⁻¹⁵

(m²/s)
3(10) −15
6(10)⁻¹⁵

Temperature and specific moisture concentration fields for the layer are computed by solving (23) and (24) through the finite element method. Normalized crack tip temperature variations are provided in Figure 4. This figure shows the normalized crack tip temperature as functions of the crack inclination angle and normalized crack length . It is seen that crack tip temperature is an increasing function of the inclination angle. This is the expected result since the crack tip gets closer to the boundaries at the higher temperature as becomes larger. The figure also indicates that at a given inclination angle, crack tip temperature is a decreasing function of . The trends observed for the crack tip specific moisture concentration are very similar to those presented in Figure 4.

Figure 4: Normalized crack tip temperature versus inclination angle and . , , and .

In order to be able to verify the developed computational methods, in Table 2 we provide comparisons of the normalized mixed-mode stress intensity factors computed for the hygrothermal loading shown in Figure 3. Normalized stress intensity factors are defined as

The hygrothermal mixed-mode stress intensity factors are calculated by considering four different values of the inclination angle , and for each inclination angle the results are given for four different values of the normalized domain radius . Examining the results, it is seen that the SIFs generated by each method are independent of the domain radius. Hence, it can be deduced that the developed form of the -integral is domain independent. Furthermore, the results obtained by methods I and II are in excellent agreement, which is indicative of the high level of accuracy of the results generated by means of the developed techniques. Deformed shape of the finite element mesh under thermal loading is provided in Figure 5.

Table 2: Comparisons of the mixed-mode stress intensity factors. , , , , , , , and .  

Method I
Method II

0°

0.0098
0.0122
0.0098
0.0119

0.0099
0.0122
0.0098
0.0120

0.0099
0.0123
0.0098
0.0120

0.0099
0.0123
0.0097
0.0119

30°

0.0251
0.0080
0.0250
0.0076

0.0252
0.0080
0.0251
0.0077

0.0253
0.0080
0.0252
0.0077

0.0253
0.0081
0.0252
0.0077

45°

0.0296
−0.0006
0.0296
−0.0008

0.0298
−0.0006
0.0298
−0.0008

0.0300
−0.0006
0.0300
−0.0007

0.0301
−0.0008
0.0301
−0.0007

60°

0.0287
−0.0081
0.0288
−0.0089

0.0292
−0.0083
0.0292
−0.0090

0.0297
−0.0086
0.0296
−0.0091

0.0301
−0.0086
0.0299
−0.0092

Figure 5: Deformed shape of the finite element mesh under thermal loading. , , , , , , , , and .

Further results are provided in Figures 6 and 7, which illustrate the influences of inclination angle, crack length, and crack location upon the mixed-mode stress intensity factors. Since both methods are shown to be capable of producing highly accurate numerical results, it is deemed as sufficient to use method I in the generation of the further results presented in these figures. Normalized domain radius is set as 0.1 in the pertaining calculations.

(a)

(b)

Figure 6: Normalized mixed-mode stress intensity factors versus inclination angle and : (a) mode I SIFs; (b) mode II SIFs. , , , , , , , and .

(a)

(b)

Figure 7: Normalized mixed-mode stress intensity factors versus inclination angle and (a) mode I SIFs; (b) mode II SIFs. , , , , , , , and

Figure 6 shows the variations of the normalized mixed-mode SIFs with respect to for four different values of normalized crack length . It is seen that curves go through maximums at values close to 60°, while attains maximums at values between 0° and 30°. decreases as is increased from to . On the other hand, is found to be an increasing function of except for relatively small values of . We also note that is positive and the crack is open for all considered, whereas can be positive or negative depending on the value of the inclination angle.

In Figure 7, we present plots of normalized modes I and II SIFs as functions of the inclination angle and the relative crack location . is the height of crack mouth and stands for the height of the orthotropic FGM layer. For values approximately less than 45°, is an increasing function of ; that is, increases as is increased from 0.15 to 0.3. When is close to 60°, is not that sensitive to the variations in the relative crack location. For relatively small values of , normalized mode II stress intensity factor is an increasing function of . However, for larger values of the inclination angle drops with a corresponding increase in .

6. Closure

This study sets forth two different -integral-based methods, which can be employed to conduct fracture mechanics analysis of orthotropic functionally graded materials that are under the influence of hygrothermal stresses. In the first of these methods, both components of the -integral are evaluated, whereas the second technique requires the computation of only the first component. Proposed procedures are implemented by means of the finite element method and integrated into a general purpose finite element analysis software. Superposition technique is utilized in the extraction of the hygrothermal modes I and II stress intensity factors.

The comparisons provided do indicate that the derived form of the -integral possesses the required domain independence and that both methods are capable of producing numerical results of high accuracy. Further results are presented to illustrate the impacts of inclination angle, normalized crack length, and relative crack location on the normalized mixed-mode stress intensity factors. In general, the influence of each of these parameters on the fracture behavior is significant. Especially, the effect of the crack inclination angle is seen to be rather pronounced. Under hygrothermal stresses, both and go through maximum values as the inclination angle is increased from zero.

In studies on fracture mechanics and fatigue of orthotropic functionally graded materials, correct evaluation of the mixed-mode stress intensity factors is a basic requirement. Moreover, existence of thermal and hygroscopic effects in the examined problems makes it necessary to incorporate these loadings into the computational framework. The methods presented in this paper are shown to be effective ways of taking into account hygrothermal effects and evaluating fracture mechanics parameters and thus can be used to solve fracture and fatigue problems involving complex geometric configurations and loading conditions.

Appendix

The Entries of the Compliance, Thermal Expansion, and Moisture Expansion Matrices

The elements of the compliance, thermal expansion, and moisture expansion matrices used in (5) are derived by coordinate transformation. For the case of plane stress, these elements are given by

and for plane strain

,  , and in these equalities are derived as follows:

Nomenclature

| Area of the region around the crack tip
| Ratio of the relative displacements of the crack faces
| Crack length
| Elements of the compliance matrix in the crack tip coordinate system
| Specific moisture concentration
| Elements of the compliance matrix in the principal coordinate system
| Mass diffusivities
| Elastic moduli and shear modulus
| Height of the orthotropic functionally graded layer
| Height of the crack mouth
| -integral vector
| Mode I stress intensity factor
| Mode II stress intensity factor
| Thermal conductivities
| Length of the orthotropic functionally graded layer
| Radius of the area around the crack tip
| Radial location at which relative displacements of crack faces are calculated
| Temperature
| Displacement vector
| Mechanical strain energy density function
| Thermal expansion coefficients
| Moisture expansion coefficients
| Length of the region over which the line integral is evaluated analytically
| Total strains in the principal coordinate system
| Total strains in the crack tip coordinate system
| Crack inclination angle
| Poisson’s ratios
| Stress components in the principal coordinate system
| Stress components in the crack tip coordinate system.

Conflict of Interests

The authors declare that there is no conflict of interests in this study.

Acknowledgment

This work was supported by the Scientific and Technological Research Council of Turkey (TUBITAK) through Grant MAG-109M511.

References

N. Noda, “Thermal stresses in functionally graded materials,” Journal of Thermal Stresses, vol. 22, no. 4, pp. 477–512, 1999.

F. Watari, A. Yokoyama, M. Omori et al., “Biocompatibility of materials and development to functionally graded implant for bio-medical application,” Composites Science and Technology, vol. 64, no. 6, pp. 893–908, 2004.

T. Nomura, H. Moriguchi, K. Tsuda, K. Isobe, A. Ikegaya, and K. Moriyama, “Material design method for the functionally graded cemented carbide tool,” International Journal of Refractory Metals and Hard Materials, vol. 17, no. 6, pp. 397–404, 1999.

D. C. Pender, N. P. Padture, A. E. Giannakopoulos, and S. Suresh, “Gradients in elastic modulus for improved contact-damage resistance—part I: The silicon nitride-oxynitride glass system,” Acta Materialia, vol. 49, no. 16, pp. 3255–3262, 2001.

Y. Liu, C. Compson, and M. Liu, “Nanostructured and functionally graded cathodes for intermediate temperature solid oxide fuel cells,” Journal of Power Sources, vol. 138, no. 1-2, pp. 194–198, 2004.

S. Sampath, H. Herman, N. Shimoda, and T. Saito, “Thermal spray processing of FGMs,” MRS Bulletin, vol. 20, no. 1, pp. 27–31, 1995.

W. A. Kaysser and B. Ilschner, “FGM research activities in Europe,” MRS Bulletin, vol. 20, no. 1, pp. 22–26, 1995.

M. A. Benatta, I. Mechab, A. Tounsi, and E. A. Adda Bedia, “Static analysis of functionally graded short beams including warping and shear deformation effects,” Computational Materials Science, vol. 44, no. 2, pp. 765–773, 2008.

S. Dag, B. Yildirim, and F. Erdogan, “Interface crack problems in graded orthotropic media: Analytical and computational approaches,” International Journal of Fracture, vol. 130, no. 1, pp. 471–496, 2004.

S. Dag, B. Yildirim, and D. Sarikaya, “Mixed-mode fracture analysis of orthotropic functionally graded materials under mechanical and thermal loads,” International Journal of Solids and Structures, vol. 44, no. 24, pp. 7816–7840, 2007.

S. Dag and A. A. Ilhan, “Mixed-mode fracture analysis of orthotropic functionally graded material coatings using analytical and computational methods,” Journal of Applied Mechanics, Transactions ASME, vol. 75, no. 5, Article ID 051104, 2008.

J.-H. Kim and G. H. Paulino, “Mixed-mode fracture of orthotropic functionally graded materials using finite elements and the modified crack closure method,” Engineering Fracture Mechanics, vol. 69, no. 14–16, pp. 1557–1586, 2002.

J.-H. Kim and G. H. Paulino, “The interaction integral for fracture of orthotropic functionally graded materials: Evaluation of stress intensity factors,” International Journal of Solids and Structures, vol. 40, no. 15, pp. 3967–4001, 2003.

B. N. Rao and S. Rahman, “A continuum shape sensitivity method for fracture analysis of orthotropic functionally graded materials,” Mechanics of Materials, vol. 37, no. 10, pp. 1007–1025, 2005.

S. Dag, E. Erhan Arman, and B. Yildirim, “Computation of thermal fracture parameters for orthotropic functionally graded materials using Jk-integral,” International Journal of Solids and Structures, vol. 47, no. 25-26, pp. 3480–3488, 2010.

ANSYS, ANSYS Basic Analysis Procedures Guide, Release 5.4, ANSYS Inc., Canonsburg, Pa, USA, 1997.

J.-H. Kim and G. H. Paulino, “Mixed-mode J-integral formulation and implementation using graded elements for fracture analysis of nonhomogeneous orthotropic materials,” Mechanics of Materials, vol. 35, pp. 107–128, 2003.

P. Sollero and M. H. Aliabadi, “Fracture mechanics analysis of anisotropic plates by the boundary element method,” International Journal of Fracture, vol. 64, no. 4, pp. 269–284, 1993.

B. Yildirim, S. Dag, and F. Erdogan, “Three dimensional fracture analysis of FGM coatings under thermomechanical loading,” International Journal of Fracture, vol. 132, no. 4, pp. 369–395, 2005.