Strojniški vestnik - Journal of Mechanical Engineering 57(2011)7-8, 622-632 D01:10.5545/sv-jme.2010.138 Paper received: 17.06.2010 Paper accepted: 20.06.2011 Computation of Stress Intensity Factor in Functionally Graded Plates under Thermal Shock Mohammad Bagher Nazari1* - Mahmoud Shariati1 - Mohammad Reza Eslami2 - Behrooz Hassani1 1 Shahrood University of Technology, Iran 2 Amir-Kabir University of Technology, Iran This paper addresses the implementation of the element-free Galerkin method, which is enriched intrinsically for fracture analysis offunctionally graded materials under mode I steady-state and transient thermal loading. The stress intensity factors are evaluated by means of both equivalent domain integral and displacement correlation technique. Continuum functions and the micromechanical model are used to describe the distribution of material properties. For thermal shock analysis, the modal decomposition method, which is a semi-discretization approach, is implemented to obtain the transient temperature field. Also, few parametric analyses are performed to study the effect of material gradation on the stress intensity factors. The results imply that the magnitude of the stress intensity factor reaches its peak a short while after the thermal shock, indicating its significant role in the fracture failure. ©2011 Journal of Mechanical Engineering. All rights reserved. Keywords: functionally graded materials, element-free Galerkin method, equivalent domain integral, displacement correlation technique, thermal stresses 0 INTRODUCTION Functionally graded materials (FGMs) are a new type of advanced composites that are introduced for use in high temperature environments. The composition, micro structure 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 ./-integral is lost. They presented a path-independent form of ./-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 . -integral was *Corr. Author's Address: Shahrood University of Technology, 7 tir Square, 3619995161 Shahrood, Iran, mbnazari@yahoo.com 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. 1 GOVERNING EQUATIONS A body occupying a space Q surrounded by a surface r under external actions, body forces and prescribed thermal boundary conditions has been considered. The governing equations for static linear thermoelasticity in the domain Q are: V-CT + b = 0, V7 ^ dT -Vq + Q = pc —. dt (1a) (1b) Also, the heat flux is obtained based on the Fourier law: q = -kIVT. (2) The constitutive equation is defined as: a= C :(e-e'h ), where, e=Vi u, £,h = a (T - T0)I. (3) (4a) (4b) Here, the material properties are the forth-order Hooke tensor C, isotropic conductivity k, expansion coefficient a, density p and specific heat c. The field variables are displacement u, strain tensor s, stress tensor a, and thermal strain sth and the imposed values are heat source Q and body force b. Moreover, I is the identity second-order tensor and Vs is the symmetric gradient operator on a vector field. The boundary conditions are as follows: T = T on rT, (5a) kl VT • n = q on rq, (5b) klVT ■ n + h(T - Tm) = q on Tc, (5c) u = u on ra, (5d) a ■ n = t on r (5e) where h is the convection coefficient and n is the outward unit vector which is normal to r. 2 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: C thT + (th + Kf ) T = Fth + Ff, U = F + FY, (k + K (6a) (6b) where the dot ( ) denotes differentiation with respect to time and: cth _ Cij _ ■■ [ pcfafydQ, (7a) J Q K j = J kBthT B'fdQ. + J hfdjdr, (7b) f? = f Q^dn+ f qfrdr+ f hojdr, (7c) » o j r J f. f? =y\ e^dr J Tr (7d) where: and B? = dfa/ dx dfa/ dx2 Kh = f DBT BfdQ, J JQ J KY =r\ vTt SVjdT, F, = £ bpfdfi + j^ tyj FY =r\T su + K th*y = Mt (Fth + F'h), (13) where Kth * = m ■ Kth ■ M, Cth * = mt ■ Cth ■ M. (14a) (14b) The system ofEq. (13) contains ^uncoupled equations, Yi + si¥i = \ ,(i = 1,2,..., N), (15) Cii where st = Kf / Cf and A=M(Fth+F/). The initial condition y (0) can be obtained from r(0) = M • y(0). Depending on the complexity of the right-hand side of Eq. (15), it is solved either analytically or numerically. 5 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., self-consistent 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 x3 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 T0. Temperatures of x1 = 0 and x1 = W faces are decreased to constant temperatures T1 and T2, respectively. All other faces, including the crack surfaces, are assumed to be insulated, which results in a dimensional heat conduction problem in the x1 direction. In all cases, the calculated SIFs will be normalized by dividing to: K0 = E(0)a(0)T0-Jnaj(1 -v(0)). (16) 5.1. 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, v, a, k, pc), in the x1 direction, e.g.: E (x,) = E (0) exp(PEx1), (17) where the nonhomogeniety parameters are defined e.g., as: Pe = ¿1* rE (W) ^ E(0) (18) Fig. 1. An FGMplate 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. E(0) = k(0) = a(0) = pc(0) = 1.0, v(0) = 0.3. For the second case, the ceramic/metal FGM ZrO2/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 x1 = 0 and x1 = W faces are decreased from T0 to T1 and T2, 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 T0 and suddenly cooled down to constant temperatures T1 = 0.2 T0 and T2 = 0.5 T0 at the left and right hand side faces, respectively. The obtained results for the transient temperature distribution in the ZrO2/Ti-6Al-4V FGM versus normalized time t, as defined in Eq. (19), is depicted in Fig. 2. k(0)/p(0)c(0), T = - (19) 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 ZrO2/Ti-6Al-4V plate resulting from the transient temperature field versus the normalized time t 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]. 2 Table 1. Material properties of ZrO2 and Ti-6A\-4V Materials Young's modulus [GPa] Poisson's ratio Coefficient of thermal expansion [10-6 /K] Thermal conductivity [W/(m K)] Mass density [kg/m3] Specific heat [J/(kg K)] Z1O2 151.0 0.33 10.0 2.09 5331 456.7 Ti-6Al-4V 116.7 0.33 9.5 7.5 4420 537.0 Table 2. Normalized mode I SIF in FGP under steady-state thermal loading Material parameters Load Analysis type Normalized SIF Present Erdogan and Wu [4] KC and Kim [8] Yildirim [6] EDI DCT WPE=ln(5) WPa=ln(2) T = 0.5 T0 T2 = 0.5 T0 Plane strain 0.0124 0.0126 0.0125 0.0128 0.0128 Plane stress 0.0090 0.0088 - 0.0090 0.0090 T = 0.05 T0 T2 = 0.05 Tq Plane strain 0.0246 0.0240 0.0245 0.244 - WP^ln(5) WPa"ln(2) WPfc=ln(10) T1 = 0.2 T0 T2 = 0.5 TQ Plane strain 0.0334 0.0343 0.0335 0.0334 0.034 Plane stress 0.0234 0.0239 - 0.0235 0.024 T1 = 0.05 T0 T2 = 0.5 Tq Plane strain 0.0405 0.0411 0.0410 0.0406 - -•-1-0 i-le-4 i-le-2 ■ •i-le-l t—1 f ......... \ 1 / I / / ..... N \ S 1 \ L 1 ---------------- °'20 0.2 0.4 06 0.8 1 x/W Fig. 2. Transient temperature distribution in the FGP (ZrO2/Ti-6Al-4V) for various normalized times with T1/ T0 = 0.2 and T2/ T0 = 0.5 0.4 I) .................................................E........-......... 0 0.2 0.4 0.6 0.8 1 X Fig. 3. NormalizedKI in the ZrO2/Ti-6Al-4V plate versus normalized time and different crack lengths in plane strain condition 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 one-dimensional 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-v). Regarding the fact 0 < v < 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. 5.2. 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 power-law function is assumed to describe the material properties in the xrdirection e.g., as follows. E(x1) = E(0) + (E(W) - E(0))(x1 / W)p. (20) 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. -EDI- a/W=0.1 .........EDI- a/W-0.2 ---------EDI- a/W=0 3 o DCT »V \ \ \\\ .v,.,,a„vfc„A.„v,„„S„.„v„vi.„v„„Av..v„„f 0 0.2 04 0.6 08 1 X Fig. 4. Normalized KI in the ZrO2/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 T1 =0 while having the free convection at the other face with a convection coefficient of h=10 W/(m2K) and the ambient temperature is assumed T0. 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 x1 = 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 Table 3. Material properties of Ni and TiC Materials Young's modulus [GPa] Poisson's ratio Coefficient of thermal expansion [10-6 /K] Thermal conductivity [W/(m K)] Mass density [kg/m3] Specific heat [J/(kg K)] TiC 320 0.195 7.4 25.1 4940 134 Ni 206 0.312 13.3 90.5 8890 439.5 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. Fig. 5. Transient temperature distribution in the FGP (Ni/TiC) for various normalized times with Tj = 0 and h = 10 0.3 0.2 -EDl-a/W-O.I -.......[Til-a/W 0? .......EDI- a/W-0.3 ° DCT- a/W-0.1 DCT- a/W=0.2 ? DCT- a/W=0.3 0 0.5 1 Fig. 6. Normalized KI in the Ni/TiC plate versus normalized time and different crack lengths in plane strain condition and p = 5 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. Fig. 7. Normalized KI in the Ni/TiC plate versus normalized time for different crack lengths in plane strain condition andp = 0.2 y 0.02 .......h=0 ......h=l \ .........h=10 V -h=oo V—....... .........................i................................. Fig. 8. The effect of thermal boundary condition at xj = W on the variation of normalized KI 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 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 pE is increased. These results indicate that for all values of pE, 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 pE. 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 pv, 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 one-dimensional temperature distribution for the plane strain and plane stress cases are given as [4]: