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

Abstract

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.

1 Introduction

Thermal barrier coatings (TBCs), in particular functionally graded coatings (FGCs), are used in different engineering structures, and the most important one is the application in components where protection against elevated temperatures is required. An example of this application is gas turbine engine parts in which TBCs protect metal parts from overheating and melting [1]. Functionally graded materials (FGMs) are special types of composite materials, the properties of which gradually change along a spatial coordinate. This is achieved by creating variations in material compositions and/or structures. The application of FGC as thermal barrier coatings helps to reduce thermal residual stresses and increases the fracture resistance of the coating.

However, under very high temperatures and complex mechanical loadings, complicated diffusion, oxidation and mechanical interaction processes occur in FGCs structures. As a result of oxidative processes, a thermal growing oxide (TGO) layer is formed [1,2,3]. As reported in the literature [4,5,6], the oxide layer can serve as 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^{\circ }$ 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.

2 Description of the problem

2.1 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,\ldots , 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 $\alpha _{2n}= -\alpha _{\mathrm {2}{n}{\mathrm {+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_{k}^{\mathrm {0}}$, $y_{k}^{\mathrm {0}})$ and the inclination angles $\alpha _{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_{k}^{\mathrm {0}}=x_{k}^{\mathrm {0}}+\textit{iy}_{k}^{\mathrm {0}}$ (i is the imaginary unit).

The FGC/TGO/H structure is cooled by $\Delta T$, $\Delta 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.

Fig. 1

Geometry of the problem: a FGC/TGO/H with a system of cracks, b curved weak layer with cracks, and c the crack system with crack numbering

Fig. 2

Geometry of the problem: a coordinate systems connected with cracks and b edge and internal cracks with fracture angles

2.2 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 used:

$$\begin{aligned} \alpha _{t} (y)=\alpha _{t1} \exp (\varepsilon (y+h)), E(y)=E_{1} \exp (\omega (y+h)), -h\le y\le 0. \end{aligned}$$

(1)

Here, $\alpha _{t}$ is the coefficient of thermal expansion, and E is Young’s modulus with non-homogeneity parameters $\varepsilon$ and $\omega$, respectively. $E_{\mathrm {1}}$ and $\alpha _{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 $\varepsilon h$ and $\omega h$ (h is the thickness of the FGC) are obtained from Eq. (1) as follows:

$$\begin{aligned} \varepsilon h= & {} \ln (\alpha _{t2} /\alpha _{t1} ), \quad \alpha _{t2} =\left. {\alpha _{t} (y)} \right| _{y=0}, \quad \alpha _{t1} =\left. {\alpha _{t} (y)} \right| _{y=-h}, \end{aligned}$$

(2)

$$\begin{aligned} \omega h= & {} \ln (E_{2} /E_{1} ), \quad E_{2} =\left. {E(y)} \right| _{y=0} , \quad E_{1} =\left. {E(y)} \right| _{y=-h}. \end{aligned}$$

(3)

Since the exponential law was chosen for the material parameters $\alpha _{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:

$$\begin{aligned} K_{Ic} (y)=K_{Ic1} \exp (\gamma (y+h)), \end{aligned}$$

(4)

where $\gamma$ is the inhomogeneous parameter of the fracture toughness. The non-dimensional value $\gamma h$ is obtained as

$$\begin{aligned} \gamma h=\ln (K_{Ic2} /K_{Ic1} ),K_{Ic2} =K_{Ic} (y)_{\left| {y=0} \right. }, K_{Ic1} =K_{Ic} (y)_{\left| {y=-h} \right. }. \end{aligned}$$

(5)

In the local coordinate system ($x_{n}$, $y_{n})$ connected with the nth crack, the function $K_{Icn}$ is written as

$$\begin{aligned} K_{Icn} (x_{n} )=K_{Ic1} \exp (\gamma (h+y_{n}^{0} -x_{n} \sin \beta _{n} )). \end{aligned}$$

(6)

Similar expressions are also written for $\alpha _{t}$ and E in the local coordinates $(x_{n}, \ y_{n})$.

2.3 Loadings on the crack faces

With changing the temperature, e.g., an FGM/H structure is cooled on $\Delta 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]:

$$\begin{aligned} \sigma _{xx}^{T} (y)=[\alpha _{t} (y)-\alpha _{t1} ]\Delta TE(y),\quad \sigma _{xx}^{e} (y)=[E(y)/E_{1} -1]\sigma _{xx}^{0} , \sigma _{xx}^{0} =p, \end{aligned}$$

(7)

$\alpha _{t\mathrm {1}}$ and $E_{\mathrm {1}}$ are, respectively, the thermal expansion coefficient and Young’s modulus of a homogeneous substrate material and at the interface, $\alpha _{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

$$\begin{aligned} p_{n} =\sigma _{n} -i\tau _{n} =p(1-\exp (2i\beta _{n} ))/2=pf(\beta _{n} )\; (n = 1, 2,\ldots , N). \end{aligned}$$

(8)

In the common case of FGMs, the full load on the nth crack consists of $p_{n}$, $\sigma _{n}^{T}$ and $\sigma _{n}^{e}$ (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:

$$\begin{aligned} p_{n} +\sigma _{n}^{e} +\sigma _{n}^{T}= & {} \exp (\omega (h+y_{n}^{0} -x_{n} \sin \beta _{n} ))\left[ pf(\beta _{n} ) \right. \nonumber \\&\left. +\,Q\,\exp (\varepsilon (h+y_{n}^{0} -x_{n} \sin \beta _{n} ))-Q \right] \; (n = 1, 2,\ldots , N),\nonumber \\&Q=\alpha _{\text{ t1 }} \Delta TE_{1}. \end{aligned}$$

(9)

It is assumed that $p =Q$, or an additional loading parameter p/Q should be considered.

3 Solution of the problem

3.1 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”.

3.2 Stress intensity factors

The stress intensity factors at tips of the nth crack are obtained from the following formulas:

$$\begin{aligned} K_{nI}^{\pm } -iK_{nII}^{\pm } =\mp \lim \limits _{\eta \rightarrow \pm 1} \sqrt{\pi a_{n} } \sqrt{1-\eta ^{2}} \,{g}'_{n} (\eta ) \end{aligned}$$

(10)

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

$$\begin{aligned} K_{In} -iK_{IIn} =p\sqrt{\pi a_{n} } (k_{In} -ik_{IIn} ). \end{aligned}$$

(11)

3.3 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 $\phi$ (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:

$$\begin{aligned} \phi _{n}= & {} 2\arctan \left[ {\left( {K_{In} -\sqrt{K_{In}^{2} +8K_{IIn}^{2} } } \right) / {4K_{IIn} }} \right] , \end{aligned}$$

(12)

$$\begin{aligned} K_{n}^{eq}\equiv & {} \cos ^{3}(\phi _{n} /2)\left( {K_{In} -3K_{IIn} \tan (\phi _{n} /2)} \right) =K_{Ic,\,tip} \text{ or } K_{n}^{eq} =K_{Ic,\,tip}. \end{aligned}$$

(13)

Using a single crack subjected to a load p normal to the crack line as a reference crack with the stress intensity factor

$$\begin{aligned} K^{0}=p\sqrt{\pi a}, \end{aligned}$$

(14)

the corresponding critical load is obtained as

$$\begin{aligned} p_{0} ={K_{Ic1} } / {\sqrt{\pi a} }, \end{aligned}$$

(15)

where $a=\mathop {\max }\nolimits _{n=1,\ldots ,N} a_{n}$.

By substituting (11) into condition (13), the critical loads are obtained as

$$\begin{aligned} p_{cr\,n} =\frac{K_{Ic} (y)}{\sqrt{\pi a_{n} } }\frac{1}{\cos ^{3}(\phi _{n} /2)\left( {k_{In} -3k_{IIn} \tan (\phi _{n} /2)} \right) } \end{aligned}$$

or

$$\begin{aligned} \frac{p_{cr\,n} }{p_{0} }=\frac{\exp (\gamma (h+y_{n}^{0} -x_{n} \sin \beta _{n} ))}{\cos ^{3}(\phi _{n} /2)\left( {k_{In} -3k_{IIn} \tan (\phi _{n} /2)} \right) }\frac{\sqrt{a} }{\sqrt{a_{n} } }, \end{aligned}$$

(16)

where $p_{\mathrm {0}}$ is defined by Eq. (15) and $K_{Ic}$ by Eqs. (4) and (6).

The fracture angle $\phi _{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

$$\begin{aligned} P_\mathrm{cr} =\mathop {\min }\limits _{n} p_{cr\,n} /p_{0}\; (n = 1, 2, \ldots , N). \end{aligned}$$

(17)

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 $\alpha _{\mathrm {2n}}=- \alpha _{\mathrm {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\,\text {mm}$, the distance between the edge cracks is $d = 2, 4, 6\,(\text {mm})$, and $D=0.25\,\text { mm}$; the dimensionless parameters are used in the numerical calculations: h/a, d/a and D/a, where $a =\hbox { max }a_{n}$ ($n=1, 2,\ldots , 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~\text {mm}$ ($n=1, 2, 3$), internal cracks $a_{n}=0.22~\text {mm}$ ($n=4, 5,\ldots , 9$); $a_{n}$ denotes the half-length of the nth crack. In the calculations the non-dimensional size is used $a_{n}/a$, $a =\hbox { max } a_{n}$. The following designations are also used: $a_\mathrm{edge}$ and $a_\mathrm{int}$ for edge and internal cracks, respectively.

The stress intensity factors (SIFs) are used in dimensionless form: $k_{I,II}=K_{I,II}/K^{\mathrm {0}}$, where $K^{\mathrm {0}} = p(\pi a)^{\mathrm {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.58p(\pi a)^{\mathrm {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.12p(2\pi a)^{\mathrm {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.

Table 1 The coordinates of the centers of cracks

4.1 Homogeneous semi-infinite medium

4.1.1 Stress intensity factors (SIFs)

Figure 3 presents the dimensionless stress intensity factors (SIFs) $k_{I,II}$ for edge cracks depending on $\alpha , 0^{\circ } \le \alpha \le ~60^{\circ }$, and for different distances d/a. The SIFs for internal cracks are shown in Fig. 4 as functions of $\alpha$, and $k_{I,II}^{-}$ refers to the left crack tips (Fig. 4a, c, e), and $k_{I,II}^{+}$—to the right crack tips (Fig. 4b, d, f). Both $k_{I}$ and $k_{II}$ are nonzero; thus, the mixed-mode conditions are realized.

The influence of $\alpha$ (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_{II}$ 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 $\alpha$ and d/a affect $k_{I,II}$; 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 $\alpha =60^{\circ }$ (and also crack 9) have largest values for $k_{I}$ (greater than 1). That is, from these cracks, propagation can start.

Fig. 3

SIFs $k_{I}$ and $k_{II}$ for edge cracks as functions of internal crack angles $\alpha$ and for different distances d between the edge cracks; a for crack 1, b for middle crack 2 and c for crack 3. Homogeneous material

Fig. 4

SIFs $k_{I}$ and $k_{II}$ for internal cracks as functions of angles $\alpha$ and for different distances d between the cracks: a, b for crack 4 at left and right crack tips, respectively, c, d for crack 5 at left and right crack tips, e, f for crack 5 at left and right crack tips. Homogeneous material

4.1.2 Fracture angles

The fracture angles (angles of the crack propagation direction) depending on $\alpha$ 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 $\phi$ 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,\ldots , 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 $\alpha$ has little effect on $\phi$ 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_\mathrm{edge}= 1/a$ and $a_\mathrm{int}= 0.22/a$, the maximum range of change $\phi$ is about $1^{\circ }$ or 7% for $d=2$ (Fig. 5). The influence of the distance d between cracks on $\phi$ is stronger: the smaller d, the larger $\phi$; see Fig. 5 and Table 3. The middle edge crack 2 has $\phi =0^{\circ }$ for all considered parameters of the problem.

Fig. 5

Fracture angles for edge cracks as functions of angles $\alpha$ and for different distances d between the edge cracks: a for crack 1, b for middle crack 2 and c for crack 3. Homogeneous material

Table 2 Fracture angles for edge cracks

Behavior of internal cracks in the system of edge + internal cracks relative to the crack deviation path. Figure 6 shows the fracture angles $\phi$ 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 $\phi$ is highly sensitive to the inclination angle $\alpha$; see Fig. 6. For collinear cracks ($\alpha =0$) with crack sizes $a_\mathrm{edge}=$ 1/a and $a_\mathrm{int}=$ 0.22/a, and distance $d=2$, the $\phi$ value is changed from $0^{\circ }$ to $50^{\circ }$. What is more, the maximum value $\phi$ has inner cracks in the row of these collinear cracks, i.e., cracks 5–8, while the outer cracks 4 and 9 have smaller $\phi$, $0^{\circ }$ for the outer tips and $18^{\circ }$ for the inner ones; see Figs. 6 and 9a. For small internal cracks $a_\mathrm{int}=0.1/a$ (with the size for edge cracks $a_\mathrm{edge}= 1$), the results for outer cracks are opposite: $0^{\circ }$ for the inner tips and $18^{\circ }$ for the outer ones.

The distance d also affects $\phi$. The influence of d on $\phi$ 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, $\alpha =0^{\circ }$, 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^{\circ }$, 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).

Fig. 6

Fracture angles for edge cracks as functions of angles $\alpha$ and for different distances d between the cracks: a, b for crack 4 at left and right crack tips, respectively, c, d for crack 5 at left and right crack tips, e, f for crack 5 at left and right crack tips. Homogeneous material

4.2 FGC/H structures

The inhomogeneity parameters are shown in Table 3 for two FGMs and for temperatures $20\,^{\circ }\hbox {C}$ and $1000\,^{\circ }\hbox {C}$; they are cited as in [3]. For an illustrative example, the parameters of inhomogeneity are taken as follows: $\varepsilon h=-0.5$, $\omega h=-1.5$, $\gamma h=-2.3$. These values correspond to (PSZ/Ni)/Ni with material parameters that increase from the ceramic top to the metal substrate. The geometrical parameters are the same as in the previous section. The fracture angles $\phi$ as functions of $\alpha$ and for different distances d/a are presented in Fig. 7 for edge cracks, and in Fig. 8—for internal cracks. The critical loads are shown in Figs. 10 and 11. Equation (16) is used for dimensionless critical load, and Eq. (12) for fracture angles.

Table 3 The inhomogeneity parameters $\varepsilon h$, $\omega h$ and $\gamma h$ for the thermal expansion coefficient ($\alpha _{t})$, Young’s modulus (E) and fracture toughness ($K_{Ic})$, respectively

4.2.1 Fracture angles

As in the homogeneous case, the inclination angle $\alpha$ does not strongly affect $\phi$ of edge cracks, while the influence of the distance d/a on $\phi$ is noticeable; see Fig. 7 and Table 2. For middle crack 2, the fracture angle $\phi$ is approximately equal to zero. For internal cracks, both $\alpha$ and the distance d/astrongly affect the fracture angles $\phi$ (Fig. 8). However, the maximum $\phi$ values in this case are less than those in the previous homogeneous one, these maxima are equal to $25^{\circ }$ and $50^{\circ }$ for FGC and homogeneous materials, respectively. Figure 9 c, d schematically illustrates the behavior of edge and internal cracks for inclination angles $\alpha =0^{\circ }$ and $30^{\circ }$ ($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.

Fig. 7

Fracture angles $\phi$ for edge cracks as functions of angles $\alpha$ and for different distances d between the cracks: a for crack 1, b for middle crack 2 and c for crack 3. FGC/H structure

Fig. 8

Fracture angles $\phi$ for internal cracks as functions of angles $\alpha$ and for different distances d between the cracks; a, b for crack 4 at left and right crack tips, respectively; c, d for crack 5 at left and right crack tips, e, f for crack 5 at left and right crack tips. FGC/H structure

Fig. 9

Scheme for fracture angles$\phi$ for internal cracks with inclination angles $\alpha =0^{\circ }$ and $30^{\circ }$ and three edge cracks: a, b—for homogeneous medium; c, d—for an FGC/H structure. $a_\mathrm{edge}/a= 1, a_\mathrm{int}/a = 0.22$, $d/a=2$

4.2.2 Critical loads

The critical loads as functions of $\alpha$ 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 $\alpha$ and the distance d/a (Figs. 10 and 11). The weakest crack has the smallest $p_\mathrm{cr}/p_{\mathrm {0}}$. For edge cracks in Fig. 11a, the largest value for $p_\mathrm{cr}/p_{\mathrm {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_\mathrm{cr}/p_{\mathrm {0}}$ are shown for all internal cracks.

It can be seen that the values $p_\mathrm{cr}/p_{\mathrm {0}}$ for crack pairs 4 and 9, 5 and 8, 6 and 7 are equal, as mentioned earlier. The smallest $p_\mathrm{cr}/p_{\mathrm {0}}$ values are for crack 4 at right tip (and accordingly for crack 9 at left tip) for $0^{\circ }<\alpha < 30^{\circ }$, and for crack 4 at left tip (and accordingly for crack 9 at right tip) for $30^{\circ }< \alpha < 50^{\circ }$. Figure 11c shows $p_\mathrm{cr}/p_{\mathrm {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.

Fig. 10

Critical loads as functions of angles $\alpha$ and for different distances d between the edge cracks; a, b for crack 4 at left and right crack tips, respectively, c, d for crack 5 at left and right crack tips, e, f for crack 5 at left and right crack tips; $a_\mathrm{edge}/a = 1$, $a_\mathrm{int}/a = 0.22$. FGC/H structure

Fig. 11

Critical loads as function of angles $\alpha$: a for edge cracks, b for internal cracks, c for edge cracks and internal crack 4; $a_\mathrm{edge}/a = 1$, $a_\mathrm{int}/a = 0.22$, $d/a=2$. FGC/H structure

Behavior of substructure II The 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.

5 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.

Appendices

Appendix A: Singular integral equations

The boundary value problem of elasticity for a system of cracks is reduced to a system of singular integral equations [24]:

$$\begin{aligned}&\int \limits _{-a_{n} }^{a_{n} } \frac{{g}'_{n} (t)\mathrm{d}t}{t-x}+\mathop {\sum }\limits _{\begin{array}{c} k=1 \\ k\ne n \end{array}}^N \int \limits _{-a_{k} }^{a_{k} } [{g}'_{k} (t)R_{nk} (t,x)+\overline{{g}'_{k} (t)} S_{nk} (t,x)]\mathrm{d}t=\pi p_{n} (x), \vert x\vert <a_{n}, n=1, 2, \ldots , N,\nonumber \\&\quad \int \limits _{-a_{n} }^{a_{n} } {{g}'_{n} (t)\mathrm{d}t} =0\qquad \qquad (\text {for internal cracks}),\nonumber \\&{g}'_{n} (x)=\frac{2\mu }{i(\kappa +1)}\frac{\partial }{\partial x}\left( {[u_{n} ]+i[v_{n} ]} \right) . \end{aligned}$$

(A.1)

The number of equations N is equal to the number of cracks. The unknown functions ${g}'_{n} (x)$ contain the shear [$u_{n}$] and normal [$v_{n}$] displacement jumps on the nth crack line, $\mu =E/2(1+\nu )$ is the shear modulus, E is Young’s modulus, $\nu$ is Poisson’s ratio, $\kappa =3-4\nu$ for the plane strain state, and $\kappa =(3-\nu )/(1+\nu )$ for the plane stress state. An overbar $\bar{(...)}$ denotes the complex conjugate. In Eq. (A.1), the functions $p_{n}$ are known functions determined by the load on the crack lines (Eq. (9)). The regular kernels $R_{nk}$(t,x) and $S_{nk}$(t,x) contain the geometry of the problem and are written as

$$\begin{aligned} R_{nk} (t,x)= & {} (1-\delta _{nk} )K_{nk} (t,x)+\frac{e^{i\alpha _{k} }}{2}\left\{ {\frac{1}{X_{n} -\bar{{T}}_{k} }+\frac{e^{-2i\alpha _{n} }}{\bar{{X}}_{n} -T_{k} }} \right. \nonumber \\&+(\bar{{T}}_{k} -T_{k} )\left. {\left[ {\frac{1+e^{-2i\alpha _{n} }}{(\bar{{X}}_{n} -T_{k} )^{2}}-\frac{2e^{-2i\alpha _{n} }(X_{n} -T_{k} )}{(\bar{{X}}_{n} -T_{k} )^{3}}} \right] } \right\} , \end{aligned}$$

(A.2)

$$\begin{aligned} S_{nk} (t,x)= & {} (1-\delta _{nk} )L_{nk} (t,x)+\frac{e^{-i\alpha _{k} }}{2}\left[ {\frac{T_{k} -\bar{{T}}_{k} }{(X_{n} -\bar{{T}}_{k} )^{2}}+\frac{1}{\bar{{X}}_{n} -T_{k} }} \right. -\left. {e^{-2i\alpha _{n} }\frac{X_{n} -T_{k} }{(\bar{{X}}_{n} -T_{k} )^{2}}} \right] , \end{aligned}$$

(A.3)

$$\begin{aligned} T_{k}= & {} te^{i\alpha _{k} }+z_{k}^{0} \quad X_{n} =xe^{i\alpha _{n} }+z_{n}^{0}, \end{aligned}$$

(A.4)

and

$$\begin{aligned} \delta _{nk} =\left\{ {\begin{array}{l} 0\,\,\,\,\text{ for }\,\,\,n\ne k \\ 1\,\,\,\,\,\text{ for }\,\,\,n=k \\ \end{array}} \right. \end{aligned}$$

with kernels $K_{nk}$(t,x) and $L_{nk}$(t,x)

$$\begin{aligned} K_{nk} (t,x)= & {} \frac{e^{i\alpha _{k} }}{2}\left( {\frac{1}{T_{k} -X_{n} }+\frac{e^{-2i\alpha _{n} }}{\overline{T}_{k} -\overline{X}_{n} }} \right) , \end{aligned}$$

(A.5)

$$\begin{aligned} L_{nk} (t,x)= & {} \frac{e^{-i\alpha _{k} }}{2}\left( {\frac{1}{\overline{T}_{k} -\overline{X}_{n} }+\frac{T_{k} -X_{n} }{(\overline{T}_{k} -\overline{X}_{n} )^{2}}e^{-2i\alpha _{n} }} \right) , \end{aligned}$$

(A.6)

which are the same as for the system of cracks in an infinite plane. Additional terms (in addition to $K_{nk}$ and $L_{nk})$ in Eqs. (A.2) and (A.3) take into account the influence of the edge of the half plane. $\alpha _{n}$ is the inclination angle of nth crack to the x-axis; $z^{\mathrm {0}}_{n}$ is the coordinate of the center of crack in the global coordinate system (x,y); see Fig. 2a.

Appendix B: Numerical solution

The singular integral equations (A.1) are solved numerically based on the Gauss–Chebyshev quadrature. Different versions of this method are used for the solution, and the effectiveness of the method has been proven in many studies [25]. In the present work, the version described in [24] is applied.

Equation (A.1) is rewritten in dimensionless form with the non-dimensionless coordinates $\xi =t/a_{k}$ and $\eta =x/a_{n}$, where $2a_{k}$ is the length of the kth crack. The unknown function ${g}'_{n} (\eta )$ presents as

$$\begin{aligned} {g}'_{n} (\eta )=u_{n} (\eta )/\sqrt{1-\eta ^{2}}, \end{aligned}$$

(B.1)

where the function $u_{n} (\eta )$ is a bounded continuous function in the segment $[-1,1]$ and $1/\sqrt{1-\eta ^{2}}$ is the weight function, which is taking into account the square root singularities at the crack tips. For an edge crack, the function ${g}'_{n} (\eta )$ possesses a singularity less than $1/\sqrt{1+\eta }$ at the edge point $\upeta =- 1$; this condition is accounted for as $u_{n} (-1)=0$.

Using Gauss’s quadrature formulae for regular and singular integrals, the singular integral equations are reduced to the following system of $N\hbox {x}M$ (N—number of cracks, M—number of nodes) algebraic equations

$$\begin{aligned}&\frac{1}{M}\sum \limits _{m=1}^M \sum \limits _{k=1}^N \left[ {u_{k} (\xi _{m} )R_{nk} (\xi _{m},\eta _{r} )+} \right. \left. {\overline{u_{k} (\xi _{m} )} S_{nk} (\xi _{m},\eta _{r} )} \right] =\pi p_{n} (\eta _{r} ),\nonumber \\&\sum \limits _{m=1}^M (-1)^{m}u_{n} (\xi _{m} )\tan \frac{2m-1}{4M}\pi =0 \quad (\text {for edge cracks})\ \hbox {or} \sum \limits _{m=1}^M u_{n} (\xi _{m} )=0\quad (\text {for internal cracks}),\nonumber \\&\quad \quad (n=1, 2, \ldots , N; r=1, 2, \ldots , M-1),\nonumber \\&\xi _{m} =\cos \frac{2m-1}{2M}\pi (m~=1, 2,\ldots , M), \eta _{r} =\cos \frac{\pi \,r}{M} (r = 1, 2, \ldots , M-1). \end{aligned}$$

(B.2)

M is the total number of discrete points of the unknown functions $u_{n} (\eta )$ on the segment $[-1,1]$. By applying the conjugate operation to the system (B.2), additional $N\hbox {x}M$ equations are obtained, i.e., $2N\hbox {x}M$ equations should be solved, where N is the number of cracks.

The functions $u_{n} (\eta )$ are calculated by the interpolation formula:

$$\begin{aligned} u_{n} (\eta )=\frac{2}{M}\sum \limits _{m=1}^M u_{n} (\xi _{m} )\sum \limits _{r=0}^{M-1} T_{r} (\xi _{m} )T_{r} (\eta )-\frac{1}{M}\sum \limits _{m=0}^M u_{n} (\xi _{m} ) \end{aligned}$$

(B.3)

$T_{r}$ are Chebyshev polynomials of the first kind.

Taking into account Eqs. (10) and (B.1), the stress intensity factors are obtained as:

$$\begin{aligned} K_{In}^{+} -iK_{IIn}^{+}= & {} -\sqrt{\pi a_{n} } u_{n} (+1) =p_{n} \sqrt{\pi a_{n} } \frac{1}{M}\sum \limits _{m=1}^M (-1)^{m}u_{n} (\xi _{m} )\cot \frac{2m-1}{4M}\pi , \nonumber \\ K_{In}^{-} -iK_{IIn}^{-}= & {} \sqrt{\pi a_{n} } u_{n} (-1) =p_{n} \sqrt{\pi a_{n} } \frac{1}{M}\sum \limits _{m=1}^M (-1)^{M+m}u_{n} (\xi _{m} )\tan \frac{2m-1}{4M}\pi , n=1, 2, \ldots , N. \end{aligned}$$

(B.4)

Here the signs “$+$” and “−“ refer to the right and left crack tips, respectively.

The functions (B.3) written for $\eta =\pm 1$

$$\begin{aligned} u_{n} (1)= & {} \frac{1}{M}\sum \limits _{m=1}^M (-1)^{m+1}u_{n} (\xi _{m} )\cot \frac{2m-1}{4M}\pi \\ u_{n} (-1)= & {} \frac{1}{M}\sum \limits _{m=1}^M (-1)^{M+m}u_{n} (\xi _{m} )\tan \frac{2m-1}{4M}\pi \end{aligned}$$

are used in the calculation of SIFs (Eq. (B.4)).