Thermal fracture of functionally graded thermal barrier coatings with pre-existing edge cracks and multiple internal cracks imitating a curved interface

This work is devoted to the problem of thermal fracture of a functionally graded coating on a homogeneous substrate (FGC/H) with an emphasis on the analysis of a special system of cracks that simulates a curved interface. The FGC/H structure contains the pre-existing crack system in the FGC, both edge cracks (which are often seen in FGC/H structures) and internal cracks. The stress intensity factors are calculated. (Generally, both Mode I and Mode II are nonzero.) Then, using the appropriate fracture criterion for mixed-mode fracture conditions, the crack propagation direction (so-called fracture angles) and critical loads, when this propagation is initiated, are determined. The application of fracture criteria requires knowledge of the fracture toughness near the crack tips. Thus, it is assumed that the fracture toughness of an FGC, as well as other material properties, continuously varies through the thickness of the coating. For multiple cracks, it is also important to know the weakest crack that starts to propagate first, and the initial direction of this growth. Therefore, the main attention is paid to the evaluation of the fracture angles for the cracks for different parameters of the FGC/H structure. Both cases of a homogeneous semi-infinite medium with a system of cracks imitating a curved interface and FGC/H structures with identical crack systems are studied.

a source of fracture (cracking). Besides, cracks can occur as a result of initial defects or microcracks during manufacturing. Previous investigations [5][6][7][8][9] have revealed that thermal fracture of FGCs is significantly affected by a complex crack interaction mechanism, e.g., interacting cracks can enhance or suppress the propagation of each other. Therefore, the study of fracture of FGCs is important for a better understanding of the fracture resistance of graded coatings.
The interface or oxide layer between coatings and a substrate is rarely perfectly flat. This is usually a wavy, rough interface [4,10,11]. The problem of the roughness of the interface between coatings and a substrate has received a lot of attention in the literature. Different aspects of this problem were discussed in numerous publications, e.g., see [4,[9][10][11][12]. Experimental studies of thermal barrier coatings [4] have reported that a rough interface leads to a higher interface toughness, and the contribution of friction at the interface between the TBC and bond coat may be the cause of this effect. An analytical theoretical study of a wavy non-damaged interface between a thin film and a substrate is presented in [12]. The film thickness and interface roughness are in the nanometer range. The perturbation technique was used within the first-order approximation, where the ratio between the maximum deviation of the interface from the flat state and a perturbation wavelength is considered as a small parameter. The numerical analysis [12] has shown that the interface stress reduces the stress concentration factor, when the residual interface stresses are neglected. This effect decreases when the size of interface asperities (the width and depth), stiffness of the covering film and curvature radius of cavities increase. In [13], analytical solutions are presented for a periodic set of edge dislocations and point forces interacting with a planar traction-free surface of a semi-infinite elastic solid. The fundamental solutions obtained in [13] can be used for applying the boundary integral equation method to an analysis of defects such as cracks and inhomogeneities, periodically distributed at the nanometer distance from the boundary. In [14], a rate-dependent cohesive zone model was presented with a novel traction-separation law with the capability to simulate heat conduction through the interface. The applications of cohesive elements to an FE simulation to describe the cracking of a coating and the delamination of this cracked coating from the substrate material at high temperatures were demonstrated. Interface crack problems were also theoretically studied in [15] for Mode III interface cracks and in [16] for an interface crack and multiple internal cracks in infinite FGM/H bimaterials under thermo-mechanical loading. Flat interfaces were considered in [15,16]. In [11], an experiment and finite element analysis modeling was used to study the influence of a temperature gradient across an YSZ topcoat on the thermal cyclic lifetime of this TBC. A study of the effect of interface roughness on the strain energy release rate and surface cracking behavior in air plasma sprayed thermal barrier coating system was presented in [10]. The extended FEM and periodic boundary conditions were used. Predictions for the stress field and the driving force of multiple surface cracks in the film/substrate system were presented. It was shown that the interface roughness significantly affects the strain energy release rate, the interfacial stress distribution, and the crack formation patterns.
The present work is devoted to the thermal fracture problem for a functionally graded coating on a homogeneous substrate (FGM/H). The model presented in authors' previous papers [17][18][19][20] is used in this study in an extension to a special crack system. The present FGC/H structure contains a pre-existing system of cracks in the FGC. In addition to edge cracks, there are internal/interface cracks imitating a partially wavy interface. The coordinates of the centers of these internal cracks are located on a straight line parallel to the surface, but the inclination angles of these cracks vary symmetrically from 0 • to a certain angle; see Fig. 1. The Mode I and Mode II stress intensity factors are calculated. Then, using the appropriate fracture criterion for mixed-mode fracture conditions, the crack propagation direction and the critical loads at which this propagation begins are determined. For multiple cracks, it is also important to know the weakest crack that starts to propagate first, and the initial direction of this growth, which may indicate the interface roughness in the case of interface cracks. Therefore, the main attention is paid to the evaluation of the fracture angles for the cracks for different parameters of the FGC/H structure. In addition to the FGC/H structures, a homogeneous semi-infinite medium with an identical crack system is also studied as a special case. This paper is organized as follows. Section 2 presents the description of the problem, namely the geometry of the problem, the properties of functionally graded coatings and the corresponding load on the crack faces. Section 3 describes the solution of the problem, which includes the main singular integral equations, determination of stress intensity factors, fracture angles and critical loads. Section 4 presents the results for two cases: for a particular case for a homogeneous medium and for an FGC/TGO/H structure. The parametric results are analyzed with respect to stress intensity factors, fracture angles and critical loads for a special system of cracks, that is, edge cracks and internal-interface cracks, imitating a partly curved weakened interface. The conclusions are presented in Sect. 5.

Geometry of the problem and loading
A structure consisting of a functionally graded layer of thickness h on a semi-infinite homogeneous substrate is shown in Fig. 1a. For curved-rough interfaces, the thickness is determined by the distance from the top surface to the center line of the curved profile. For use in thermal barrier coatings (TBCs), the layer consists of a ceramic on the top layer and mixtures of ceramic and metal gradually transforming into metal on the substrate. TBCs operate at high temperatures at which oxidation processes can occur, resulting in the formation of an oxide layer between the coating and the underlying alloy. That is, the structure under consideration is a functionally graded coating (FGC) (with a material gradient perpendicular to the interface) on a homogeneous substrate with the thermally grown oxide (TGO) layer between them. Special arrangement of internal cracks imitates a weak interface.
In the functionally graded coating, a system of N cracks of length 2a k (k = 1, 2, . . . , N ) is located, which can be edge and/or internal cracks. In the present work, a special system of cracks is considered: three edge cracks normal to the surface and internal-interface cracks with inclination angles α 2n = −α 2n+1 (n = 2, 3, 4) in a curved weak zone (Fig. 1a, b). The numbering of cracks in the system under consideration is shown in Fig. 1c. The global coordinate system (x, y) is set at the top of the FGC surface, and the local coordinates (x k , y k ) are referring to each crack with the x-axis on the crack lines as shown in Fig. 2a. Using these coordinate systems, the positions of the cracks are determined by their midpoint coordinates (x 0 k , y 0 k ) and the inclination angles α k to the x-axis (Fig. 2a). Further, the complex variables method will be used in the formulation of the equations of the problem; therefore, it is convenient to represent the midpoint coordinates of the cracks in complex form as z 0 k = x 0 k + iy 0 k (i is the imaginary unit). The FGC/TGO/H structure is cooled by T , T > 0. (This can be cooling from operating temperatures.) This sudden cooling causes thermal residual stresses in the structure and, in addition, a tensile load p is applied to the structure parallel to the surface.

Model for mechanical and physical properties for an FGC
In functionally graded materials, the composition of the material gradually changes mainly in one direction. In particular, in our problem, the composition varies with the y-coordinate from the ceramic in the upper part of the FGC to the metal in the substrate. Consequently, the thermal and mechanical properties of an FGC also vary continuously with the thickness coordinate y. In this work, an exponential form of these properties is Here, α t is the coefficient of thermal expansion, and E is Young's modulus with non-homogeneity parameters ε and ω, respectively. E 1 and α t1 stand for the thermal and mechanical properties of the homogeneous substrate.
Poisson's ratio is assumed to be constant and is equal to the value of the homogeneous substrate. The values of the dimensionless inhomogeneity parameters εh and ωh (h is the thickness of the FGC) are obtained from Eq. (1) as follows: Since the exponential law was chosen for the material parameters α t and E, it is reasonable to use the exponential law also for the fracture toughness, e.g., see [21]. Thus, the fracture toughness of a functionally graded material can be written as follows: where γ is the inhomogeneous parameter of the fracture toughness. The non-dimensional value γ h is obtained as In the local coordinate system (x n , y n ) connected with the nth crack, the function K I cn is written as Similar expressions are also written for α t and E in the local coordinates (x n , y n ).

Loadings on the crack faces
With changing the temperature, e.g., an FGM/H structure is cooled on T , the residual stresses are arising due to mismatch in the coefficients of thermal expansion [22,23]. In the presented model, the inhomogeneity of FGMs is taken into account through continuously varying residual stresses, and these stresses are the following [21]: α t1 and E 1 are, respectively, the thermal expansion coefficient and Young's modulus of a homogeneous substrate material and at the interface, α t (y) and E(y) are defined by Eq. (1). The method of superposition is used to solve this problem, so that loads at infinity are reduced to the corresponding loads on the crack faces. Thus, the tensile load is reduced to the load p n on the crack surfaces and written in complex form as In the common case of FGMs, the full load on the nth crack consists of p n , σ T n and σ e n (Eq. (7)), where the index "n" denotes that the functions are written in the local coordinate system (x n ,y n ) connected with the nth crack: It is assumed that p = Q, or an additional loading parameter p/Q should be considered.

Singular integral equations
The boundary value problem of elasticity for a system of N cracks is reduced to a system of N singular integral equations [24] with respect to the unknown functions g n (x) containing the shear [u n ] and normal [v n ] displacement jumps on the nth crack line. These equations are written in "Appendix A" (Eq. (A.1)-(A.6)). The singular integral equations (A.1) are solved numerically based on the Gauss-Chebyshev quadrature. Different versions of this method were used for the solution, the effectiveness of which has been proven in many studies [25]. In the present work, the version described in [24] is used, and the solution scheme for this method is given in "Appendix B".

Stress intensity factors
The stress intensity factors at tips of the nth crack are obtained from the following formulas: Here, the signs "+" and "−" refer to the right and left crack tips, respectively. a n is the half-length of the crack. In general, the SIFs can be written as

Stress intensity factors, fracture angles and critical loads
For predicting the crack growth and the determination of the direction of this growing crack, the criterion of maximum circumferential stresses [26] is used. According to this criterion, the crack deflection angle φ (or the so-called fracture angle, as shown in Fig. 2b), the critical stress intensity factor and then the critical stresses are calculated from the following relations: K eq n ≡ cos 3 (φ n /2) (K I n − 3K I I n tan(φ n /2)) = K I c, ti p or K eq n = K I c, ti p .
Using a single crack subjected to a load p normal to the crack line as a reference crack with the stress intensity factor the corresponding critical load is obtained as where a = max n=1,...,N a n .
By substituting (11) into condition (13), the critical loads are obtained as p cr n = K I c (y) √ πa n 1 cos 3 (φ n /2) (k I n − 3k I I n tan(φ n /2)) or p cr n p 0 = exp(γ (h + y 0 n − x n sin β n )) where p 0 is defined by Eq. (15) and K I c by Eqs. (4) and (6). The fracture angle φ n is shown in Fig. 2b. First, the angle of the crack propagation (fracture angle) in Eq. (12) is obtained using the results of the calculated stress intensity factors (Eq. (10)). Next, the local fracture stability is evaluated by Eq. (13). Then, the critical loads are obtained near the crack tips (Eq. (16)). Finally, the weakest crack or crack tip is defined from the condition P cr = min n p cr n / p 0 (n = 1, 2, . . . , N ).
4 Results: stress intensity factors, fracture angles and critical loads for a system of interacting cracks Figure 1a shows the geometry, which is used for numerical example. Three edge cracks are normal to the surface, and the internal cracks have inclination angles α 2n = −α 2n+1 (n = 2, 3, 4) in the curved weak zone. The crack patterns, consisting of edge and internal cracks, have been observed in experiments reported in the literature, e.g., see [5,8]. Table 1 introduces the coordinates of the crack centers. According to this definition of coordinates, when the distance d between the edge cracks changes, the distance between the pairs of cracks 4 and 5, 6 and 7, 8 and 9 also changes. The simplest repeating "crack pattern" in our crack system consists of one edge crack and two internal cracks. There are three such subsystems: pattern I includes cracks 1, 4 and 5, II-cracks 2, 6 and 7, and III-cracks 3, 8 and 9; see Fig. 1c. Pattern II is an inner one with two neighboring patterns I and III. The behavior of pattern II may reflect the behavior of these cracks in a periodic crack system with a repeating pattern II. The results for these subsystems will be done in the following sections.
The results are presented for the following parameters: the width of the coating is h = 2.5 mm, the distance between the edge cracks is d = 2, 4, 6 (mm), and D = 0.25 mm; the dimensionless parameters are used in the numerical calculations: h/a, d/a and D/a, where a = max a n (n = 1, 2, . . . , 9). Results are also obtained for other values for cracks and distances, which are not shown here, but are mentioned in some cases in the discussion. In the following figures, the dimensionless distance d/a is denoted by d. The sizes of cracks are: edge cracks a n = 1 mm (n = 1, 2, 3), internal cracks a n = 0.22 mm (n = 4, 5, . . . , 9); a n denotes the half-length of the nth crack. In the calculations the non-dimensional size is used a n /a, a = max a n . The following designations are also used: a edge and a int for edge and internal cracks, respectively.
The stress intensity factors (SIFs) are used in dimensionless form: k I,I I = K I,I I /K 0 , where K 0 = p(πa) 1/2 is the SIF for a single crack (Eq. (14)). A single edge crack of half-length a, normal to the surface of the layer, has SIF equal to K I = 1.58 p(πa) 1/2 . This definition for SIFs is more convenient for the considered mixed system of internal and edge cracks than the commonly used one for the SIF for an edge crack: K I = 1.12 p(2πa) 1/2 , where 2a is the full length of the edge crack.
Further, two cases are considered: a particular case of a homogeneous semi-infinite medium and a case of a FGC/H structure. In both cases, the values of stress intensity factors and fracture angles are approximately the same for the following cracks: for crack 4 and 9, for crack 5 and 8, for crack 6 and 7. The critical loads are also approximately the same for the listed cracks. Therefore, the results are presented not for all internal cracks, but mostly for cracks 4, 5 and 6. Figure 3 presents the dimensionless stress intensity factors (SIFs) k I,I I for edge cracks depending on α, 0 • ≤ α ≤ 60 • , and for different distances d/a. The SIFs for internal cracks are shown in Fig. 4 as functions of α, Table 1 The coordinates of the centers of cracks Edge cracks 1 2 3

Stress intensity factors (SIFs)
Internal cracks 4 5 6 7 8 9  (Fig. 4a, c, e), and k + I,I I -to the right crack tips (Fig. 4b, d, f). Both k I and k I I are nonzero; thus, the mixed-mode conditions are realized.
The influence of α (the inclination angle of internal cracks) on the SIFs for edge cracks is negligible (Fig. 3). The distance d/a between the cracks affects SIFs: the larger the distance, the greater the k I , and the value k I tends to k I = 1.58, while k I I tends to zero, to the values for a single crack. Besides, k I less than 1.58 for all parameters, that is, the shielding effect is observed. For internal cracks, both α and d/a affect k I,I I ; see Fig. 4. The weakest crack can be defined by the largest value for k I . As seen in Figs. 3 and 4, outer edge cracks 1 and 3 and outer internal crack 4 for α = 60 • (and also crack 9) have largest values for k I (greater than 1). That is, from these cracks, propagation can start.

Fracture angles
The fracture angles (angles of the crack propagation direction) depending on α and for different d are shown in Fig. 5 and in Table 2 for edge cracks 1, 2, 3 and in Fig. 6 for internal cracks 4, 5, 6 (remind that for cracks 4 and 9, 5 and 8, 6 and 7 the values of φ are approximately the same).
Behavior of edge cracks in the system of edge + internal cracks relative to the crack deviation path. If all cracks in the system in Fig. 1 are small cracks (e.g., a n = 0.1/a, n = 1, 2, . . . , 9, these figures are not shown here) the edge cracks have zero or close to zero fracture angles. Thus, if they begin to propagate, the propagation path will be straight. However, as the edge cracks become larger, the fracture angles for outer edge cracks (crack 1 and 3) become larger (Fig. 5). The angle α has little effect on φ for edge cracks. This influence is slightly more pronounced for large cracks, both edge and internal, as would be expected due to the stronger interaction of large cracks. For example, for a edge = 1/a and a int = 0.22/a, the maximum range of change φ is about 1 • or 7% for d = 2 (Fig. 5). The influence of the distance d between cracks on φ is stronger: the smaller d, the larger φ; see Fig. 5 and Table 3. The middle edge crack 2 has φ = 0 • for all considered parameters of the problem.
Behavior of internal cracks in the system of edge + internal cracks relative to the crack deviation path. Figure 6 shows the fracture angles φ for internal cracks at left and right crack tips, and Fig. 9a, b schematically illustrates the behavior of internal cracks for some parameters of the problem. Note that the fracture angle φ is highly sensitive to the inclination angle α; see Fig. 6. For collinear cracks (α = 0) with crack sizes a edge = 1/a and a int = 0.22/a, and distance d = 2, the φ value is changed from 0 • to 50 • . What is more, the maximum value φ has inner cracks in the row of these collinear cracks, i.e., cracks 5-8, while the outer cracks 4 and 9 have smaller φ, 0 • for the outer tips and 18 • for the inner ones; see Figs. 6 and 9a. For small internal cracks    Table 3 The inhomogeneity parameters εh, ωh and γ h for the thermal expansion coefficient (α t ), Young's modulus (E) and fracture toughness (K I c ), respectively 1/a (with the size for edge cracks a edge = 1), the results for outer cracks are opposite: 0 • for the inner tips and 18 • for the outer ones. The distance d also affects φ. The influence of d on φ is greatly pronounced for large cracks (and, accordingly, with stronger interaction); see Fig. 6. For small cracks, a = 0.1/a, the effect of d on the fracture parameters is negligible.
Behavior of cracks in the system of edge + collinear internal cracks, weak cracks and deviation path. For collinear internal cracks, α = 0 • , the SIF k I at outer crack tips (crack 4 left tip and crack 9 right tip) has smaller values in comparison to inner crack tips for internal cracks and also for edge cracks. It means that propagation rarely starts from these tips. More dangerous are outer edge cracks 1 and 3, right tip for crack 4 and left for 5, right tip for crack 8 and left tip for crack 9, so these internal cracks will grow towards each other.
A deviation of cracks from the straight path indicates that cracks can cause some roughness. For example, an initially flat interface tends to become rougher due to this crack deviation path. For the collinear cracks in our system, the magnitude of fracture angles changes from zero to approximately 50 • , as can be seen in Fig. 6 and in the scheme in Fig. 9a. The instability of collinear cracks relative to the deviations of the crack path has been observed in many experimental and theoretical studies, e.g., see [27][28][29]. It was noticed, when two collinear Mode I cracks were growing towards each other, they do not merge tip to tip, but instead repel each other [27]. The origin of this effect has been discussed by several authors. In [28], it was theoretically proved that the straight crack path is unstable, i.e., that tip to tip coalescence will not take place. This was done by considering a periodic array of approximately collinear but slightly curved cracks. In our problem, collinear cracks are loaded parallel to the cracks; in the classical case, the cracks are loaded perpendicular to the cracks (Mode I), and in both these cases, the interaction of cracks causes strong perturbations of the crack path, as can be seen from the change in the fracture angles (Fig. 9).

FGC/H structures
The inhomogeneity parameters are shown in Table 3 for two FGMs and for temperatures 20 • C and 1000 • C; they are cited as in [3]. For an illustrative example, the parameters of inhomogeneity are taken as follows:

Fracture angles
As in the homogeneous case, the inclination angle α does not strongly affect φ of edge cracks, while the influence of the distance d/a on φ is noticeable; see Fig. 7 and Table 2. For middle crack 2, the fracture angle φ is approximately equal to zero. For internal cracks, both α and the distance d/astrongly affect the fracture angles φ (Fig. 8). However, the maximum φ values in this case are less than those in the previous homogeneous one, these maxima are equal to 25 • and 50 • for FGC and homogeneous materials, respectively. Figure 9 c, d schematically illustrates the behavior of edge and internal cracks for inclination angles α = 0 • and 30 • (d/a = 2). It can be seen that the deviation of the edge cracks is approximately the same as for the homogeneous case (Fig. 9a, c). However, the internal cracks have a smaller fracture angle, so the crack path for internal cracks is smoother in FGM than in a homogeneous medium for both collinear and oblique cracks; see Fig. 9.

Critical loads
The critical loads as functions of α are presented in Fig. 10 for internal cracks at left (Fig. 10a, c, e) and right (Fig. 10b, d, f) crack tips for d/a = 2, 4, 6. The other figures refer to the following cases: Fig. 11a for edge cracks, Fig. 11b for internal cracks, and Fig. 11c for edge cracks and internal crack 4, all are for d/a = 2. The results for critical loads take into account the variation of fracture toughness in accordance with the law in Eq. (4).
The critical loads are influenced by both the angle α and the distance d/a (Figs. 10 and 11). The weakest crack has the smallest p cr / p 0 . For edge cracks in Fig. 11a, the largest value for p cr / p 0 is for crack 2 and the smallest one for cracks 1 and 3. Thus, the fracture can start from the edge cracks 1 and 3. In Fig. 11b, the critical loads p cr / p 0 are shown for all internal cracks.
It can be seen that the values p cr / p 0 for crack pairs 4 and 9, 5 and 8, 6 and 7 are equal, as mentioned earlier. The smallest p cr / p 0 values are for crack 4 at right tip (and accordingly for crack 9 at left tip) for 0 • < α < 30 • , and for crack 4 at left tip (and accordingly for crack 9 at right tip) for 30 • < α < 50 • . Figure 11c shows p cr / p 0 for internal crack 4 and for edge cracks; as can be seen, the weakest cracks are edge crack with the lowest values of critical loads.
Behavior of substructure IIThe behavior of the crack pattern II (Fig. 1c) may reflect the behavior of a periodic crack system with pattern II as a period. Comparing the behavior of pattern II with I and III, one can see that this subsystem is more resistant to a possible crack propagation. In both cases, of a homogeneous material and of a FG material, these cracks are less dangerous compared to cracks in patterns I and III (Figs. 4, 10, and 11). A shielding effect is observed for crack 2. At the same time, the crack deviations (fracture angles) of internal cracks (imitating a partly curved interface) are greater than those in patterns I and III for some parameters of the problem; see Figs. 6, 8, and 9. The fracture angles for edge cracks are zero in pattern II.

Summary
Summarizing the obtained results, the following conclusions can be drawn.
The weakest cracks in the considered system are the outer edge cracks (1 and 3) and the outer internal cracks (4 and 9) in a homogeneous medium and the outer edge cracks (1 and 3) in the considered case of FGC. In the homogeneous case, these weak (critical) cracks are determined by the stress intensity factors, and in the case of the FGC-by the critical loads. The fracture criterion of maximal hoop stresses was used, and the structural variation of fracture toughness for the FGM was taken into account. If the material parameters of the FGC tend to homogeneous, critical cracks will be the same as in the homogeneous case.  The fracture angles at the crack tips determine the initial crack paths. These values are strongly influenced by geometric parameters (i.e., crack sizes, distances between the cracks as well as inclination angles). Besides, the inhomogeneity of the FGC also affects the fracture angles. In the considered FGC/H structure, the deviation of cracks from their straight path will be less than that in a homogeneous semi-infinite medium with an identical crack system.
Considering the system of internal collinear cracks, which can imitate a flat interface, it should be noted that the cracks deviate from their straight path, thus these cracks lead to some roughness of the interface. In the case of a homogeneous medium, this roughness is more pronounced than that in the case of FGCs, because of the larger fracture angles of these collinear cracks (Fig. 9a, c).
A system of inclined cracks, imitating a partly curved interface, tends to increase the roughness of the interface in the case of a homogeneous medium (Fig. 9b). However, in the case of the FGC/H system, there is a tendency towards smoothing out the interface roughness (Fig. 9d).
The detailed study of the representative crack patterns, which are observed in thermal barrier coatings, can help find ways to improve the fracture resistance of FGC/H structures.
T k = te iα k + z 0 k X n = xe iα n + z 0 n , (A.4) and δ nk = 0 for n = k 1 for n = k with kernels K nk (t,x) and L nk (t,x)