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

This paper puts forward two different Jk-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 elementmethod and integrated into a general purpose finite element analysis software. Temperature and specificmoisture 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 proposedmethods 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.


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  1 and  2 of the   -integral vector, whereas in the second method the first component  1 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.

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  1 and  2 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.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  1 and  2 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 Δ = 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  1 ,  2 , and  12 .When the results for hygroscopic loading are to be calculated, these quantities need to be replaced by Δ, ℎ 1 , ℎ 2 , and ℎ 12 .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  1 - 2 .  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 s 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.
Although (6) is to be evaluated in the local coordinate system  1 - 2 , 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 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.

𝐽 𝑘 -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 mixedmode SIFs for both thermal and hygroscopic types of loading.
The first method makes use of both components of the  integral,  1 and  2 , whereas in the second method we utilize the first component  1 and the asymptotic displacement fields.

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  2 involves both area and line integrals.The expression of  2 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  0 11 is the constant term in the asymptotic expansion of the stress component  11 and  tip 11 is the value of  11 at the crack tip. 11 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  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  2 given by (11b) contains  0 11 ,  I , and  II , which are unknowns.
In order to be able to evaluate  2 , we define a new variable Ĵ2 as follows: For a given loading condition, Ĵ2 is calculated for two different values of .Denoting these  values by  1 ,  2 , the corresponding Ĵ2 values by Ĵ1 2 , Ĵ2 2 and considering (11b),  2 is expressed as Once  1 is calculated using (11a) and  2 using (15), mixedmode stress intensity factors can be evaluated by the following equalities that relate the components of the   -integral to  I and  II [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  I as an arbitrary positive number and  II as zero.

Method II.
The second method we developed is based on the use of the first component of the   -integral  1 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  1 and  2 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),  I is derived as In this second method, once  1 is determined through (11a),  I and  II can be computed by the use of ( 22) and (21), respectively.

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  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,  1 and  2 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  1 ,  2 , and   , respectively, as /2000, /1000, and /8000, where  is the length of the inclined crack.

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. 0 and  0 values used in hygrothermal loading are also the reference values for the temperature and the specific moisture concentration.In all analyses  0 and  0 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.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  0 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.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.
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.
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.
Figure 6 shows the variations of the normalized mixedmode 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 ℎ  /ℎ.

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   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)

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

Figure 2 :
Figure 2: Integration domains used in the evaluation of the  integral.
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.

𝑘
are the principal thermal conductivities.The governing partial differential equation for the specific moisture concen-

Table 1 :
Material properties at  = 0 and  = ℎ and the exponents.
are derived by coordinate transformation.For the case of plane stress, these elements are given by  11 =  11 ,  12 =  12 ,  16 = 2 16 , (A.1a)  22 =  22 ,  26 = 2 26 ,  66 = 2 ( 66 +  69 ) , (A.1b)  1 =  1 ,  2 =  2 ,  12 = 2 6 .  ,   , and   in these equalities are derived as follows:  − sin 2 ]  ) cos 2   − sin 2 ]  ) sin 2  Elements of the compliance matrix in the principal coordinate system   ,   : Mass diffusivities   ,   ,   : Elastic moduli and shear modulus ℎ: H e i g h t o f t h e o r t h o t r o p i c f u n c t i o n a l l y graded layer ℎ  : H e i g h t o f t h e c r a c k m o u t h   :   -integral vector  I : Mode I stress intensity factor  II : Mode II stress intensity factor   ,   : Th e r m a lc o n d u c t i v i t i e s : Length of the orthotropic functionally graded layer : Radius of the area around the crack tip   : R a d i a l l o c a t i o n a t w h i c h r e l a t i v e displacements of crack faces are calculated 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   ,   ,   : Stresscomponentsintheprincipal coordinate system  11 ,  22 ,  12 : Stresscomponentsinthecracktip coordinate system.