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

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: 𝜀_𝑥𝑥𝜀_𝑦𝑦2𝜀_𝑥𝑦=𝑐₁₁𝑐₁₂0𝑐₁₂𝑐₂₂000𝑐₆₆𝜎_𝑥𝑥𝜎_𝑦𝑦𝜎_𝑥𝑦+𝑑₁𝑑₂0𝑒Δ𝑇+𝑒¹₂0Δ𝑐. (1) 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 𝑐₁₁=1𝐸_𝑥,𝑐₁₂𝜈=−^𝑥𝑦𝐸_𝑥,𝑐₂₂=1𝐸_𝑦,𝑐₆₆=1𝐺_𝑥𝑦,𝑑₁=𝛼_𝑥,𝑑₂=𝛼_𝑦,𝑒₁=𝛽_𝑥,𝑒₂=𝛽_𝑦, (2a)(2b) and for plane strain𝑐₁₁=1−𝜈_𝑧𝑥𝜈_𝑥𝑧𝐸_𝑥,𝑐₁₂𝜈=−^𝑥𝑦+𝜈_𝑥𝑧𝜈_𝑧𝑦𝐸_𝑥,𝑐₂₂=1−𝜈_𝑦𝑧𝜈_𝑧𝑦𝐸_𝑦,𝑐₆₆=1𝐺_𝑥𝑦,𝑑₁=𝜈_𝑧𝑥𝛼_𝑧+𝛼_𝑥,𝑑₂=𝜈_𝑧𝑦𝛼_𝑧+𝛼_𝑦,𝑒₁=𝜈_𝑧𝑥𝛽_𝑧+𝛽_𝑥,𝑒₂=𝜈_𝑧𝑦𝛽_𝑧+𝛽_𝑦. (3a)(3b)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: 𝜈_𝑦𝑥𝐸_𝑦=𝜈_𝑥𝑦𝐸_𝑥,𝜈_𝑧𝑥𝐸_𝑧=𝜈_𝑥𝑧𝐸_𝑥,𝜈_𝑧𝑦𝐸_𝑧=𝜈_𝑦𝑧𝐸_𝑦. (4)

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 𝜀₁₁𝜀₂₂2𝜀₁₂=𝑎₁₁𝑎₁₂𝑎₁₆𝑎₁₂𝑎₂₂𝑎₂₆𝑎₁₂𝑎₂₆𝑎₆₆𝜎₁₁𝜎₂₂𝜎₁₂+𝑔₁𝑔₂𝑔₁₂ℎΔ𝑇+ℎ¹₂ℎ₁₂Δ𝑐. (5) 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 Δ𝑐=0, and the results corresponding to hygroscopic loading can be calculated by taking Δ𝑇=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 Δ𝑇 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: 𝐽_𝑘=_𝐴𝜎_𝑖𝑗𝑢_𝑖,𝑘−𝑊𝛿_𝑘𝑗𝑞_,𝑗𝑑𝐴−_𝐴𝑊_,𝑘expl−𝑞𝑑𝐴_Γ𝑐𝑊⁺−𝑊⁻𝛿_2𝑘𝑞𝑑𝑠,𝑖,𝑗,𝑘=1,2. (6) 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 𝐴; (𝑊_,𝑘)_expl 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 𝑞𝑥₁,𝑥₂=1−𝑥²₁+𝑥²₂𝑅, (7) 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 (𝑊_,𝑘)_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 𝑊 in the principal coordinate system is derived as 𝐸𝑊=^2𝑥𝜀_𝑥𝑥−𝛼_𝑥Δ𝑇+𝜈_𝑥𝑦𝐸_𝑥𝐸_𝑦𝜀_𝑦𝑦−𝛼_𝑦Δ𝑇2𝐸_𝑥−𝜈²_𝑥𝑦𝐸_𝑦×𝜀_𝑥𝑥−𝛼_𝑥Δ𝑇+2𝐺_𝑥𝑦𝜀²_𝑥𝑦+𝜈_𝑥𝑦𝐸_𝑥𝐸_𝑦𝜀_𝑥𝑥−𝛼_𝑥Δ𝑇+𝐸_𝑥𝐸_𝑦𝜀_𝑦𝑦−𝛼_𝑦Δ𝑇2𝐸_𝑥−𝜈²_𝑥𝑦𝐸_𝑦×𝜀_𝑦𝑦−𝛼_𝑦Δ𝑇, (8) and for plane strain 𝑊 reads𝐸𝑊=^𝑥𝜀_𝑥𝑥×𝐸2𝐷1−^𝑦𝜈_𝑥𝑧𝜈²_𝑧𝑦𝐸_𝑥𝜈_𝑧𝑥𝐸_𝑥𝜀_𝑥𝑥+𝜈_𝑥𝑦+𝜈_𝑥𝑧𝜈_𝑧𝑦𝐸_𝑦𝜀_𝑦𝑦+2𝐺_𝑥𝑦𝜀²_𝑥𝑦+𝐸_𝑥𝐸_𝑦𝜀_𝑦𝑦×𝜈2𝐷_𝑥𝑦+𝜈_𝑥𝑧𝜈_𝑧𝑦𝜀_𝑥𝑥+1−𝜈_𝑧𝑥𝜈_𝑥𝑧𝜀_𝑦𝑦+𝐸_𝑥𝜈_𝑧𝑥2𝜈_𝑥𝑧𝛼²_𝑧(Δ𝑇)²,𝐷=𝐸_𝑥1−𝜈_𝑧𝑥𝜈_𝑥𝑧−𝐸_𝑦𝜈²_𝑧𝑦𝜈_𝑥𝑧𝜈_𝑧𝑥+𝜈²_𝑥𝑦+2𝜈_𝑥𝑦𝜈_𝑥𝑧𝜈_𝑧𝑦,𝜀_𝑥𝑥=𝜀_𝑥𝑥−𝜈_𝑧𝑥𝛼_𝑧+𝛼_𝑥Δ𝑇,𝜀_𝑦𝑦=𝜀_𝑦𝑦−𝜈_𝑧𝑦𝛼_𝑧+𝛼_𝑦Δ𝑇. (9a)(9b)(9c)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𝜕𝑊𝜕𝑥_𝑘expl=𝜕𝑊𝜕𝐸_𝑥𝜕𝐸_𝑥𝜕𝑥_𝑘+𝜕𝑊𝜕𝐸_𝑦𝜕𝐸_𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝜈_𝑥𝑦𝜕𝜈_𝑥𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝐺_𝑥𝑦𝜕𝐺_𝑥𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝛼_𝑥𝜕𝛼_𝑥𝜕𝑥_𝑘+𝜕𝑊𝜕𝛼_𝑦𝜕𝛼_𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕()𝜕()Δ𝑇Δ𝑇𝜕𝑥_𝑘,(𝑘=1,2)forplanestress,𝜕𝑊𝜕𝑥_𝑘expl=𝜕𝑊𝜕𝐸_𝑥𝜕𝐸_𝑥𝜕𝑥_𝑘+𝜕𝑊𝜕𝐸_𝑦𝜕𝐸_𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝜈_𝑧𝑦𝜕𝜈_𝑧𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝜈_𝑥𝑧𝜕𝜈_𝑥𝑧𝜕𝑥_𝑘+𝜕𝑊𝜕𝜈_𝑥𝑦𝜕𝜈_𝑥𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝜈_𝑧𝑦𝜕𝜈_𝑧𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝐺_𝑥𝑦𝜕𝐺_𝑥𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝛼_𝑥𝜕𝛼_𝑥𝜕𝑥_𝑘+𝜕𝑊𝜕𝛼_𝑦𝜕𝛼_𝑦𝜕𝑥_𝑘+𝜕𝑊𝜕𝛼_𝑧𝜕𝛼_𝑧𝜕𝑥_𝑘+𝜕𝑊𝜕()𝜕()Δ𝑇Δ𝑇𝜕𝑥_𝑘,(𝑘=1,2)forplanestrain. (10a)(10b)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 (i.e.,Δ𝑇orΔ𝑐) 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 𝐾_I and 𝐾_II. 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:𝐽₁=_𝐴𝜎_𝑖𝑗𝑢_𝑖,1−𝑊𝛿_1𝑗𝑞_,𝑗−𝑑𝐴_𝐴𝑊_,1expl𝑞𝑑𝐴𝑖,𝑗=1,2,𝐽₂=_𝐴𝜎_𝑖𝑗𝑢_𝑖,2−𝑊𝛿_2𝑗𝑞_,𝑗𝑑𝐴−_𝐴𝑊_,2expl−𝑞𝑑𝐴_0𝑅−𝛿𝑊⁺−𝑊⁻𝑞𝑑𝑠−2𝛿𝜋𝐶𝑖,𝑗=1,2, (11a)(11b)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 𝐶=−2𝑎^tip₁₁𝜎⁰₁₁𝛼₁𝜂₂+𝛼₂𝜂₁𝐾_I+𝜂₁+𝜂₂𝐾_II, (12) where 𝜎⁰₁₁ is the constant term in the asymptotic expansion of the stress component 𝜎₁₁ and 𝑎^tip₁₁ is the value of 𝑎₁₁ at the crack tip. 𝑎₁₁ is given by (A.1a) for plane stress and by (A.2a) for plane strain in the appendix. 𝛼_𝑖 and 𝜂_𝑖 are the real and imaginary parts of the two roots of the characteristic equation 𝑎^tip₁₁𝑠⁴−2𝑎^tip₁₆𝑠³+2𝑎^tip₁₂+𝑎^tip₆₆𝑠²−2𝑎^tip₂₆𝑠+𝑎^tip₂₂=0, (13) whose imaginary parts are positive. 𝑎^tip_𝑖𝑗 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 𝜎⁰₁₁, 𝐾_I, and 𝐾_II, which are unknowns. In order to be able to evaluate 𝐽₂, we define a new variable 𝐽₂̂ as follows: 𝐽₂̂=_𝐴𝜎_𝑖𝑗𝑢_𝑖,2−𝑊𝛿_2𝑗𝑞_,𝑗𝑑𝐴−_𝐴𝑊_,2expl−𝑞𝑑𝐴_0𝑅−𝛿𝑊⁺−𝑊⁻𝑞𝑑𝑠. (14) 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 𝐽₂=1𝛿₂/𝜋−𝛿₁/𝜋𝛿₂𝜋𝐽¹₂̂−𝛿₁𝜋𝐽²₂̂. (15) 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 𝐾_I and 𝐾_II [17]:𝐽₁=𝑎^tip₁₁2𝑐₁₁𝐾²_I+𝑐₁₂𝐾_I𝐾_II+𝑐₂₂𝐾²_II,𝐽₂𝑎=−^tip114𝑑₁₁𝐾²_I+𝑑₁₂𝐾_I𝐾_II+𝑑₂₂𝐾²_II, (16a)(16b)where𝑐₁₁=𝛼²₁𝜂₂+𝛼²₂𝜂₁+𝜂²₁𝜂₂+𝜂₁𝜂²₂,𝑐₁₂=2𝛼₁𝜂₂+𝛼₂𝜂₁,𝑐₂₂=𝜂₁+𝜂₂,𝑑₁₁=𝛼₁𝜂₂+𝛼₂𝜂₁,𝑑₁₂=4𝜂²₁𝜂₂+𝜂₁𝜂²₂−𝛼₁𝛼₂𝜂₁+𝜂₂,𝑑₂₂=2𝛼²₁𝜂₂+𝛼²₂𝜂₁+𝛼₁𝛼₂𝜂₁+𝜂₂−𝛼₁𝜂₂−𝛼₂𝜂₁. (17a)(17b)(17c)(17d)(17e)(17f)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 𝐾_I as an arbitrary positive number and 𝐾_II 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]:𝑢₁(𝑟,𝜋)−𝑢₁(𝑟,−𝜋)=22𝑟𝜋𝐷₁₁𝐾_I+𝐷₁₂𝐾_II,𝑢₂(𝑟,𝜋)−𝑢₂(𝑟,−𝜋)=22𝑟𝜋𝐷₂₁𝐾_I+𝐷₂₂𝐾_II, (18a)(18b)where𝐷₁₁𝑠=Im²𝑃₁₁−𝑠₁𝑃₁₂𝑠₁−𝑠₂,𝐷₁₂𝑃=Im¹¹−𝑃₁₂𝑠₁−𝑠₂,𝐷₂₁𝑠=Im²𝑃₂₁−𝑠₁𝑃₂₂𝑠₁−𝑠₂,𝐷₂₂𝑃=Im²¹−𝑃₂₂𝑠₁−𝑠₂,𝑃_1𝑖=𝑎^tip₁₁𝑠²_𝑖+𝑎^tip₁₂−𝑎^tip₁₆𝑠_𝑖(𝑖=1,2),𝑃_2𝑖=𝑎^tip₁₂𝑠_𝑖+𝑎^tip₂₂𝑠_𝑖−𝑎^tip₂₆(𝑖=1,2), (19a)(19b)(19c)(19d)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: 𝐴_𝑟=𝑢₂𝑟_𝑎,𝜋−𝑢₂𝑟_𝑎,−𝜋𝑢₁𝑟_𝑎,𝜋−𝑢₁𝑟_𝑎.,−𝜋 (20) From (18a), (18b), and (20) it then follows that 𝐾_I𝐾_II=𝐴_𝑟𝐷₁₂−𝐷₂₂𝐷₂₁−𝐴_𝑟𝐷₁₁=𝐵. (21) Substituting this result into (16a), 𝐾_I is derived as 𝐾_I=2𝐽₁𝑎^tip₁₁𝑐₁₁+𝑐₁₂/𝐵+𝑐₂₂/𝐵². (22) In this second method, once 𝐽₁ is determined through (11a), 𝐾_I and 𝐾_II 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 𝜕𝑘𝜕𝑥^𝑥𝑥,𝑦𝜕𝑇+𝜕𝜕𝑥𝑘𝜕𝑦^𝑦𝑥,𝑦𝜕𝑇𝜕𝑦=0, (23) 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 𝜕𝐷𝜕𝑥^𝑥𝑥,𝑦𝜕𝑐+𝜕𝜕𝑥𝐷𝜕𝑦^𝑦𝑥,𝑦𝜕𝑐𝜕𝑦=0. (24)𝑐 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 𝐾_I and 𝐾_II 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 𝑎/2000, 𝑎/1000, and 𝑎/8000, 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: 𝑀(𝜉)=𝑀₀+𝑀_ℎ−𝑀₀𝜉ℎ^𝑝, (25) where 𝑀₀ and 𝑀_ℎ are the values of the material property at 𝜉=0 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 𝜉=0 and 𝜉=ℎ and the symbols designating the exponents are provided in Table 1.

Table 1 

Material properties at 𝜉=0 and 𝜉=ℎ and the exponents.

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 𝑎/𝐿. 𝜂𝑥= 𝜂𝑦=2.5, ℎ𝑐/ℎ=0.25, and ℎ/𝐿=0.4.

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 𝐾_I𝑛=𝐾_I𝛼(0)𝐸_𝑥(0)𝑇₀√,𝐾𝜋𝑎_II𝑛=𝐾_II𝛼(0)𝐸_𝑥(0)𝑇₀√.𝜋𝑎 (26) 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. 𝛾_𝑥=𝛾_𝑦=𝛾_𝑥𝑦=2, 𝜔_𝑥𝑦=𝜔_𝑥𝑧=𝜔_𝑧𝑥=𝜔_𝑧𝑦=2, 𝜆_𝑥=𝜆_𝑦=𝜆_𝑧=2.5, 𝜅_𝑥=𝜅_𝑦=𝜅_𝑧=1.5, 𝜂_𝑥=𝜂_𝑦=2.5, 𝜌_𝑥=𝜌_𝑦=1.5, ℎ_𝑐/ℎ=0.25,ℎ/𝐿=0.4, and 𝑎/𝐿=0.1.  

Figure 5 

Deformed shape of the finite element mesh under thermal loading. 𝛾𝑥=𝛾𝑦=𝛾𝑥𝑦=2, 𝜔𝑥𝑦=𝜔𝑥𝑧=𝜔𝑧𝑥=𝜔𝑧𝑦=2, 𝜆𝑥=𝜆𝑦=𝜆𝑧=2.5, 𝜅𝑥=𝜅𝑦=𝜅𝑧=1.5, 𝜂𝑥=𝜂𝑦=2.5, 𝜌𝑥=𝜌𝑦=1.5, ℎ𝑐/ℎ=0.25, ℎ/𝐿=0.4, and 𝑎/𝐿=0.1.

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. 𝛾𝑥=𝛾𝑦=𝛾𝑥𝑦=2, 𝜔𝑥𝑦=𝜔𝑥𝑧=𝜔𝑧𝑥=𝜔𝑧𝑦=2, 𝜆𝑥=𝜆𝑦=𝜆𝑧=2, 𝜅𝑥=𝜅𝑦=𝜅𝑧=2, 𝜂𝑥=𝜂𝑦=2, 𝜌𝑥=𝜌𝑦=2, ℎ𝑐/ℎ=0.25, and ℎ/𝐿=0.4.

(a)

(b)

Figure 7 

Normalized mixed-mode stress intensity factors versus inclination angle 𝜙 and ℎ𝑐/ℎ: (a) mode I SIFs; (b) mode II SIFs. 𝛾𝑥=𝛾𝑦=𝛾𝑥𝑦=2, 𝜔𝑥𝑦=𝜔𝑥𝑧=𝜔𝑧𝑥=𝜔𝑧𝑦=2, 𝜆𝑥=𝜆𝑦=𝜆𝑧=2, 𝜅𝑥=𝜅𝑦=𝜅𝑧=2, 𝜂𝑥=𝜂𝑦=2, 𝜌𝑥=𝜌𝑦=2, ℎ/𝐿=0.4, and 𝑎/𝐿=0.1.

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 𝐾_I𝑛 curves go through maximums at 𝜙 values close to 60°, while 𝐾_II𝑛 attains maximums at 𝜙 values between 0° and 30°. 𝐾_I𝑛 decreases as 𝑎/𝐿 is increased from 1/12 to 1/8. On the other hand, 𝐾_II𝑛 is found to be an increasing function of 𝑎/𝐿 except for relatively small values of 𝜙. We also note that 𝐾_I𝑛 is positive and the crack is open for all 𝜙 considered, whereas 𝐾_II𝑛 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°, 𝐾_I𝑛 is an increasing function of ℎ_𝑐/ℎ; that is, 𝐾_I𝑛 increases as ℎ_𝑐/ℎ is increased from 0.15 to 0.3. When 𝜙 is close to 60°, 𝐾_I𝑛 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 𝐾_II𝑛 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 𝐾_I and 𝐾_II 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𝑎₁₁=𝑆₁₁,𝑎₁₂=𝑆₁₂,𝑎₁₆=2𝑆₁₆,𝑎₂₂=𝑆₂₂,𝑎₂₆=2𝑆₂₆,𝑎₆₆=2𝑆₆₆+𝑆₆₉,𝑔₁=𝛼₁,𝑔₂=𝛼₂,𝑔₁₂=2𝛼₆.ℎ₁=𝛽₁,ℎ₂=𝛽₂,ℎ₁₂=2𝛽₆, (A.1a)(A.1b)(A.1c)(A.1d)and for plane strain𝑎₁₁=𝑆₁₁−𝑆²₁₃𝑆₃₃,𝑎₁₂=𝑆₁₂−𝑆₁₃𝑆₂₃𝑆₃₃,𝑎₁₆=2𝑆₁₆−𝑆₃₆𝑆₁₃𝑆₃₃,𝑎₂₂=𝑆₂₂−𝑆²₂₃𝑆₃₃,𝑎₂₆=2𝑆₂₆−𝑆₂₃𝑆₃₆𝑆₃₃,𝑎₆₆=2𝑆₆₆+𝑆₆₉−2𝑆²₃₆𝑆₃₃,𝑔₁=𝛼₁−𝑆₁₃𝑆₃₃𝛼₃,𝑔₂=𝛼₂−𝑆₂₃𝑆₃₃𝛼₃,𝑔₁₂=2𝛼₆−𝑆₃₆𝑆₃₃𝛼₃,ℎ₁=𝛽₁−𝑆₁₃𝑆₃₃𝛽₃,ℎ₂=𝛽₂,ℎ₁₂=2𝛽₆−𝑆₃₆𝑆₃₃𝛽₃. (A.2a)(A.2b)(A.2c)(A.2d)𝑆_𝑖𝑗,  𝛼_𝑖, and 𝛽_𝑖 in these equalities are derived as follows:𝑆₁₁=1𝐸_𝑥cos²𝜙−sin²𝜙𝜈_𝑥𝑦cos²𝜙+sin²𝜙𝐸_𝑦−𝜈_𝑥𝑦𝐸_𝑥cos²𝜙sin²𝜙+cos²𝜙sin²𝜙𝐺_𝑥𝑦,𝑆₁₂=1𝐸_𝑥cos²𝜙−sin²𝜙𝜈_𝑥𝑦sin²𝜙+sin²𝜙𝐸_𝑦−𝜈_𝑥𝑦𝐸_𝑥cos²𝜙cos²𝜙−cos²𝜙sin²𝜙𝐺_𝑥𝑦,𝑆₁₆=1𝐸_𝑥𝜈_𝑥𝑦sin²𝜙−cos²𝜙+cos𝜙sin𝜙sin²𝜙𝐸_𝑦−𝜈_𝑥𝑦𝐸_𝑥cos²𝜙+1cos𝜙sin𝜙2𝐺_𝑥𝑦cos³𝜙sin𝜙−cos𝜙sin³𝜙,𝑆₂₂=1𝐸_𝑥sin²𝜙−𝜈_𝑥𝑦cos²𝜙sin²𝜙+cos²𝜙𝐸_𝑦−𝜈_𝑥𝑦𝐸_𝑥sin²𝜙cos²𝜙+cos²𝜙sin²𝜙𝐺_𝑥𝑦,𝑆₂₆=1𝐸_𝑥𝜈_𝑥𝑦cos²𝜙−sin²𝜙+cos𝜙sin𝜙cos²𝜙𝐸_𝑦−𝜈_𝑥𝑦𝐸_𝑥sin²𝜙+1cos𝜙sin𝜙2𝐺_𝑥𝑦cos𝜙sin³𝜙−cos³𝜙sin𝜙,𝑆₆₆=1𝐸_𝑥𝜈_𝑥𝑦+1cos²𝜙sin²𝜙+𝜈_𝑥𝑦𝐸_𝑥+1𝐸_𝑦cos²𝜙sin²1𝜙+2𝐺_𝑥𝑦cos⁴𝜙+sin⁴𝜙,𝑆₆₉=1𝐸_𝑥𝜈_𝑥𝑦+1cos²𝜙sin²𝜙+𝜈_𝑥𝑦𝐸_𝑥+1𝐸_𝑦cos²𝜙sin²𝜙−cos²𝜙sin²𝜙𝐺_𝑥𝑦,𝑆₁₃𝜈=−^𝑥𝑧𝐸_𝑥cos²𝜈𝜙−^𝑦𝑧𝐸_𝑦sin²𝜙,𝑆₂₃𝜈=−^𝑥𝑧𝐸_𝑥sin²𝜈𝜙−^𝑦𝑧𝐸_𝑦cos²𝜙,𝑆₃₃=1𝐸_𝑧,𝑆₃₆=𝜈_𝑥𝑧𝐸_𝑥−𝜈_𝑦𝑧𝐸_𝑦cos𝜙sin𝜙,𝛼₁=cos²𝜙𝛼_𝑥+sin²𝜙𝛼_𝑦,𝛼₂=sin²𝜙𝛼_𝑥+cos²𝜙𝛼_𝑦,𝛼₃=𝛼_𝑧,𝛼₆=cos𝜙sin𝜙𝛼_𝑦−𝛼_𝑥,𝛽₁=cos²𝜙𝛽_𝑥+sin²𝜙𝛽_𝑦,𝛽₂=sin²𝜙𝛽_𝑥+cos²𝜙𝛽_𝑦,𝛽₃=𝛽_𝑧,𝛽₆=cos𝜙sin𝜙𝛽_𝑦−𝛽_𝑥. (A.3a)(A.3b)(A.3c)(A.3d)(A.3e)(A.3f)(A.3g)(A.3h)(A.3i)(A.3j)(A.3k)(A.3l)(A.3m)(A.3n)(A.3o)(A.3p)(A.3q)(A.3r)(A.3s)

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
𝐾I:	Mode I stress intensity factor
𝐾II:	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
𝜀11,𝜀22,𝜀12:	Total strains in the crack tip coordinate system
𝜙:	Crack inclination angle
𝜈𝑥𝑦,𝜈𝑥𝑧,𝜈𝑧𝑥,𝜈𝑧𝑦:	Poisson’s ratios
𝜎𝑥𝑥,𝜎𝑦𝑦,𝜎𝑥𝑦:	Stress components in the principal coordinate system
𝜎11,𝜎22,𝜎12:	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.