Computation of Stress Intensity Factor in Functionally Graded Plates under Thermal Shock

*Corr. Author’s Address: Shahrood University of Technology, 7 tir Square, 3619995161 Shahrood, Iran, mbnazari@yahoo.com 622 Computation of Stress Intensity Factor in Functionally Graded Plates under Thermal Shock Nazari, M.B. ‒ Shariati, M. ‒ Eslami, M.R. ‒ Hassani, B. Mohammad Bagher Nazari1,* ‒ Mahmoud Shariati1 ‒ Mohammad Reza Eslami2 ‒ Behrooz Hassani1 1 Shahrood University of Technology, Iran 2 Amir-Kabir University of Technology, Iran


INTRODUCTION
Functionally graded materials (FGMs) are a new type of advanced composites that are introduced for use in high temperature environments.The composition, microstructure and/or crystal structure of the FGMs change gradually, forming a non-homogeneous material with continuously varying thermomechanical properties.In recent years, FGMs have been used widely in other applications.According to the experimental studies of Kawasaki and Watanabe [1], when sudden cooling is applied to ceramic/ metal FGMs, some edge cracks are created on the ceramic surface.Therefore, examining the surface crack problem in FGMs under thermal loading,especially thermal shock, is important in failure analysis of these materials.
Jin and Noda [2] derived the general form of the thermoelastic crack-tip fields in FGMs.They assumed that the material properties are continuous and piecewise differentiable function of spatial position and some of them are not zero at the crack-tip.According to their study, the variation of material properties does not affect the order of singularity of thermoelastic crack-tip fields.Kishimoto et al. [3] showed that in the presence of thermal loading, the path independency of original J-integral is lost.They presented a path-independent form of J-integral included extra term to regard the thermal effect.Analytical approaches including the perturbation method and singular integral equations have been used to consider thermal fracture of FGMs [4] and [5].It is important to know that using analytical approaches is limited to some simple problems or especial conditions.For example, Noda and Guo [5] have studied the edge crack problem in FGMs under thermal shock using the perturbation method.For the sake of simplification, they assumed that the Poisson's ratio is constant.Yildirim [6] and Dag [7] developed an equivalent domain integral to compute the mode-I stress intensity factor (SIF) under steady-state and transient thermal loading in isotropic and orthotropic FGMs, respectively.These analyses were performed by using very fine meshes of regular elements in HEAT2D and FRAC2D software.KC and Kim [8] and [9] used the interaction integral to evaluate the mixed-mode SIFs under steady-state thermal loading.Chen [10] used the interaction integral in conjunction with element-free Galerkin (EFG) method to compute SIFs for an interface crack in orthotropic functionally graded coating under steady-state thermal loading.These results were obtained by using first-order polynomial basis functions, which lead to a fine node arrangement.Also, Chen reported the value of J-integral was

Computation of Stress Intensity Factor in Functionally Graded Plates under Thermal Shock
not completely path-independent and results were unreliable for small integral domain size.
The EFG method provides an efficient and robust framework of analyzing fracture mechanics problems.This method has been implemented for fracture analysis of cracks in FGMs under mechanical loading e.g.[11] or steady-state thermal stresses [10].In this paper, the EFG method is applied in both steady-state and transient thermal fracture of FGMs.The transient thermal loading is imposed in the form of thermal shock.
This paper is organized as follows.Section 1 presents the thermoelastic governing equations.Section 2 provides the EFG discretization form of governing equations.Section 3 explains the use of the equivalent domain integral for thermal fracture of FGMs.Section 4 describes the modal decomposition technique to obtain the transient temperature field.Section 5 presents the obtained numerical results of thermal SIF as well as parametric analyses and the relevant aspects of the results are discussed.Finally, in Section 6 conclusions are drawn.

GOVERNING EQUATIONS
A body occupying a space Ω surrounded by a surface Γ under external actions, body forces and prescribed thermal boundary conditions has been considered.The governing equations for static linear thermoelasticity in the domain Ω are: Also, the heat flux is obtained based on the Fourier law: The constitutive equation is defined as: where, Here, the material properties are the forthorder Hooke tensor C  , isotropic conductivity k, expansion coefficient α, density ρ and specific heat c.The field variables are displacement u, strain tensor ε, stress tensor σ, and thermal strain ε th and the imposed values are heat source Q and body force b.Moreover, I is the identity second-order tensor and ∇ s is the symmetric gradient operator on a vector field.The boundary conditions are as follows: where h is the convection coefficient and n is the outward unit vector which is normal to Γ.

ELEMENT-FREE GALERKIN METHOD IN THERMOELASTICITY
We implement the EFG method to solve governing partial differential equations (PDEs) of 2D thermoelastic problems.This method needs only a set of nodes to construct the discretized model.In EFG, using moving least square (MLS) approximation leads stability in function approximation and applying the Galerkin procedure leads to a stable and well-behaved system of discretized equations.Here, the EFG discretization in the space dimension only is used and the Kantorovitch semi-discretization process is followed.According to the EFG method, the final discrete equations can be obtained as: where the dot ( . ) denotes differentiation with respect to time and: where: and where In the enriched EFG method, the singularity problems due to the presence of a crack is alleviated by enrichment functions.In the intrinsic enrichment, the standard basis (usually polynomials) vector is enriched by including the near-tip asymptotic displacement field [12]: , , , cos , sin , sin sin , cos sin where r and θ are the usual crack-tip polar coordinates.

EQUIVALENT DOMAIN INTEGRAL FOR THERMAL FRACTURE
The J-integral is an energy-based method which is widely used to calculate SIFs.The J-integral originally was derived in the form of a contour integral [13]: where Γ A is an arbitrary contour enclosing the crack-tip and n j is the j th component of the outward unit vector normal to Γ A .For the sake of simplyfying the calculation, it is suitable to convert this contour form into an equivalent domain integral.Defining a smooth weight function q and applying the divergence theorem, the equivalent domain form of the J-integral is derived as [7]: where A is the area inside the contour Γ A .The first integral contains W ,1 , i.e., the partial derivatives of W with respect to x 1 .It should be noted that in FGMs the temperature field and material properties are dependent on spatial coordinates.
In linear elastic fracture mechanics, J-integral is equal to the energy release rate and the relationship between the energy release rate and the mode I SIF is given by: where E E tip tip * = for plane stress and E tip tip ( ) 1 2  −ν for plane strain.E tip and ν tip are Young's modulus and Poisson's ratio, respectively, evaluated at the crack-tip.

TRANSIENT HEAT CONDUCTION PROBLEM
To obtain the temperature field, the firstorder matrix differential equation (6a) should be solved.Among many methods, the modal decomposition technique [14] was chosen.Modal decomposition is an analytical approach to solve systems of ordinary differential equations (ODEs) without the introduction of additional approximations.Based on the modal decomposition procedure, a coupled system of ODEs is turned into uncoupled equations by using eigenvectors.The complete solution of Eq. (6a) can be expressed as a linear combination of all eigenvectors of the system where M is an N×N square matrix whose columns are the eigenvectors.Substituting the above definition into Eq.(6a) and premultiplying it by M T , the uncoupled system of equation is obtained, where The system of Eq. ( 13) contains N uncoupled equations, where and Λ = M T (F th + F γ th ).The initial condition ψ (0) can be obtained from T(0) = M • ψ(0).Depending on the complexity of the right-hand side of Eq. ( 15), it is solved either analytically or numerically.

NUMERICAL RESULTS AND DISSCUSSION
In this section, the calculation of the mode I SIF for an edge crack in functionally graded plate (FGP) under thermal stresses is considered.In addition, a few parametric analyses are performed to study the effect of the gradation of material properties on the stress intensity factor.The distribution of material properties is determined by means of continuum functions e.g., exponential function or micromechanics models e.g., selfconsistent model.Examples are presented here: • An edge cracked plate: exponential gradation.
• FGP with an edge crack: power law gradation.
• Edge crack in an FGP: micromechanics model.The FGP of length W and height H with a crack of length a, as depicted in Fig. 1a, is considered.The thickness (in the x 3 direction) of the plate is assumed to be quite thin for plane stress analysis and large enough for plane strain analysis.The crack is aligned parallel to the direction of material property gradation.Initially, the FGP is at a uniform stress-free temperature T 0 .Temperatures of x 1 = 0 and x 1 = W faces are decreased to constant temperatures T 1 and T 2 , respectively.All other faces, including the crack surfaces, are assumed to be insulated, which results in a dimensional heat conduction problem in the x 1 direction.In all cases, the calculated SIFs will be normalized by dividing to:

An Edge Cracked Plate with Exponentially Gradation
Fig. 1a shows an unconstrained FGP with an edge crack of length a, Fig. 1b presents the complete node arrangement of the FGP which consists of 1695 regular nodes and 40 crack-tip nodes, with a total of 1735.Fig. 1c shows the crack-tip node arrangement.In this case, two different types of functionally graded materials are considered with exponentially varying thermomechanical properties (E, ν, α, k, ρc), in the x 1 direction, e.g.: where the nonhomogeniety parameters are defined e.g., as: Fig. 1

. An FGM plate with an edge crack; a) geometry, b) complete node arrangement, c) crack-tip node arrangement
The values of the nonhomogeneity parameters for the first material are selected arbitrarily (academic materials) as they follow to provide conditions for which the references solutions are available.
For the second case, the ceramic/metal FGM ZrO 2 /Ti-6Al-4V material with properties of Table 1 is assumed.
For the sake of comparison, two different cases of the thermal boundary conditions are considered in the steady-state analysis.In the third case, a transient analysis is also carried out for different temperatures at the left and right sides of the plate.
In order to verify the implementation of DCT and EDI approaches in the framework of EFG method, comparisons of the calculated SIFs and the available reference solutions are first presented.In this case, the temperature of x 1 = 0 and x 1 = W faces are decreased from T 0 to T 1 and T 2 , respectively.Table 2 compares the normalized SIFs with the results provided by Erdogan and Wu [4], KC and Kim [8] and Yildirim [6].The obtained solutions are in good agreement with the references.It is interesting to note that our model is comprised of 1735 nodes, while the 2D mesh discretization in KC and Kim [8] consists of 966 elements and 2937 nodes in the framework of the finite element method.
Since the surface crack is usually created during cooling, the FGP problem subjected to a cooling shock is considered here.To consider the thermal shock, it is assumed that the FGP is initially at a uniform stress-free temperature T 0 and suddenly cooled down to constant temperatures T 1 = 0.2 T 0 and T 2 = 0.5 T 0 at the left and right hand side faces, respectively.The obtained results for the transient temperature distribution in the ZrO 2 /Ti-6Al-4V FGM versus normalized time τ, as defined in Eq. ( 19), is depicted in Fig. 2.
According to these results, the temperature gradient near the plate edges is considerably large at the early times after imposing the thermal shock.
Figs. 3 and 4 present normalized SIFs in the ZrO 2 /Ti-6Al-4V plate resulting from the transient temperature field versus the normalized time τ and the normalized crack length a/W for plane strain and plane stress cases, respectively.As shown in these figures, the SIF quickly increases to a peak value that is drastically larger than the steady value and then decreases rapidly to the corresponding steady value for all crack lengths.In addition, the magnitude of SIF decreases as the normalized crack length a/W becomes larger in both transient and steady states that are in agreement with the results that have recently been reported by Noda and Guo [5].As the final point, the magnitude of SIF for the plane strain is larger than plane stress.Noda et al. [14] have derived thermal stresses analytically for a homogeneous isotropic strip under onedimensional transient temperature distribution.These results indicate that the thermal stresses for the plane strain case are equal to those of the plane stress multiplied by a factor of 1/(1-ν).Regarding the fact 0 < ν < 0.5, this factor is greater than one, which implies a larger SIF for the plane strain in comparison to the plane stress problem, which can be noticed from Figs. 2 and 3.

FGP with an Edge Crack with Power Law Gradation
A Ni/TiC plate with the configuration of the first example is considered here and a powerlaw function is assumed to describe the material properties in the x 1 -direction e.g., as follows.
The exponent p is a positive constant used as an adjusting parameter to obtain certain distribution for material properties.As the exponent p can be chosen independently from the comprised materials, this function is significantly flexible and hence widely used in practice for the analysis of the FGMs.In the proportional material properties, the exponent p is assumed the same value for all material properties while it can be selected differently for non-proportional materials.

Fig. 4. Normalized K I in the ZrO 2 /Ti-6Al-4V plate versus normalized time for different crack lengths in plane stress condition
Moreover, here different thermal boundary conditions are imposed on the uncracked face of FGP.To apply a thermal shock, the cracked face is assumed to be quenched to a constant temperature of T 1 = 0 while having the free convection at the other face with a convection coefficient of h=10 W/(m 2 K) and the ambient temperature is assumed T 0 .The transient temperature distribution in the Ni/TiC plate is presented in Fig. 5 for the proportional case with p = 5.The effect of the convection boundary condition at the x 1 = W face on the temperature distribution is more apparent at the steady-state.Figs. 6 and 7 show the transient thermal SIF versus crack lengths for the proportional case with p = 5 and p = 0.2, respectively.As can be seen, the variation of the thermal SIF is completely different for these cases.In the ceramic-riched case (p = 5), at the beginning of the thermal shock the SIF increases to a peak value and declines to its minimum quickly and then increase gradually to a steady-state value.
However, in the metal-riched case (p = 0.2) the SIF increases quickly to a peak value and then decreases rapidly until the crack is closed.The corresponding time of the crack closure increases as the crack length is increased.In this example, the crack closure occurred in steady-state for all crack lengths.The effect of the thermal boundary condition applied on the uncracked face for the linear proportional material i.e., p = 1, is illustrated in the Fig. 8. Here, the h = 0 corresponds to the insulated thermal boundary condition and h = ∞ corresponds to a known temperature boundary condition.According to these results, while the value of the SIF is independent of the type of the thermal boundary condition applied on the uncracked face, the steady-state value is completely dependent.Moreover, a greater value for the steady-state SIF is obtained for the case of constant temperature at both faces.Now, the effect of the material gradation is studied and some parametric analyses are carried out to assess their effect on the SIFs.In all cases, it is assumed that the exponent p gets different values for the special property and Computation of Stress Intensity Factor in Functionally Graded Plates under Thermal Shock p = 0.2 for other material properties.Figs. 9 and 10 present the effect of variation in FGP elastic properties, i.e.Young's modulus and Poisson's ratio, on the SIF for the crack length a/W = 0.3 and the plane strain problem.According to Fig. 9, the magnitude of SIF, especially its peak value, increases significantly as the parameter p E is increased.These results indicate that for all values of p E , the peak and the crack closure time occur roughly simultaneously.This can be explained by the fact that the transient temperature distribution is independent of the variations of the parameter p .We believe that the effect of Young's modulus is responsible for the slight difference between the peak time and the steady time.
The influence of Poisson's ratio gradation on the SIF is shown in Fig. 10.It can be seen that, by a decrease in p ν , i.e. for metal-riched whose greater Poisson's ratio, the magnitude of SIF and the crack closure time increase.
The analytical solutions for thermal stress distribution in an uncracked FGP under onedimensional temperature distribution for the plane strain and plane stress cases are given as [4]: and respectively, where C 1 and C 2 are unknown coefficients determined from the force and moment boundary conditions in the x 2 direction.From Eq. (21), it is observed that the thermal stresses are an increasing function of Young's modulus.

for other material properties
The effects of the gradation of the thermal properties on the SIF during the shock period are shown in Figs.11, 12 and 13.Fig. 11 depicts the normalized SIF versus normalized time for various values of the exponent p for the thermal expansion coefficient, i.e., p α .According to this figure, the peak value of SIF for the case p = 0.2 is significantly greater than others.Also, depending on the p α value, the trend of SIF might be completely different.For example, the crack is closed for the p α = 0.2 and p α = 1 cases, while the SIF for p α = 5 increases gradually to a steady-state value after it peaks and reduces to a local minimum.

Edge Crack in an FGP with Micromechanics Model
Prediction of the effective macroscopic properties is one of the basic issues in composite material theory.For FGMs, as the graded composites, some micromechanics models of composites have been developed.Among many micromechanics models extended for FGMs, selfconsistent method (SCM) was used.Zuiker [16] has pointed out that the SCM provides a simple and initial estimate of effective properties which is a benefit for relate optimal property distributions.Moreover, in this method the properties are determined independently of the phase of the inclusion and the matrix.This is significant for FGMs in which the volume fraction of the constituent phases varies in a wide range.For twophase FGMs, the volume fraction of the ceramic is assumed in the form of a power function, i.e.V c = 1-(x 1 /L) p , in which L is the material gradation length and the exponent p is known as the gradient index.Here x 1 = 0 corresponds to pure matrix phase (ceramic) and x 1 = L to pure inclusion material (metal).For a two-phase composite, the effective materials are determined from [16], We consider an edge crack in an unconstrained FGP of length W and height H = 8 W. To consider the thermal shock, it is assumed that only the cracked face of the FGP is suddenly cooled down to constant temperature T 1 = 0 from the stress-free temperature T 0 .Fig. 14 presents the transient temperature distribution in the FGP.Here, it is assumed that ΔT = T(x 1 ,t) -T 0 .

Fig. 6 .
Fig. 6.Normalized K I in the Ni/TiC plate versus normalized time and different crack lengths in plane strain condition and p = 5

Fig. 7 . 2 Fig. 8 .
Fig. 7. Normalized K I in the Ni/TiC plate versus normalized time for different crack lengths in plane strain condition and p = 0.2

Fig. 9 .
Fig. 9. Normalized K I versus normalized time for different p E ; plane strain with a/W = 0.3, p = 0.2 for other material properties

Fig. 10 .Fig. 11 .
Fig. 10.Normalized K I versus normalized time for different p ν ; plane strain with a/W = 0.3, p = 0.2 for other material properties

Fig. 12
Fig. 12 illustrates the SIF variation in terms of time for various values of conductivity parameter p k .It can be seen that an increase of the parameter p k causes a delay in the occurrence of the peak value of SIF and steady-state.This delay is not surprising since the diffusivity k/ρc is an increasing function of the conductivity and the p k = 0.2 correspond to metal-riched composition

Fig. 12 .Fig. 13 .
Fig. 12. Normalized K I versus normalized time for different p k .Plane strain with a/W = 0.3, p = 0.2 for other material properties

Fig. 14 .
Fig. 14.Transient temperature distribution in the FGP for various normalized times with T 1 =0

Table 2 .
Normalized mode I SIF in FGP under steady-state thermal loading

Table 3 .
Material properties of Ni and TiC