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

Authors: Serra Topal; Serkan Dag
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper puts forward two differentJk-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 Jk-integral, whereas the second method employs the first component J1 and the asymptotic crack tip displacement fields. Plane orthotropic hygrothermoelasticity is the basic theory behind the Jk-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 Jk-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.


---

## Body

## 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 Jk-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 differentJk-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 J1 and J2 of the Jk-integral vector, whereas in the second method the first component J1 is utilized in conjunction with the asymptotic crack tip displacement fields. The formulation of the Jk-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 Jk-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 Section2, we outline the formulation of Jk-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 theJk-Integral for Orthotropic FGMs Subjected to Hygrothermal Loading

Figure1 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. x and y axes are aligned along the principal directions of orthotropy; and the x1 and x2 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 x. A 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:(1)[εxxεyy2εxy]=[c11c120c12c22000c66][σxxσyyσxy]+[d1d20]ΔT+[e1e20]Δc.
In this equation, εij and σij, respectively, stand for strain and stress; ΔT=T-Tr and Δc=c-cr, where T is the temperature, c is the specific moisture concentration, Tr is the reference temperature, and cr designates the reference specific moisture concentration. For plane stress, entries of the matrices are given by (2a)c11=1Ex,c12=-νxyEx,c22=1Ey,c66=1Gxy,(2b)d1=αx,d2=αy,e1=βx,e2=βy, and for plane strain(3a)c11=1-νzxνxzEx,c12=-νxy+νxzνzyEx,c22=1-νyzνzyEy,c66=1Gxy,(3b)d1=νzxαz+αx,d2=νzyαz+αy,e1=νzxβz+βx,e2=νzyβz+βy.In (2a), (2b), (3a), and (3b), Ei, νij, and Gxy are material parameters of plane orthotropic elasticity, and αi and βi 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:
(4)νyxEy=νxyEx,νzxEz=νxzEx,νzyEz=νyzEy.Constitutive relations in the crack tip coordinate system comprising the axesx1 and x2 are derived by considering coordinate transformation rules valid for the strain and stress tensors. These constitutive relations are obtained as
(5)[ε11ε222ε12]=[a11a12a16a12a22a26a12a26a66][σ11σ22σ12]+[g1g2g12]ΔT+[h1h2h12]Δc.
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  Δc=0, and the results corresponding to hygroscopic loading can be calculated by taking ΔT=0. 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 ΔT and thermal expansion matrix entries g1, g2, and g12. When the results for hygroscopic loading are to be calculated, these quantities need to be replaced by Δc, h1, h2, and h12. 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 Figure1, the Jk-integral is derived to be in the following form:
(6)Jk=∬A{σijui,k-Wδkj}q,jdA-∬A(W,k)explqdA-∫Γc(W+-W-)δ2kqds,(i,j,k=1,2).
This equation is valid in the local coordinate system x1-x2.ui here stands for the displacement vector; W designates the mechanical strain energy density function; δkj is Kronecker delta; q is a piecewise smooth function that is equal to unity at the crack tip and zero on the circumference of area A; (W,k)expl denotes the explicit derivative of W; and s is the arc length. Γc is the straight line that is initiated at the point where A intersects crack faces and terminates at the origin. The outcome of the analysis is independent of the size and shape of area A used around the crack tip. In the present study, A is specified as a circular area with its center lying at the crack tip. The q-function utilized is in the form
(7)q(x1,x2)=1-x12+x22R,
where R is the radius of area A. The integration domains A and Γc are depicted in Figure 2.Figure 2
Integration domains used in the evaluation of theJk-integral.Although (6) is to be evaluated in the local coordinate system x1-x2, the scalar quantities W and (W,k)expl 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 W in the principal coordinate system is derived as
(8)W=Ex2(εxx-αxΔT)+νxyExEy(εyy-αyΔT)2(Ex-νxy2Ey)×(εxx-αxΔT)+2Gxyεxy2+νxyExEy(εxx-αxΔT)+ExEy(εyy-αyΔT)2(Ex-νxy2Ey)×(εyy-αyΔT),
and for plane strain W reads(9a)W=Exε-xx2D×{(1-Eyνxzνzy2Exνzx)Exε-xx+(νxy+νxzνzy)Eyε-yy}+2Gxyεxy2+ExEyε-yy2D×{(νxy+νxzνzy)ε-xx+(1-νzxνxz)ε-yy}+Exνzx2νxzαz2(ΔT)2,(9b)D=Ex(1-νzxνxz)-Ey(νzy2νxzνzx+νxy2+2νxyνxzνzy),(9c)ε-xx=εxx-(νzxαz+αx)ΔT,ε-yy=εyy-(νzyαz+αy)ΔT.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(10a)(∂W∂xk)expl=∂W∂Ex∂Ex∂xk+∂W∂Ey∂Ey∂xk+∂W∂νxy∂νxy∂xk+∂W∂Gxy∂Gxy∂xk+∂W∂αx∂αx∂xk+∂W∂αy∂αy∂xk+∂W∂(ΔT)∂(ΔT)∂xk,ddddd(k=1,2)forplanestress,(10b)(∂W∂xk)expl=∂W∂Ex∂Ex∂xk+∂W∂Ey∂Ey∂xk+∂W∂νzy∂νzy∂xk+∂W∂νxz∂νxz∂xk+∂W∂νxy∂νxy∂xk+∂W∂νzy∂νzy∂xk+∂W∂Gxy∂Gxy∂xk+∂W∂αx∂αx∂xk+∂W∂αy∂αy∂xk+∂W∂αz∂αz∂xk+∂W∂(ΔT)∂(ΔT)∂xk,ddddd(k=1,2)forplanestrain.Derivatives of W 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.Jk-Integral-Based Methods

In this section, we elucidate two differentJk-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 (i.e.,ΔTorΔc) 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 Jk-integral, J1 and J2, whereas in the second method we utilize the first component J1 and the asymptotic displacement fields.


### 3.1. Method I

In this method, we use the relations between the components of theJk-integral and the mixed-mode stress intensity factors KI and KII. Equation (6) indicates that the first component of the Jk-integral comprises solely area integrals, while the second component J2 involves both area and line integrals. The expression of J2 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 Jk-integral are written as follows:(11a)J1=∬A{σijui,1-Wδ1j}q,jdA-∬A(W,1)explqdA(i,j=1,2),(11b)J2=∬A{σijui,2-Wδ2j}q,jdA-∬A(W,2)explqdA-∫0R-δ(W+-W-)qds-2δπC(i,j=1,2),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 C in (11b) is derived in the form
(12)C=-2a11tipσ110{(α1η2+α2η1)KI+(η1+η2)KII},
where σ110 is the constant term in the asymptotic expansion of the stress component σ11 and a11tip is the value of a11 at the crack tip. a11 is given by (A.1a) for plane stress and by (A.2a) for plane strain in the appendix. αi and ηi are the real and imaginary parts of the two roots of the characteristic equation
(13)a11tips4-2a16tips3+(2a12tip+a66tip)s2-2a26tips+a22tip=0,
whose imaginary parts are positive. aijtip in this equation are the crack-tip values of aij. The expressions of aij are provided in the appendix for both plane stress and strain.J
1 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 J2 given by (11b) contains σ110, KI, and KII, which are unknowns. In order to be able to evaluate J2, we define a new variable J^2 as follows:
(14)J^2=∬A{σijui,2-Wδ2j}q,jdA-∬A(W,2)explqdA-∫0R-δ(W+-W-)qds.
For a given loading condition, J^2 is calculated for two different values of δ. Denoting these δ values by δ1, δ2, the corresponding J^2 values by J^21,  J^22 and considering (11b), J2 is expressed as
(15)J2=1δ2/π-δ1/π{δ2πJ^21-δ1πJ^22}.
Once J1 is calculated using (11a) and J2 using (15), mixed-mode stress intensity factors can be evaluated by the  following equalities that relate the components of the Jk-integral to KI and KII [17]:(16a)J1=a11tip2{c11KI2+c12KIKII+c22KII2},(16b)J2=-a11tip4{d11KI2+d12KIKII+d22KII2},where(17a)c11=α12η2+α22η1+η12η2+η1η22,(17b)c12=2(α1η2+α2η1),(17c)c22=η1+η2,(17d)d11=α1η2+α2η1,(17e)d12=4{η12η2+η1η22-α1α2(η1+η2)},(17f)d22=2{α12η2+α22η1+α1α2(η1+η2)}-α1η2-α2η1.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 KI as an arbitrary positive number and KII as zero.


### 3.2. Method II

The second method we developed is based on the use of the first component of theJk-integral J1 and the asymptotic crack tip displacement fields. Referring to Figure 2, asymptotic relative displacements of the crack faces can be written as follows [18]:(18a)u1(r,π)-u1(r,-π)=22rπ{D11KI+D12KII},(18b)u2(r,π)-u2(r,-π)=22rπ{D21KI+D22KII},where(19a)D11=Im{s2P11-s1P12s1-s2},D12=Im{P11-P12s1-s2},(19b)D21=Im{s2P21-s1P22s1-s2},D22=Im{P21-P22s1-s2},(19c)P1i=a11tipsi2+a12tip-a16tipsi(i=1,2),(19d)P2i=a12tipsi+a22tipsi-a26tip(i=1,2),and s1 and s2 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 ra in the following form:
(20)Ar=u2(ra,π)-u2(ra,-π)u1(ra,π)-u1(ra,-π).
From (18a), (18b), and (20) it then follows that
(21)KIKII=ArD12-D22D21-ArD11=B.
Substituting this result into (16a), KI is derived as
(22)KI=2J1a11tip(c11+(c12/B)+(c22/B2)).
In this second method, once J1 is determined through (11a), KI and KII can be computed by the use of (22) and (21), respectively.


## 3.1. Method I

In this method, we use the relations between the components of theJk-integral and the mixed-mode stress intensity factors KI and KII. Equation (6) indicates that the first component of the Jk-integral comprises solely area integrals, while the second component J2 involves both area and line integrals. The expression of J2 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 Jk-integral are written as follows:(11a)J1=∬A{σijui,1-Wδ1j}q,jdA-∬A(W,1)explqdA(i,j=1,2),(11b)J2=∬A{σijui,2-Wδ2j}q,jdA-∬A(W,2)explqdA-∫0R-δ(W+-W-)qds-2δπC(i,j=1,2),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 C in (11b) is derived in the form
(12)C=-2a11tipσ110{(α1η2+α2η1)KI+(η1+η2)KII},
where σ110 is the constant term in the asymptotic expansion of the stress component σ11 and a11tip is the value of a11 at the crack tip. a11 is given by (A.1a) for plane stress and by (A.2a) for plane strain in the appendix. αi and ηi are the real and imaginary parts of the two roots of the characteristic equation
(13)a11tips4-2a16tips3+(2a12tip+a66tip)s2-2a26tips+a22tip=0,
whose imaginary parts are positive. aijtip in this equation are the crack-tip values of aij. The expressions of aij are provided in the appendix for both plane stress and strain.J
1 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 J2 given by (11b) contains σ110, KI, and KII, which are unknowns. In order to be able to evaluate J2, we define a new variable J^2 as follows:
(14)J^2=∬A{σijui,2-Wδ2j}q,jdA-∬A(W,2)explqdA-∫0R-δ(W+-W-)qds.
For a given loading condition, J^2 is calculated for two different values of δ. Denoting these δ values by δ1, δ2, the corresponding J^2 values by J^21,  J^22 and considering (11b), J2 is expressed as
(15)J2=1δ2/π-δ1/π{δ2πJ^21-δ1πJ^22}.
Once J1 is calculated using (11a) and J2 using (15), mixed-mode stress intensity factors can be evaluated by the  following equalities that relate the components of the Jk-integral to KI and KII [17]:(16a)J1=a11tip2{c11KI2+c12KIKII+c22KII2},(16b)J2=-a11tip4{d11KI2+d12KIKII+d22KII2},where(17a)c11=α12η2+α22η1+η12η2+η1η22,(17b)c12=2(α1η2+α2η1),(17c)c22=η1+η2,(17d)d11=α1η2+α2η1,(17e)d12=4{η12η2+η1η22-α1α2(η1+η2)},(17f)d22=2{α12η2+α22η1+α1α2(η1+η2)}-α1η2-α2η1.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 KI as an arbitrary positive number and KII as zero.


## 3.2. Method II

The second method we developed is based on the use of the first component of theJk-integral J1 and the asymptotic crack tip displacement fields. Referring to Figure 2, asymptotic relative displacements of the crack faces can be written as follows [18]:(18a)u1(r,π)-u1(r,-π)=22rπ{D11KI+D12KII},(18b)u2(r,π)-u2(r,-π)=22rπ{D21KI+D22KII},where(19a)D11=Im{s2P11-s1P12s1-s2},D12=Im{P11-P12s1-s2},(19b)D21=Im{s2P21-s1P22s1-s2},D22=Im{P21-P22s1-s2},(19c)P1i=a11tipsi2+a12tip-a16tipsi(i=1,2),(19d)P2i=a12tipsi+a22tipsi-a26tip(i=1,2),and s1 and s2 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 ra in the following form:
(20)Ar=u2(ra,π)-u2(ra,-π)u1(ra,π)-u1(ra,-π).
From (18a), (18b), and (20) it then follows that
(21)KIKII=ArD12-D22D21-ArD11=B.
Substituting this result into (16a), KI is derived as
(22)KI=2J1a11tip(c11+(c12/B)+(c22/B2)).
In this second method, once J1 is determined through (11a), KI and KII can be computed by the use of (22) and (21), respectively.


## 4. Numerical Implementation

The methods described in Section3 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
(23)∂∂x(kx(x,y)∂T∂x)+∂∂y(ky(x,y)∂T∂y)=0,
where T denotes the temperature; x and y are the principal coordinates of orthotropy shown in Figure 1; and kx and ky are the principal thermal conductivities. The governing partial differential equation for the specific moisture concentration is
(24)∂∂x(Dx(x,y)∂c∂x)+∂∂y(Dy(x,y)∂c∂y)=0.c in this equation is specific moisture concentration and Dx and Dy 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 Jk-integral and the modes I and II stress intensity factors KI and KII 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 inJk-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, δ1 and δ2 values used in (15) and ra 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 δ1, δ2, and ra, respectively, as a/2000, a/1000, and a/8000, where a 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 Figure3. The figure illustrates both the geometry and hygrothermal boundary conditions. A graded orthotropic layer of length L contains an inclined edge crack whose length is denoted by a. The axes x and y are the principal axes of orthotropy, and the angle of inclination between the crack plane and x-axis is symbolized by ϕ. We suppose that the medium is in a state of plane strain. T0 and c0 values used in hygrothermal loading are also the reference values for the temperature and the specific moisture concentration. In all analyses T0 and c0 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 h. 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:
(25)M(ξ)=M0+(Mh-M0)(ξh)p,
where M0 and Mh are the values of the material property at ξ=0 and ξ=h, respectively, and p 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 ξ=0 and ξ=h and the symbols designating the exponents are provided in Table 1.Table 1
Material properties atξ=0 and ξ=h and the exponents.
ξ
=
0
ξ
=
h
Exponent
E
x
(GPa)
80
80
γ
x
E
y
(GPa)
80
40
γ
y
G
x
y
(GPa)
32
26.7
γ
x
y
ν
x
y
0.25
0.25
ω
x
y
ν
x
z
0.25
0.25
ω
x
z
ν
z
x
0.25
0.3
ω
z
x
ν
z
y
0.25
0.3
ω
z
y
α
x
(1/°C)
3(10)−5
3(10)−5
λ
x
α
y
(1/°C)
3(10)−5
6(10)−5
λ
y
α
z
(1/°C)
3(10)−5
4.5(10)−5
λ
z
β
x
0.33
0.33
κ
x
β
y
0.33
0.66
κ
y
β
z
0.33
0.495
κ
z
k
x
(W/m·°C)
0.5
0.5
η
x
k
y
(W/m·°C)
0.5
1
η
y
D
x
(m2/s)
3(10)−15
3(10)−15
ρ
x
D
y
(m2/s)
3(10)−15
6(10)−15
ρ
yTemperature 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 a/L. 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 a/L. 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 a/L. ηx=ηy=2.5, hc/h=0.25, and h/L=0.4.In order to be able to verify the developed computational methods, in Table2 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
(26)KIn=KIα(0)Ex(0)T0πa,KIIn=KIIα(0)Ex(0)T0πa.
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 R/a. 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 Jk-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.γx=γy=γxy=2, ωxy=ωxz=ωzx=ωzy=2, λx=λy=λz=2.5, κx=κy=κz=1.5, ηx=ηy=2.5, ρx=ρy=1.5, hc/h=0.25,h/L=0.4,  and  a/L=0.1.
ϕ
R
/
a
Method I
Method II
K
I
n
K
II
n
K
I
n
K
II
n
0°
0.1
0.0098
0.0122
0.0098
0.0119
0.2
0.0099
0.0122
0.0098
0.0120
0.3
0.0099
0.0123
0.0098
0.0120
0.4
0.0099
0.0123
0.0097
0.0119
30°
0.1
0.0251
0.0080
0.0250
0.0076
0.2
0.0252
0.0080
0.0251
0.0077
0.3
0.0253
0.0080
0.0252
0.0077
0.4
0.0253
0.0081
0.0252
0.0077
45°
0.1
0.0296
−0.0006
0.0296
−0.0008
0.2
0.0298
−0.0006
0.0298
−0.0008
0.3
0.0300
−0.0006
0.0300
−0.0007
0.4
0.0301
−0.0008
0.0301
−0.0007
60°
0.1
0.0287
−0.0081
0.0288
−0.0089
0.2
0.0292
−0.0083
0.0292
−0.0090
0.3
0.0297
−0.0086
0.0296
−0.0091
0.4
0.0301
−0.0086
0.0299
−0.0092Figure 5
Deformed shape of the finite element mesh under thermal loading.γx=γy=γxy=2, ωxy=ωxz=ωzx=ωzy=2, λx=λy=λz=2.5, κx=κy=κz=1.5, ηx=ηy=2.5, ρx=ρy=1.5, hc/h=0.25, h/L=0.4, and a/L=0.1.Further results are provided in Figures6 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 R/a is set as 0.1 in the pertaining calculations.Normalized mixed-mode stress intensity factors versus inclination angleϕ and a/L: (a) mode I SIFs; (b) mode II SIFs. γx=γy=γxy=2, ωxy=ωxz=ωzx=ωzy=2, λx=λy=λz=2, κx=κy=κz=2, ηx=ηy=2, ρx=ρy=2, hc/h=0.25, and h/L=0.4.
(a)
(b)Normalized mixed-mode stress intensity factors versus inclination angleϕ and hc/h: (a) mode I SIFs; (b) mode II SIFs. γx=γy=γxy=2, ωxy=ωxz=ωzx=ωzy=2, λx=λy=λz=2, κx=κy=κz=2, ηx=ηy=2, ρx=ρy=2, h/L=0.4, and a/L=0.1.
(a)
(b)Figure6 shows the variations of the normalized mixed-mode SIFs with respect to ϕ for four different values of normalized crack length a/L. It is seen that KIn curves go through maximums at ϕ values close to 60°, while KIIn attains maximums at ϕ values between 0° and 30°. KIn decreases as a/L is increased from 1/12 to 1/8. On the other hand, KIIn is found to be an increasing function of a/L except for relatively small values of ϕ. We also note that KIn is positive and the crack is open for all ϕ considered, whereas KIIn can be positive or negative depending on the value of the inclination angle.In Figure7, we present plots of normalized modes I and II SIFs as functions of the inclination angle ϕ and the relative crack location hc/h. hc is the height of crack mouth and h stands for the height of the orthotropic FGM layer. For ϕ values approximately less than 45°, KIn is an increasing function of hc/h; that is, KIn increases as hc/h is increased from 0.15 to 0.3. When ϕ is close to 60°, KIn 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 hc/h. However, for larger values of the inclination angle KIIn drops with a corresponding increase in hc/h.


## 6. Closure

This study sets forth two differentJk-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 Jk-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 theJk-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 KI and KII 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.



---

*Source: 315176-2013-09-14.xml*