Theoretical modeling and analysis of thermal fracture of semi-infinite functionally graded materials with edge cracks

The present investigation is devoted to a problem of the interaction of two edge cracks inclined arbitrary to the boundary of a non-homogeneous half-plane, which is a functionally graded layer on a homogeneous substrate. The functionally graded properties vary exponentially in thickness direction. One cycle of cooling from sintering temperature is considered. An approach based on integral equations is used and a solution is obtained, then the stress intensity factors are calculated and direction of the initial crack propagation is evaluated by using the maximum circumferential stress criterion. Influence of geometrical and material (inhomogeneity) parameters on the fracture characteristics is investigated. This study can serve as a part of the modeling of the fracture process in FGM coatings under cyclic heating–cooling thermal loading.


Introduction
In different engineering applications, e.g., nuclear energy, aerospace, energy conversions, thermal barrier coating are used to protect metallic or composite components from extremely high temperatures [1][2][3]. Last years for these purposes the so-called functionally graded materials (FGMs) are used. FGMs are composite materials with continuously varying properties in one direction. The application of FGM coatings can reduce bimaterial mismatch at interfaces between the coating and the substrate and prevent delamination and debonding along interfaces. However, cracks can initiate from initial defects or microcracks appear during manufacturing or service. Experimental results [4] showed that when functionally graded plates are subjected to thermal shock, multiple cracks often occur on the ceramic surface during cooling-heating cycles. This fracture process begins from formation of a single crack from initial defects and then a system of edge cracks is formed. Therefore, the study of fracture of FGM coatings is important for a better understanding of the fracture processes in FGM structures and to improve their fracture resistance.
Numerous papers are devoted to different problems of modelling and analysis of fracture processes in FGMs, references can be found in the review papers [2,5]. Different methods are widely used for modeling of FGMs and structures under thermal and mechanical loadings, among them FE methods [6][7][8], the boundary integral methods [9][10][11][12][13][14] and their modifications [15,16]. In spite of many available solutions the problems of interaction of arbitrary located cracks in FGMs are still important.
An approximation method for determining stress intensity factors for a periodic system of edge cracks in an FGM coating in a semi-infinite medium was introduced in [9]. The method is based on singular integral equations. The validity of this approach is discussed in [17] and good accuracy is demonstrated for some gradient (inhomogeneity) parameters of FGMs and crack lengths.
The present work is devoted to the theoretical modeling of fracture of a FGM coating on a semiinfinite homogeneous substrate under thermal and mechanical loading. One cycle of cooling from sintering temperature is considered. It is supposed that two edge cracks arbitrary inclined to the boundary are located in the FGM. The FGM properties are presented by exponential functions. The method of singular integral equations is used and approach similar to the presented in [9] is applied. It is supposed that the inhomogeneity of material is revealed in nonhomogeneous residual stresses on crack surfaces. An example of accounting such residual forces can be found in [9] where a semi-infinite functionally graded material (FGM) with edge cracks is considered. In [18] the influence of an additional loading, which varies with a coordinate along the crack line (and can be considered as a residual stress), on both stress intensity factors Mode I and Mode II was considered in the problem for two parallel cracks under shear loading (pure Mode II) corresponding to the loading in a Compact Shear specimen. In the present investigation the interaction of two edge cracks is considered and this study can serve as a part of the modeling of the fracture process in FGM coatings and further formation of a system of cracks under cyclic thermal loading.

Formulation of the problem
The present investigation is devoted to a problem of the interaction of two edge cracks (with length 2a n , n = 1, 2) inclined arbitrary to the boundary of a nonhomogeneous half-plane (Fig. 1). Cartesian coordinates (x,y) have x-axis along the boundary of the halfplane, local coordinate systems (x n , y n ) are attached to each crack. The crack position is determined by the crack midpoint coordinate z 0 is the imaginary unity) and an inclination angle b n to the boundary, i.e. to the x-axis (Fig. 1).
A functionally graded material (FGM) is located in the region 0 y À h with width h. The Poisson's ratio m is assumed to be a constant because the effect of its variation on the crack-tip stress intensity factors is negligible [19,20]. The remaining thermo-mechanical properties, i.e. the Young's modulus E(y) and the coefficient of thermal expansion a t (y), depend on the y-coordinate only and are modeled by the exponential function [12-14, 19, 20].
For arbitrary located cracks in a half-plane the system of singular integral equations is written as [21,22].
and for two cracks we have two equations. N is the number of cracks, i.e. N = 2 is in our case. The unknown functions are the derivative of displacement jumps on the crack faces; [u n ] and [v n ] are shear and vertical displacement jumps, respectively, on the n-th crack line, l ¼ E=2ð1 þ mÞis the shear modulus, E-the Young's modulus, m-the Poisson's ratio, j ¼ 3 À 4m-for the plane strain state and j ¼ ð3 À mÞ=ð1 þ mÞ-for the plane stress state. The functions p n are determined by the applied load. The regular kernels R nk (t,x) and S nk (t,x) contain geometry of the problem and are cited in the ''Appendix'' in Eqs. (22)- (26).
Equations (1) and (2) are for common case of two arbitrary cracks in a half plane. For edge cracks we should put z 0 n ¼ x 0 n À ia n sin b n . The FGM is cooled from sintering temperature. The FGMs inhomogeneity is accounted via continuously varying residual stresses arising due to mismatch in the coefficients of thermal expansion. The additional stresses p* are the following [9]: These functions will be added to the right side of Eq. (1). a t0 is the thermal expansion coefficient of a homogeneous material (in the region y\ À h) and E 0 is the Young's modulus of this material.
The thermal expansion coefficient of the FGM layer is presented in the exponential form where e is the inhomogeneity parameter of this coefficient. The Young's modulus is with the inhomogeneity parameter d. This exponential model describes the smooth variation of material properties of the FGM in the y-axis direction. For example, if a t1 is decreased with increasing y-coordinate (from y ¼ Àh to y = 0), then the inhomogeneity parameter e is negative, and this case can correspond to a ceramic/metal FGM layer on a metal substrate with gradual transition from a metal at y ¼ Àh to a ceramic at the upper part of the FGM layer (a ceramic t1 \a metal t1 ). The relation between the global coordinates (x,y) and the local coordinate systems (x k ,y k ) can be written in the complex form as follows k is the origin coordinate of the system (x k ,y k ) in the global system. In the local coordinate system (x k ,y k ) connected with each crack the coefficient a t1 Eq. (4) possesses the form a t1 ðx k ; y k Þ ¼ a t0 e eðhþy 0 k Þ e e 1 x k þe 2 y k ; e 1 ¼ e sinðÀb k Þ; e 2 ¼ e cos b k and on the crack lines, where y k = 0, we will have Similar expressions can be written for the Young's modulus.
If we suppose that the Young's moduli of materials in the FGM have approximately same values, then in this case-r e xx ðyÞ ¼ 0 in Eq. (3). It means that the material is elastically homogeneous. Examples of such materials are the following: ceramic/ceramic TiC/SiC, MoSi 2 /Al 2 O 3 and MoSi 2 /SiC, and also ceramic/metal FGMs, e.g., zirconia/nickel and zirconia/steel. For this special case we will investigate the influence of the inhomogeneity parameter e on the fracture characteristics of the materials. For a fully non-homogeneous material we will have two inhomogeneity parameters e and d.
Equations (1) and (2) are rewritten in dimensionless form with the non-dimension coordinates n ¼ t=a and g ¼ x=a, where 2a is a length of the crack (here we suppose that a 1 ¼ a 2 ¼ a). In the considered case of the edge cracks functions g 0 n ðgÞ are bounded in the edge point, i.e. at the point g = -1. At the other tip of the crack, for g = 1, the functions g 0 n ðgÞ as well as the stresses have the square root singularity. The stress intensity factors (SIFs) at the internal tips of the edge cracks are obtained as 3 Solution

Numerical solution
The system of Eq. (1) is solved by the method of mechanical quadrature [21][22][23] which is based on the Chebyshev polynomials. The solution for edge cracks is presented in the following form Here u n ðgÞ are regular functions on the segment [-1,1] and 1= ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À g 2 p is the weight function. Condition that the functions g 0 n ðgÞ are bounded at the edge point g ¼ À1 (or have a singularity less than 1= ffiffiffiffiffiffiffiffiffiffiffi 1 þ g p ) is the following [22] u n ðÀ1Þ ¼ 0: That is the exact singularity at the edge points is not taking into account, but the result obtained with this assumption has shown good accuracy [22].
Using Gauss's quadrature formulae for the regular and the singular integrals the integral equations are reduced to the system of N 9 M (N = 2-number of cracks, M-number of nodes) algebraic equations M is the total number of the discrete points of the unknown functions u n ðgÞ on the interval (-1,1). After solution of the algebraic system (10) and (11) the functions u n ðgÞ are calculated by the interpolation formula: Here T r are the Chebyshev polynomials of the first kind. Setting g ¼ 1 in Eq. (13), it is obtained and for g ¼ À1 we have ðÀ1Þ Mþm u n ðn m Þ tan 2m À 1 4M p This equation and the condition (9) yield Eq. (11).
Applying the conjugate operation to the system (9) additional N 9 M equations are obtained, i.e. 2N 9 M equations should be solved, for two cracks we have 4 9 M equations. Inserting Eq. (14) into the formula (8) and then into Eq. (7) the SIFs are obtained K In À iK IIn ¼ À ffiffiffiffi ffi a n p u n ðþ1Þ ðn ¼ 1; 2Þ: ð15Þ

Validation
To validate this model and verify computational results, we consider some numerical examples and compare our results for SIFs with other available in the literature. The validation of the approximation method similar to the present method with respect to stress intensity factors was discussed in [17]. It was shown that the error depends on the gradient of the profile of Young's modulus of FGMs and crack lengths. For a small crack length the error remains within acceptable limits even for a large gradient (a large inhomogeneity parameter) at the crack tip. However, for a large crack length the gradient (inhomogeneity parameter) should be small at the crack tip for the error to be small. In our modeling the inhomogeneity parameters are used in the range from -1 to 1. Comparison with published in the literature solutions for SIFs can be done for the homogeneous halfplane. The values for SIF for oblique edge cracks can be found in [10,[21][22][23][24][25]. Detailed analysis of these problems was done for one crack in [10] and for periodic edge cracks of unequal cracks in a semiinfinite tensile sheet in [25] for b = 45°. Unfortunately, the results for SIFs for two oblique edge cracks by Nisitani 1977 cited in the handbook of Murakami [24] are not clear for comparison with our results (besides, the original paper by Nisitani 1977 is not available). The comparison have been done with results for SIFs for a single inclined edge crack cited in [10] and with SIFs for periodic edge cracks cited in [25]. Table 1 shows that with increasing the distance d between the cracks the SIF tends to the value for a single edge crack cited in [10]. It should be noted that for the small slope angles b = 15°and 30°and for the distance for d = 10 the values of SIFs are close to the values of SIFs for a single edge crack [10]. With increasing the angle b larger distances d should be taken to achieve the same result, e.g., for b = 90°good result is for the distance d = 100. It means that for small inclination angles (b = 15°and 30°) the interaction between a half-plane boundary and cracks is stronger than the interaction between the cracks. Comparison of the values for SIFs for two interacting edge cracks with the SIFs [25] for a crack in a periodic system of edge cracks inclined on the angle b = 45°shows that the values of SIFs are similar for d = 10 and differ considerably for close located cracks with d = 1, see Table 1.
For all angles b the SIF k I increases with increasing the distance d between the cracks and tends to the value for a single crack, while with decreasing d the SIF k I decreases, i.e. the so-called shielding effect is observed which is known for the parallel cracks under tension normal to the crack lines. Behavior of k II is more complicated than k I . For the angles b = 15°, 30°, 45°and 60°k II increases with increasing the distance between the cracks, and for b = 75°and 90°k II decreases.

Material parameters
For analysis of the fracture development in an FGM/ homogeneous half plane with pre-existing two edge cracks the parameters of the materials should be chosen. Besides, a model for the functional gradation has to be selected. In this study the exponential form for FGMs is used, Eqs. (4)- (6). Then, on the basis of this model and real material combinations of the  Functionally graded materials are used in thermal barrier coating to protect details from high temperatures as well as from wear and corrosion. The materials for protecting from high temperatures should have a low thermal conductivity and at the same time they are desired to have a thermal expansion coefficient close to that of the material for the protected substrate. Consider some actual material combinations ceramic/ceramic and ceramic/metal which can be used in the model. The parameters of these materials, which are available in the literature, are presented in the Table 2 [26]. The Young's modules of these materials are similar; it means that these FGMs are elastically homogeneous, in this case-r e xx ðyÞ ¼ 0 in Eq. (3). For this special case we can investigate the influence of the inhomogeneity parameter e on the fracture characteristics of the material, such as the stress intensity factors at crack tips and the fracture angles which determine the direction of further propagation of cracks.
From Eq. (4) a t1 =a t0 ¼ expðeðy þ hÞÞ and the inhomogeneity parameter e is written as Equation (16) in combination with data in Table 2 shows that the values of e for a thick FGM layer are not large and, accordingly, the actual variation of the residual stresses with coordinate y is not very strong. Figure 2 illustrates the variations of exponential functions expðeðy=h þ 1ÞÞ with the non-dimensional  Fig. 2b. Here e denotes the non-dimensional eh. The value of the exponential function (and, hence, the other values, containing this function, e.g., the thermal expansion coefficient and residual stresses) increases by 35 % for e = 0.3 and decreases by the same value for the negative e equals to e = -1.

Stress intensity factors and fracture angles
To study the effects of inhomogeneity parameters of FGMs on the SIFs at the tips of the edge cracks, as well as on the other fracture characteristic, such as fracture angles, some examples for the problem with respect to loading, material parameters and geometry are considered.
The direction of the initial crack propagation (fracture angle) is evaluated by using the maximum circumferential stress criterion (Cherepanov 1963; Erdogan and Sih 1963; Panasyuk and Berezhnitskij 1964, see for the references [27] ) and is written as For cracks in pure Mode II loading (K I = 0) the fracture angle is calculated as j/ 0 j % 70:5 . For elastically homogeneous materials we can use this criterion without any assumptions. In the case of an elastically non-homogeneous material it is supposed that the material is elastically homogeneous in a vicinity of the crack tips. We will suppose that the cracks have equal lengths a 1 = a 2 = a and the same inclination angles b 1 = b 2 = b to the boundary. These assumptions are made for simplicity of the parametric analysis. Other crack geometries will be studied in other future works.

Two edge cracks in a homogeneous half-plane
For this case we have only tensile loading p parallel to the boundary of the half-plane, which corresponds to the function in the integral Eqs. (1) and (9), here a n ¼ Àb n . The non-dimensional SIFs (k I;II ¼ K I;II =r ffiffiffiffiffi  Fig. 3 and k I \ 1 for all b and d, hence k I is smaller than k I for a solitary crack and the shielding effect is observed. The SIF k II decreases with increasing b from 60°to 120°and changes the sign at b & 105°for the crack 1 and at b & 75°for the crack 2. It means also that at these angles we have pure mode conditions. Figure 4 shows that the fracture angles increases with increasing b and change orientation (sign) at b & 105°for the crack 1 it corresponds to the change of the sign of k II , which characterizes the fracture angles in mixed-mode conditions.
For close located cracks (d = 1, 1.5, 2) the influence of d on k I and k II is small (Fig. 3) but for d = 1, 5, 10 this influence is stronger (Fig. 5). For d = 10 and b & 90°k I [ 1, i.e. greater than k I for a single crack in an infinite plane and close to the value k I for a single edge crack.

Elastically homogeneous material
It was mentioned above that E = E 0 and r e xx ðyÞ ¼ 0 for this case and, hence, Eq. (3) is written as If we suppose that p = Q, i.e. mechanical and thermal loadings are equal, then Otherwise the new parameter p/Q should be considered. In the right side of Eqs. (1) and (10) the function (18) with p determined by (20) should be used. For elastically inhomogeneous materials under tension (without thermal loading, i.e. for r T xx ðyÞ ¼ 0) the load (3) will be determined by the same expressions (19) and (20), where instead of the thermal load Q will be the mechanical load p. In the numerical results the inhomogeneity parameter e is used in the non-dimensional form ea (here a is half size of a crack). Figures 7,8,9,10 show the results of calculation of SIFs and the fracture angles (the maximum circumferential stress criterion was used) for two edge cracks. For this case of thermally non-homogeneous materials we have only one inhomogeneity parameter ea of the thermal expansion coefficient. The SIFs are presented in the non-dimensional form k I;II ¼ K I;II =Q ffiffiffiffiffi 2a p . The calculation were performed for the non-dimensional distances d = 1, 1.5, 2.5 and the non-dimensional h/a = 4. Two values of the inhomogeneity parameter are used ea = -1 and 0.5. The result is obtained on the basis of the solution of the system (10)-(11) for a special case where at the right side of Eq. (1) and accordingly in Eq. (9) is the function (18) with (20). Results for ea = 0.5 for the case, where the thermal expansion coefficient is larger in the upper part of the FGM layer, are presented in Fig. 9. The values of k I are much larger than in the previous case for ea = -1 and for the homogeneous case. The variation of k I with b is small, the changes of k II with b are larger and for both k I and k II are almost linear. The variation of / (Fig. 9) is similar to the previous cases.

Inhomogeneous material, thermal and mechanical loadings
For this case For simplicity of the parametric analysis we assume that the inhomogeneity parameters e and d (of the thermal expansion coefficient and the Young's modulus, correspondingly) are equal and, besides, p = Q as it was in the previous case for thermally inhomogeneous materials. Figures 11 and 12 show the results of calculation of the non-dimensional SIFs (k I;II ¼ K I;II = Q ffiffiffiffiffi 2a p ) for this case of loading. Other parameters are the same as in the previous cases, i.e. d = 1, 1.5, 2.5 and h/a = 4. The results for ea = -0.3, where the thermal expansion coefficient is larger in the upper FGM, are presented in Fig. 11 and they have similar trends as in the thermally inhomogeneous case (Fig. 7, ea = -1) with very small values of k I and k II . Figure 12 shows the values of k I and k II for ea = 0.3 and these results are similar to the previous case shown in Fig. 9 (for thermally inhomogeneous materials for ea = 0.5) with slightly different magnitudes of k I and k II .

Conclusions
Theoretical modeling of thermal fracture of a semiinfinite FGM is presented for the case of one cycle of thermo-mechanical loading and can be used as a part of a study of the fracture process in FGM coatings under cyclic heating-cooling thermal loading.
Influence of geometrical and material (inhomogeneity) parameters on the fracture characteristics is investigated. Strong influence of the inhomogeneity parameter ea of the material on SIFs is observed. If ea is negative, which corresponds to smaller values of the thermal expansion coefficient in the upper part of the FGM layer, the SIFs k I and k II at the edge crack tips are small and have much smaller values in comparison to the homogeneous material. For positive ea, the case where magnitudes of the thermal expansion coefficient have larger values in the upper part of the FGM layer, the SIF k I is much larger than the values for the negative ea and for a homogeneous material. The magnitude of SIF k II is also larger, but not so much as for k I . The fracture angles have similar tendencies for all considered cases and depend mainly on the inclination angle of the cracks to the boundary. The maximum circumferential stress criterion was applied and probably other criteria should be used in order to take into account the inhomogeneity material properties. For example, the strain energy density criterion [28] includes material parameters and can get more results with respect to the fracture angle prediction.