- Email: [email protected]

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings Biao Li a, Xueling Fan a,⇑, Hiroshi Okada b, Tiejun Wang a,⇑ a b

State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, China Department of Mechanical Engineering, Faculty of Science and Technology, Tokyo University of Science, 2641 Noda, Chiba 278-8510, Japan

a r t i c l e

i n f o

Article history: Received 21 February 2017 Received in revised form 17 November 2017 Accepted 24 November 2017 Available online xxxx Keywords: Thermal barrier coatings Surface crack Delamination Finite element analysis J-integral

a b s t r a c t The introduction of dense vertical cracks in ceramic top coat has been recognized as an efficient way to improve the strain tolerance and service durability of thermal barrier coatings (TBCs). However, further growth of these pre-existing cracks may cause the delamination of the coatings to form various failure modes. The motivation of this work is to numerically study the mechanisms that govern the formation of failure modes in dense vertically cracked TBCs. The effects of vertical crack density and material properties on the fracture behaviors are discussed. An energy-based criterion is employed to deal with the surface crack deflection/penetration behavior at a bimaterial interface, and an interaction integral method for discontinuous materials is utilized to evaluate the fracture parameters. The failure maps of TBCs are constructed by examining the continuous variations of crack tip driving forces over the whole fracture process. It is verified that the interfacial fracture resistance could be enhanced by introducing high-density vertical cracks. At the same time, these cracks tend to penetrate into the bond coat that results in debonding at bond coat/substrate interface or at both interfaces. In addition to benefit in alleviating the interfacial mismatch stress, it is found that the high-density vertical cracks give toughening increases to the coating system by releasing more energy from multiple paths. Moreover, the effects of mixed-mode interface toughness and initially residual stress on the formations of failure mode are discussed. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Thermal barrier coatings (TBCs) have been widely used in industrial gas turbines and aircraft engines to protect hot components against high-temperature oxidation and hot corrosion, thereby improving the lifetime of these components and allowing higher operation temperature [1,2]. In general, TBCs are multilayer system which includes an insulating ceramic top coat, a metallic bond coat, and a structural load carrying superalloy substrate [3,4]. Yttria stabilized zirconia (YSZ) is considered as an efficient ceramic top coat material owing to its low thermal conductivity, relatively high coefficient of thermal expansion and favorable environment durability. The metallic bond coat is typically made of MCrAlY (M = Ni, Co or both) to improve the adhesion between top coat and substrate as well as providing oxidation protection to the substrate. Nevertheless, due to the large temperature gradient and thermal expansion mismatch among different layers, high stresses could be developed in the coatings and at the interfaces during thermal cycling, leading to the cracking and spalling of ⇑ Corresponding authors. E-mail addresses: [email protected] (X. Fan), [email protected] (T. Wang). https://doi.org/10.1016/j.engfracmech.2017.11.037 0013-7944/Ó 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

2

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Nomenclature ad, ap adi, api

putative interface crack length and penetration crack length in judging surface crack deflection/penetration putative interface crack length and penetration crack length in judging crack onset

stiffness tensor at crack tip (i, j, k, l = 1, 2) C tip ijkl Ei Young’s modulus (i = 1, 2) F, U, Kd, k nodal force vector, displacement vector, stiffness matrix, and stiffness parameter of discrete fracture element Gd, Gp energy release rates of deflection crack and penetration crack I, J interaction integral, J-integral aux K aux I ; K II ; K I ; K II auxiliary and actual mode I and II stress intensity factors of continuous material crack K, K1, K2 complex stress intensity factor of interface crack, and the mode 1 and mode 2 components q weight function of J-integral r, h local polar coordinates Stip ; Sijkl compliance tensors at the crack tip and at a material point (i, j, k, l = 1, 2) ijkl T 0 ; DT uaux i ; ui a, b th ath TC ; aBC ;

initial temperature, and temperature change referred to initial temperature displacements in the auxiliary and actual fields (i = 1, 2) Dundurs’ parameters ath sub coefficients of thermal expansion of top coat, bond coat and substrate

Cd ; Cp dij di

fracture toughnesses of deflection crack and penetration crack Kronecker delta (i, j = 1, 2) crack flank displacements (i = 1, 2) strains in the auxiliary and actual fields (i, j = 1, 2)

kw ; w

parameter adjusting the contribution of mode mixity, and phase angle

BC/Sub HH TC/BC

bond coat/substrate interface He and Hutchinson’s method top coat/bond coat interface

eaux ij ; eij th etij ; em ij ; eij ; eonset total strain, mechanical strain, thermal strain, and delamination onset strain tip e ; jtip ; l oscillation index, Kolosov constants at crack tip, and characteristic length of interface crack

li ; ltip shear modulus (i = 1, 2), and shear modulus of crack tip material mi ; mtip Poisson’s ratio (i = 1, 2), and Poisson’s ratio of crack tip material rij ; raux stresses in the actual and auxiliary fields (i, j = 1,2) ij

the coatings. In particular, the vertical surface cracking may occur in the cooling stage of thermal cycles, in which there exist large temperature gradient and tensile stress across the coating thickness due to transient thermal shock on the outer surface of the top coat. The tensile stresses may drive the surface crack extending towards the interface. When the crack tip touches the interface, it may further penetrate into the next layer or deflect along the interface. Actually, it is the propagation and coalescence of interface cracks that eventually lead to spallation of the coatings. Therefore, the surface cracking and the induced interfacial delamination are recognized as the dominant failure mechanisms in TBCs [5–7]. To improve the service durability, it is desirable to delay the growth of interface cracks. Fortunately, many experiments have reported the positive effect of surface crack in improving the coating strain tolerance and service durability [8–11]. The thermal shock resistance of TBCs can be significantly improved by introducing a large amount of vertical surface cracks in the top coat [9–12]. These cracks are considered beneficial in reducing the tensile stiffness of the coating, and thus relieving the mismatch stresses. Taylor et al. [9] experimentally found that the thermal cycling life of initially surface cracked coating is superior to initially intact coating. Guo et al. [11] showed that the vertically cracked coatings have an excellent thermal shock resistance. Zhou and Kokini [5,12] observed that the interfacial fracture resistance of the coatings is enhanced with the increase of surface crack density. Moreover, issues of dense vertically cracked TBCs have been extensively studied by numerical methods. Wu et al. [13] analyzed the interfacial stresses in coating-substrate systems considering the periodic surface cracks, and it showed that peak values of the interfacial stresses decrease with the increase of surface crack density. By using fracture mechanics method, Chen et al. [14] found that the energy release rate (ERR) at the interface crack tip decreases with the increase of surface crack density. The similar conclusions have also been obtained in Ref. [15–17]. The failure map for interface crack initiation from the root of surface crack was proposed by Mei et al. [18]. Influence of coating thickness and interfacial adhesion parameters on the interaction between surface cracking and interfacial delamination were discussed by Zhu et al. [19]. Additionally, multiple surface cracking problems incorporating various influence parameters, interface roughness and thermally grown oxide have been systematically investigated by Fan et al. [6,7,20].

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

3

Nevertheless, due to difficulties in dealing with crack initiation, branching and multiple crack propagation in a crack propagation analysis, most of the previous works investigated the failure mechanisms of TBCs by analyzing the driving force for stationary crack. There are few results about the continuous variations of fracture parameters over the whole failure process, which can straightforwardly understand the actual fracture behavior with the variation of applied load. In addition, aforementioned numerical works considered the interfacial delamination only occurring at top coat/bond coat interface (TC/BC interface), in which the bond coat as well as bond coat/substrate interface (BC/Sub interface) were assumed to be strong and tough enough. However, experiments performed by Qian et al. [21], Li et al. [22] and Chen et al. [23] have shown that delamination could arise not only at TC/BC interface but also at the BC/Sub interface (see Fig. 1(a) and (b)). In our experiments, it is observed that the delaminations could simultaneously emanate from both TC/BC interface and BC/Sub interface, as shown in Fig. 1(c), in which there exists vertical surface crack throughout the coating thickness along with interface cracks at both weak interfaces. Generally, it is known that new cracks are able to initiate from weak locations along with the growth of primary crack, leading to the simultaneous existence of multiple cracks in the coatings. For example, when the coatings are loaded by tension, the surface crack in top coat may deflect and propagate along the TC/BC interface. At the same time, a new vertical crack may arise in the bond coat to induce the formation of BC/Sub interface crack. As the applied load increases, these two interface cracks may competitively propagate until either of them reaching the critical length. Therefore, it should be pointed out that apart from spallation at TC/BC interface, interfacial delamination at BC/Sub interface as well as at both interfaces are common fracture modes of TBCs. However, the latter two types of failure modes are neglected in previous investigations. To further improve the durability of TBCs by introducing dense vertical cracks, it is desirable to understand how the surface crack density affects the formation of failure modes, and how the fracture parameters evolve over the whole failure process. Moreover, the coatings are commonly fabricated by thermal spraying at high-temperature. Residual stresses and initial thermal strains may be introduced into the coatings after cooling down to room temperature due to the thermal expansion mismatches among the coating layers. The residual stresses may significantly affect the driving forces at the crack tips so as to

Fig. 1. Failure modes of thermal barrier coatings with a primary surface crack along with interfacial delaminations at (a) BC/Sub interface [23]; (b) TC/BC interface [23]; (c) both interfaces.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

4

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

influence the formation of fracture modes. Thus, the role of residual stresses on the formation of failure modes in the dense vertically cracked TBCs need to be clarified. Unfortunately, little work has addressed above issues. The motivation of this paper is to numerically investigate the mechanisms governing the formation of fracture modes for dense vertically cracked TBCs. We focus on understanding the variations of fracture parameters over the whole fracture process. In Section 2, the problem is formulated. The issue of a surface crack deflection/penetration at a bimaterial interface is presented, and an interaction integral method for discontinuous materials is introduced. In Section 3, the numerical procedures and finite element models are presented. Section 4 discusses the selections of key numerical parameters in the analyses. In Section 5, the variations of fracture parameters are calculated, and effects of surface crack density, elastic mismatch, interfacial properties and residual stresses on the formation of fracture modes are discussed in detail. The conclusions and remarks are summarized in Section 6. 2. Statement of problem 2.1. Problem description In general, failure of TBCs is a three-dimensional fracture problem, while it can be characterized by two-dimensional plane strain model, as shown in Fig. 2. To further simplify the problem, a three-layer unit cell plane strain model is constructed (see Fig. 3). The w is width of the unit cell model and also the spacing of the periodical surface cracks. The hTC and hBC are thicknesses of top coat and bond coat, respectively. Three common failure modes are denoted as failure modes A, B and C. For failure mode A, the vertical surface crack deflects along TC/BC interface and propagates until the top coat spalls off. For failure mode B, the surface crack penetrates into the bond coat, and then deflects and propagates along BC/ Sub interface without debonding at TC/BC interface. For failure mode C, interfacial delaminations appear at both the interfaces simultaneously. Obviously, the analyses of failure modes A and B are relatively simple because they only concern with the surface crack deflection or penetration at TC/BC interface, whereas it is a complex process relating to the new crack initiation, propagation and branching for mode C. 2.2. Crack deflection/penetration at a bimaterial interface For plane problem of a surface crack impinging the bimaterial interface in normal direction with its crack tip located at the interface, namely ap = 0 or ad = 0 in Fig. 3, the traction ahead of crack tip in material #2 (bond coat) can be characterized as [24]

rxx ð0; yÞ ¼ kI ð2pyÞk ;

ð1Þ

where kI is amplitude factor proportional to the applied load, k is real and determined by

cos kp ¼

2ðb aÞ a þ b2 ð1 kÞ2 þ : 1þb 1 b2

ð2Þ

Fig. 2. Schematic representations of failure modes for dense vertically cracked TBCs: (a) failure mode A; (b) failure mode B; (c) failure mode C.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

5

Fig. 3. Three-layer unit cell plane strain model for TBCs: competition between surface crack deflection and penetration at TC/BC interface.

In Eq. (2), Dundurs’ parameters a and b are two pairs of nondimensional elastic mismatch parameters. Under the assumption of plane strain, they are expressed as

a¼

E2 ð1 m21 Þ E1 ð1 m22 Þ ; E2 ð1 m21 Þ þ E1 ð1 m22 Þ

b¼

1 2

l2 ð1 2m1 Þ l1 ð1 2m2 Þ ; l2 ð1 m1 Þ þ l1 ð1 m2 Þ

ð3Þ

ð4Þ

where Ei , mi , and li ¼ Ei =½2ð1 þ mi Þ (i = 1, 2) are the Young’s modulus, the Poisson’s ratio, and shear modulus. The subscripts 1 and 2 are material numbers which represent top coat and bond coat, respectively. It is seen in Eq. (1) that the traction ahead of crack tip relates to the exponent k which is a function of elastic mismatch between materials at each side of the interface. Stress field near the crack tip that touches the interface has the singularity of k, namely r rk , in which r is the distance from the crack tip [25]. For bond coat stiffer than the top coat (a>0), k < 1=2, while for bond coat more compliant than the top coat (a < 0), k > 1=2. Particularly for identical materials, k becomes 1/2, yielding r r1=2 for isotropic homogeneous materials. The stress concentration region dominated by the k singularity is commonly referred to as K-dominant field [26,27]. As the k decreases (increasing a), the intensity of the stress concentration as well as the size of the K-dominant field diminish. The K-dominant field almost vanishes when a>0.5 [26]. The concept of K-dominant field is helpful in understanding the dependence of surface crack deflection/penetration behavior on the presumed putative crack length, as will be discussed in Section 4. Competition between the surface crack deflection and penetration at a bimaterial interface is essentially the choice of lower toughness crack path. It is decided by the nature of stress field as well as the toughness of the components, which are governed by many factors such as geometrical, material and interfacial properties. To solve this problem, He and Hutchinson [24] (hereafter denoted by HH) proposed to compare the ERR at the tip of a small kink crack deflected along the interface, Gd, with the ERR at the tip of a small kink crack penetrated into the next material, Gp. The energy-based crack deflection criterion is expressed as

Cd G d < ; Cp G p

ð5Þ

where Cd denotes the interfacial fracture toughness, and Cp is the fracture toughness of the penetrated material. According to the criterion, crack deflection occurs when the ERR ratio Gd/Gp is greater than the fracture toughness ratio Cd/Cp. Otherwise, the surface crack would penetrate across the interface. The Gd and Gp can be evaluated by using various numerical techniques after the putative crack lengths ad and ap are given. Thereafter, the criterion in Eq. (5) could be utilized to decide the further deflection/penetration behavior of the surface crack. Note that the crack deflection problem in this work is restricted to the doubly-deflected crack case (crack deflection on both sides of the interface), since He et al. [28] and Martinez and Gupta [29] have theoretically demonstrated that doublydeflection is easier than singly-deflection case (crack deflection only on one side of the interface) and should be considered as the dominant crack deflection mode. Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

6

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

2.3. Interaction integral method Understanding the crack propagation behavior in TBCs relies on the accurate calculations of fracture parameters at the crack tips. Among the available methods, J-integral has attracted a great deal of interest owing to its convenience in characterizing the crack tip fields. Based on J-integral, the interaction integral method [30,31] was derived by using the composition of two admissible states, which are known as actual fields and auxiliary fields, to obtain the mode I and mode II stress intensity factors (SIFs) separately. The interaction integral method is an efficient tool to determine the mixed-mode fracture parameters both for homogeneous crack and interface crack types. Over the years, this method has been developed in solving many problems, such as two-dimensional mixed-mode fracture [32,33], interfacial fracture [34,35], nonplanar threedimensional crack propagation [36,37], and nonhomogeneous materials fracture [38,39]. Of particular note is the work by Yu et al. [40–42], who proposed an interaction integral method to consider the situation where discontinuous material properties are enclosed by the integral contour. In their works, domain form of the interaction integral was derived for evaluating the mixed-mode SIFs when the material interfaces are included in the integral domain. Actually, this situation may arise in the fracture problem of this work, in which the integral domain of a vertical crack tip may contain material interface when it extending towards the interface. Therefore, the enhanced interaction integral method in Refs. [40–42] is employed to evaluate the fracture parameters. 2.3.1. Continuous material without interface For homogeneous or nonhomogeneous material without interface on the integral contours, as shown in Fig. 4(a), the standard path-independent J-integral proposed by Rice [43] is expressed as

Z

J ¼ lim

C0 !0

C0

ðxd1j rij ui;1 Þnj dC;

ð6Þ

where rij and ui are stresses and displacements, respectively. C0 is an arbitrary integral contour enclosing the crack tip, dij is Kronecker delta and nj is the unit outward vector normal to the contour C0. The indices i and j denote the components of a variable in local coordinate system originating at the crack tip (i, j = 1, 2). When the thermal strain is included, the strain energy density x is expressed as

1 2

1 2

x ¼ rij emij ¼ rij etij eth ij ;

ð7Þ

t th th th where em is the thermal expansion ij is the mechanical strain, eij is the total strain, eij ¼ a DTdij denotes thermal strain, a coefficient, and DT ¼ T T 0 denotes the temperature change referred to the initial temperature T0. When the crack faces are assumed to be traction-free, the J-integral can be written in the form on a closed contour CA, given as [44]

I

J ¼ lim

C0 !0 C A

ðrij ui;1 xd1j Þmj qdC;

ð8Þ

þ where mj is the outward normal vector normal to the contour CA, with CA ¼ C1 þ C c þ C0 þ Cc (as shown in Fig. 4(a)). q is an arbitrary weight function chosen with a value varying from 1 on C0 to 0 on C1. Considering two independent equilibrium states of a crack body, which are respectively denoted as actual state and auxiliary state. Superposition of these two states results in another equilibrium state, for which the integral generates an interaction term [40]

Fig. 4. Schematic illustrations of the integral domain for homogeneous or nonhomogeneous material: (a) continuous material without interface, where þ integral domain A is enclosed by contour CA , with CA ¼ C1 þ C c þ C0 þ Cc ; (b) discontinuous material containing an interface, where integral domain A = þ A1 + A2 enclosed by the contour CA . The parts A1 and A2 are enclosed by contours CA1 and CA2 , respectively. CA ¼ C11 þ C13 þ C12 þ C c þ C0 þ Cc , CA1 ¼ C11 þ Cþint þ C12 þ Cc þ C0 þ Cþc and CA2 ¼ C13 þ Cint .

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

7

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

I I ¼ lim

C0 !0 C A

aux aux m ðraux ij ui;1 þ rij ui;1 rik eik d1j Þmj qdC;

ð9Þ

where rij , em ij and ui are stress, mechanical strain, and displacement for the actual state; and displacement for the auxiliary field. Generally, the auxiliary fields are defined as follows:

aux raux and uaux are stress, strain, ij , eij i

I II 1=2 1=2 uaux ¼ K aux ; h; ltip ; jtip þ K aux ; h; ltip ; jtip ; I fi r II f i r i

ð10Þ

I 1=2 II 1=2 raux ¼ K aux ; h þ K aux ;h ; ij I g ij r II g ij r

ð11Þ

eaux ¼ Sijkl ðxÞraux ij ij

ð12Þ

ði; j; k; l ¼ 1; 2Þ;

where K aux and K aux are the auxiliary mode I and II SIFs, r and h are local polar coordinates. Kolosov constants jtip ¼ 3 4mtip I II for plane strain and jtip ¼ ð3 mtip Þ=ð1 þ mtip Þ for plane stress. ltip and mtip are shear modulus and Poisson’s ratio at the crack I

II

tip, respectively. Detailed expressions of representative functions f i , f i , g Iij and g IIij are referred to [40,44]. For nonhomogeneous material, the selected auxiliary fields satisfy equilibrium and constitutive relations, but violates compatibility condiaux aux aux ¼ Sijkl ðxÞraux – ðuaux tion, namely raux ij;j ¼ 0 (without body force), eij i;j þ uj;i Þ=2, where Sijkl ðxÞ is the compliance tensor kl , and eij at an arbitrary material point in the integral domain [44]. The equivalent domain integral method is commonly used to calculate the interaction integral, which is obtained by applying divergence theorem to the contour integral of Eq. (9), given by

Z

I¼ A

Z

aux aux m ðraux ij ui;1 þ rij ui;1 rik eik d1j Þq;j dA þ

A

aux aux m ðraux ij ui;1 þ rij ui;1 rik eik d1j Þ;j qdA ¼ I1 þ I 2

ð13Þ

Using the equilibrium condition without body force (rij;j ¼ 0 and raux ij;j ¼ 0), the second integral I2 in Eq. (13) can be written as

I2 ¼ ¼ ¼

Z

aux aux m raux qdA ij ui;1 þ rij ui;1 rik eik d1j ;j

ZA

r

ZA A

aux ij;j ui;1

þr

aux ij ui;j1

þr

aux ij;j ui;1

aux m aux m þ rij uaux i;j1 rij;1 eij rij eij;1 qdA

ð14Þ

aux aux m aux m raux ij ui;j1 þ rij ui;j1 rij;1 eij rij eij;1 qdA

According to the relationship of strain and displacement, the first integrand is derived as

1 2

t aux raux raux ui;j1 þ uj;i1 ¼ raux emij;1 þ eth ij ui;j1 ¼ ij ij eij;1 ¼ rij ij;1 :

ð15Þ

Substituting Eq. (15) into Eq. (14), one obtains

I2 ¼ ¼

Z

m aux th aux aux m aux m raux ij eij;1 þ rij eij;1 þ rij ui;j1 rij;1 eij rij eij;1 qdA

ZA A

ð16Þ

aux m aux th rij uaux i;j1 rij;1 eij þ rij eij;1 qdA

aux In order to facilitate the derivation, a temporary strain tensor eaux0 ¼ Stip ij ijkl rkl is used. According to the definitions of aux-

iliary fields, the

eaux0 satisfies the relations [40] ij

aux aux aux eaux0 ¼ Stip ij ijkl rkl ¼ ui;j þ uj;i

. 2;

ð17Þ

tip where C tip ijkl and Sijkl are stiffness tensor and compliance tensor at the crack tip, respectively. Thus, the first integrand in Eq.

(16) can be written as

1 2

tip aux aux rij uaux rij uaux ¼ rij eaux0 i;j1 ¼ i;j1 þ uj;i1 ij;1 ¼ rij Sijkl rkl;1 :

ð18Þ

m aux The second integrand in Eq. (16) can be replaced by raux ij;1 eij ¼ rij;1 Sijkl ðxÞrkl . Substituting Eq. (16) into Eq. (13), the interaction integral can be expressed as

I ¼ Ih þ Inh þ Ith Z aux aux m ¼ raux ij ui;1 þ rij ui;1 rik eik d1j q;j dA A Z Z aux rij ðStip þ ijkl Sijkl ðxÞÞrkl;1 qdA þ A

A

ð19Þ

raux ij

th ath ;1 DT þ a T ;1 qdA

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

8

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

It is seen in Eq. (19) that the Ih is identical with that for homogeneous material crack. The Inh takes into account contributions of the nonhomogeneity of the material, and Ith accounts for the thermal deformation within the nonhomogeneous material. For pure mechanical loading, the Ith vanish. Moreover, for the case of homogeneous material crack, the Inh and Ith vanish, leading to the interaction integral reduces to the conventional one (I = Ih). It should be noted that the integral contour is chosen arbitrarily and no special assumption has been made. Hence, the obtained interaction integral is path-independent (also domain-independent). 2.3.2. Discontinuous materials with interface For integral domain containing discontinuous material properties (or bimaterial interface), as shown in Fig. 4(b), the domain A is divided into parts A1 and A2 which are enclosed by contours CA1 and CA2 , respectively. It leads to A = A1 + þ þ A2 , CA1 ¼ C11 þ Cþ int þ C12 þ Cc þ C0 þ Cc , CA2 ¼ C13 þ Cint , and CA ¼ C11 þ C13 þ C12 þ Cc þ C0 þ Cc . For the case of discontinuous material, the integral contour CA can also be written as CA ¼ CA1 þ CA2 Cþ int Cint . By aux aux aux m defining P ¼ rij ui;1 þ rij ui;1 rik eik d1j , the integral in Eq. (9) can be expressed as

I

I ¼ lim

C0 !0 C A

Pmj qdC

I

¼ lim

C0 !0 C A1

I

¼ lim

C0 !0 C A1

where Iint ¼

R

Pmj qdC þ lim

I

C0 !0 C A2

Pmj qdC þ lim

I

C0 !0 C A2

ðPmj qÞdC Cþ int

þ

R

Pmj qdC þ

Z Cþ int

ðPmj qÞdC þ

Z C int

ðPmj qÞdC

ð20Þ

Pmj qdC þ Iint

ðPmj qÞdC C int

is the line integral along the interface. Guo et al. [42] proved that Iint only relates

to the variables associated with thermal condition, given as

Z

Iint ¼

Cþ int

th raux ath ii A2 aA1 DT cos c qdC;

ð21Þ

th where c denotes the angle between the x1 axis and the normal of the interface, ath A1 and aA2 are thermal expansion coefficients of domain A1 and A2, respectively. Recalling the domain form of the interaction integral in Eq. (19), Eq. (20) can be further written in the domain form as

I ¼ Ih þ Inh þ Ith þ Iint Z aux aux m ¼ raux ij ui;1 þ rij ui;1 rik eik d1j q;j dA A Z Z aux th rij ðStip raux ath þ ij ;1 DT þ a T ;1 qdA þ Iint ijkl Sijkl ðxÞÞrkl;1 qdA þ A

ð22Þ

A

For linear elastic material, the relationship among J-integral, energy release rate, interaction integral, and SIFs are

1 ðK 2 þ K 2II Þ; E0tip I

ð23Þ

2 K I K aux þ K II K aux ; I II E0tip

ð24Þ

J¼G¼

I¼

where G denotes the energy release rate, KI and KII are the modes I and II SIFs, E0tip ¼ Etip for plane stress and E0tip ¼ Etip =ð1 m2tip Þ for plane strain problems. The Etip and mtip denote the Young’s modulus and the Poisson’s ratio at the crack tip. The auxiliary field can be set as either pure mode I or mode II so that the mixed mode SIFs can be calculated separately. ¼ 1, K aux ¼ 0 and K aux ¼ 0, K aux ¼ 1, and evaluating I equals II and III , respectively, one obtains the actual SIFs K I Selecting K aux I II I II and K II as

KI ¼

E0tip II ; 2

K II ¼

E0tip III : 2

ð25Þ

The expression in Eq. (22) can be reduced to Eq. (19) when the integral domain does not contain any discontinuous material properties. Note that above interaction integral is implemented under the assumption of inverse square root singularity of stress. Calculation of interaction integral is invalid for the case where the crack tip is very close to or just touches the bimaterial interface. Fortunately, the singularity, k, differs only slightly from -1/2 when a > 0 even for crack tip touches the interface, as given in Eq. (1). For the problem in this work, the calculation accuracy is sufficient if we always ensure a short distance between the vertical crack tip and the interface, and skip the case when crack tip just locates on the interface. Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

2.3.3. Interface crack Considering a bimaterial interface crack illustrated in Fig. 5, the complex SIF is K ¼ K 1 þ iK 2 , where i ¼ tractions ahead of the crack tip are given by [45]

9

pﬃﬃﬃﬃﬃﬃﬃ 1. The in-plane

tip

Kr ie ðr22 þ ir12 Þh¼0 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; 2pr

ð26Þ tip

where r and h are local polar coordinates at the crack tip, rie equals zero when elastic mismatch disappears, given by

etip ¼

1 ln 2p

¼ cosðetip ln rÞ þ i sinðetip ln rÞ. etip is the oscillation index and it

!

tip tip jtip 1 =l1 þ 1=l2 ; tip tip j2 =l2 þ 1=ltip 1

in which Kolosov constants

ð27Þ

tip tip tip tip tip for plane stress and jtip jtip m ¼ 3 mm = 1 þ mm m ¼ 3 4mm for the plane strain. lm and mm

are shear modulus and the Poisson’s ratio. m denotes the material number (m = 1, 2). The tractions can also be written in the forms

r22

h i tip Re Kr ie ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; 2pr

r12

h i tip Im Kr ie ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; 2pr

ð28Þ

where Re[] and Im[] denote the real and imaginary parts of the complex number. The crack flank displacements at a distance r behind the tip are tip

d2 þ id1 ¼

8 Kr ie ð1 þ 2ietip Þ coshðpetip Þ E

rﬃﬃﬃﬃﬃﬃﬃ r ; 2p

where dj ¼ uj ðr; h ¼ pÞ uj ðr; h ¼ pÞ with j = 1, 2.

ð29Þ 1 E

¼ 12

1 E01

þ E10 , E0m ¼ Em for plane stress and E0m ¼ Em =ð1 m2m Þ for plane 2

strain problems. Em and mm are the Young’s modulus and the Poisson’s ratio for material m (m = 1, 2). In order to characterize the relative proportion of shear to normal tractions at a fixed distance, l, ahead of the crack tip, the phase angle is commonly used, which is defined as

w ¼ tan1

r12 r22

r¼l

h tip i1 ie Im Kl 1 @ h i A: ¼ tan ietip Re Kl 0

ð30Þ

The choice of characteristic length l is somewhat arbitrary since the phase angle can be easily transformed from one characteristic length to another [45]. Herein, it is taken as top coat thickness l = hTC = 300 lm, as done by many previous works. The interaction integral method in Eq. (19) is still valid to evaluate the SIFs of interface crack [41]. In this case, the auxiliary aux displacement uaux , auxiliary stress raux are, respectively, defined as i ij , and auxiliary strain eij

uaux ¼ i

8 I tip f ðr;h;etip ;j Þ pﬃﬃﬃﬃ > aux aux > < 4li tip coshðpmetip Þ 2rp for K 1 ¼ 1; K 2 ¼ 0 m

pﬃﬃﬃﬃ > f II r;h;etip ;jtip m Þ > r : i ðtip for K aux ¼ 0; K aux ¼1 1 2 tip 2p 4lm coshðpe

ði ¼ 1; 2;

m ¼ 1; 2Þ;

ð31Þ

Þ

Fig. 5. Schematic illustrations of the integral domain of interface crack.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

10

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

aux aux raux ¼ C tip =2 ij ijkl uk;l þ ul;k

eaux ¼ Sijkl ðxÞraux ij kl

ði; j; k; l ¼ 1; 2Þ:

ð32Þ

ði; j; k; l ¼ 1; 2Þ;

ð33Þ I

II

where detailed expressions of representative functions f i , f i are referred to [35,41]. C tip ijkl is the stiffness tensor at the crack tip, which is defined as

C tip ijkl

¼

8 < C tip x1 ¼ 0; x2 ¼ 0þ for x2 > 0 ijkl : C tip ðx1 ¼ 0; x2 ¼ 0 Þ ijkl

for x2 < 0

ð34Þ

By obtaining the derivatives of actual and auxiliary fields, one can obtain the interaction integral from Eq. (19). For interface crack, the relationship among J-integral, interaction integral, and SIFs are given as

1 K 21 þ K 22 ; E cosh2 ðpetip Þ

ð35Þ

2 K 1 K aux þ K 2 K aux 1 2 : E cosh2 ðpetip Þ

ð36Þ

J¼G¼

I¼

Similar to that of homogeneous material crack, the actual SIFs K 1 and K 2 for interface crack could be obtained by letting ¼ 1, K aux ¼ 0 and K aux ¼ 0, K aux ¼ 1, respectively. They are written as K aux 1 2 1 2

E cosh ðpetip Þ II ; 2 2

K1 ¼

E cosh ðpetip Þ III : 2 2

K2 ¼

ð37Þ

It should be pointed out that the interface crack SIFs are denoted by K 1 and K 2 , but not by K I and K II , since tension and shear effects are inseparable in the vicinity of interface crack tips. K 1 and K 2 are not the conventional mode I and mode II SIFs [41]. 2.4. Fracture criterion The energy release rate based fracture criterion is employed to gauge the critical condition for crack growth. A crack extends if

G P CC ;

ð38Þ

where CC is the mixed-mode fracture toughness. Due to the symmetry of the model, the vertical crack in bond coat extends under pure mode I and the fracture occurs when the ERR is greater than the mode I fracture toughness. However, the interfacial fracture is inherently of mixed mode due to the asymmetry of the materials at each side of the interface. Several phenomenological functions have been proposed to characterize the mixed-mode interface toughness [46]. The one parameter mixed-mode interface toughness function is employed, given by

CC ¼ CðwÞ ¼ C1 1 þ tan2 ð1 kw Þw ;

ð39Þ

where C is the mode 1 interfacial fracture toughness, kw is a parameter adjusting the contribution of mode mixity. The inter1

face toughness function CðwÞ is symmetric in phase angle w. The ideally brittle interface with kw ¼ 1 indicates the CC ¼ C1 for all mode combinations. For kw – 1, the toughness CðwÞ rapidly increases as w increases towards 90°. 3. Finite element analysis 3.1. Discrete fracture element According to experimental observations, the crack paths in TBCs may be along the TC/BC interface, the BC/Sub interface or perpendicularly in top coat and bond coat. The prior knowledge of potential crack paths makes it possible to use the discrete fracture elements to simulate the crack propagation without remeshing. By connecting node pair of different material elements, the discrete fracture element acts like a stiff spring when the material or the interface is intact, while it is opened by setting its stiffness to be zero to simulate the crack propagation. This method is also recognized as nodal release approach. As illustrated in Fig. 6, the discrete fracture element has two nodes, the relationship between the nodal displacement and nodal force is given by

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

11

Fig. 6. Construction of discrete fracture element.

8 9 9 8 F 1x > > kðu1 u2 Þ > > > > > > > > > > < F 1y = < kðv 1 v 2 Þ = F¼ ; ¼ > > kðu2 u1 Þ > > F 2x > > > > > > > : : ; ; > F 2y kðv 2 v 1 Þ

ð40Þ

where F is nodal force vector in the global coordinate system, k denotes the stiffness parameter and a large value for intact material/interface is assumed (k = 1 1010 N/mm herein). The above expression can also be written as

2

0

k

k

0

0

k

k

0

where Kd is stiffness matrix, U ¼ ½ u1 respectively.

v1

F ¼ Kd U;

k

6 0 6 with Kd ¼ 6 4 k 0

0

3

k 7 7 7; 0 5

ð41Þ

k u2

v 2 T

is nodal displacement vector in the global coordinate system,

3.2. Numerical procedures The formations of failure modes A and B simply result from the surface crack deflection or penetration at the TC/BC interface. In contrast, the mode C should further deal with new crack onset at other weak locations after deflection/penetration occurs. Therefore, the numerical procedure mainly includes three aspects: (1) surface crack deflection or penetration at the TC/BC interface; (2) cracks propagation along the paths; (3) determination of new crack initiation at weak locations. The energetic method is employed to determine the crack initiation, which is implemented by inserting a short putative crack at the weak location and crack onset occurs when ERR at the putative crack tip is greater than the fracture toughness, as shown in Fig. 7. To avoid the interaction between putative cracks and main cracks, their fracture parameters are separately calculated under the same applied load. It means that the putative cracks are removed when evaluating fracture parameters for the main cracks. After the putative crack is judged to extend, it is denoted as a new part of the main crack. The numerical procedures are illustrated in Fig. 8. The main steps of the implementation are described as below: Step 1: Initiate the FE model. If the effect of residual stress is considered, determine the initially residual stress fields by the sequentially coupled thermo-mechanical method. Otherwise, the model is initially stress-free; Step 2: Impose a small load increment; Step 3: Judge surface crack deflection or penetration: Firstly, doubly deflect the pre-existing surface crack along TC/BC interface by an equal length ad and evaluate ERR Gd; thereafter, remove the deflected putative interface crack, which is implemented by recovering large stiffness for the discrete fracture elements; Afterwards, insert the putative penetration crack into bond coat by a length ap, evaluate ERR Gp, and then remove the putative penetration crack; Compare the ERR ratio with toughness ratio. If Gd/Gp > Cd/Cp, the crack deflection occurs and the interface cracks are marked as the main cracks. Otherwise, the crack penetration occurs and the vertical crack is marked as the main crack. Then, insert the interface crack or vertical crack according to the judgment; Step 4: Update the current crack tips information;

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

12

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 7. Putative cracks for determining crack onset at weak locations after the surface crack deflects or penetrates.

Fig. 8. Numerical procedure for fracture analysis of TBCs.

Step 5: Judge main crack propagation: evaluate the ERRs at the main crack tips and advance the main cracks if the ERRs are greater than the corresponding fracture toughnesses; Step 6: Judge whether new cracks already initiated at weak locations in the previous loading history. If yes, move to Step 8; otherwise, move to Step 7; Step 7: Determine cracks initiation at current load states: insert putative cracks at weak locations, where positions of the putative crack are illustrated in Fig. 7, and then evaluate the ERRs. If the ERRs are greater than the corresponding fracture toughnesses, crack onset occurs and the newly initiated cracks are marked as one of the main cracks. Otherwise, the putative crack fails to initiate, and remove the putative cracks; Step 8: Judge whether it is the end of the analysis. If yes, the analysis is finished. Otherwise, impose a load increment and return to Step 4. 3.3. Finite element models The feature of the periodical distributions of surface cracks makes it feasible to use the unit cell model in the analyses. All the layers are considered to be homogeneous, isotropic and linear elastic. Different surface crack densities are implemented by varying the width of the unit cell model (w = 250 lm, 500 lm, 750 lm and 1000 lm). The top coat is usually made of YSZ and deposited by air plasma spraying (APS) method for gas turbine application. The elastic modulus of the top coat, ETC, is around 20 GPa. Nevertheless, the bond coat may be fabricated by various techniques to Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

13

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 1 Geometrical and material properties used in the numerical analysis. Layer

Thickness

Elastic modulus

Poisson’s ratio

Top coat Bond coat Substrate

300 lm 100 lm 5 mm

20 GPa 20–240 GPa 210 GPa

0.1 0.3 0.3

Table 2 Elastic modulus of bond coat used in the numerical analyses and corresponding Dundurs’ parameters between top coat and bond coat (ETC is kept constant to be 20 GPa in all analyses). Case

1

2

3

4

5

6

7

EBC EBC/ETC

20 GPa 1 0.04 0.09

40 GPa 2 0.37 0.21

80 GPa 4 0.63 0.31

120 GPa 6 0.73 0.35

160 GPa 8 0.79 0.37

200 GPa 10 0.83 0.38

240 GPa 12 0.86 0.39

a b

obtain the improved oxidation resistance, which results in the variations of elastic modulus. Therefore, the effect of material properties is investigated by varying the elastic modulus of the bond coat, EBC, from 20 GPa to 240 GPa, whereas other material properties are kept constants. The material properties are summarized in Table 1 [47–49], and Dundurs’ parameters are listed in Table 2. Moreover, the interfacial fracture toughness of TC/BC interface, CTC/BC, is reported in the range of 10–100 J/ m2 [50–52], and it is specified to be CTC/BC = 50 J/m2 in this work. The tensile load and boundary condition are depicted in Fig. 3. Typical finite element mesh is shown in Fig. 9, in which fine meshes are used around the weak interfaces and very fine meshes are employed around the root of surface crack for the purpose of accurately capturing crack deflection/penetration behavior. The four-point isoparametric plane strain elements are utilized. The potential crack paths are placed along the interfaces and vertically throughout the coatings. The discrete fracture elements in the top coat are initially opened to represent the pre-existing surface crack. To ensure sufficient numerical accuracy in evaluating the fracture parameters, the integral domain is chosen that it contains at least four layers of elements around the crack tip and the crack face in the integral domain is straight.

Fig. 9. Typical finite element model for the unit cell model and illustration of the potential crack path.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

14

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

The numerical analyses are carried out by using commercial finite element package ABAQUS. The plane strain element and discrete fracture element are implemented by the user-defined element subroutine UEL. In the analyses, cracks can automatically propagate along their paths, branch into interface, and initiate at weak locations without remeshing. The applied displacement increment is set to a small value when the calculated ERR is greater than the fracture toughness. The evaluations of fracture parameters and controlling of automatic crack propagation is implemented in ABAQUS subroutine URDFIL.

4. Parametric studies 4.1. Putative crack length for surface crack deflection/penetration When applying the crack deflection criterion in Eq. (5), the lengths of putative crack extensions ad and ap need to be defined for numerically evaluating the ERRs at tips of the putative cracks. By dimensional considerations, HH analytically derived that the Gd/Gp ratio of a semi-infinite surface crack who has infinitesimal putative extensions ad and ap can be expressed as [28]

12k Gd ad ¼ f ða; b; rres Þ ; Gp ap

ð42Þ

where f ða; b; rres Þ is a function of Dundurs’ parameters and residual stresses rres . The relation in Eq. (42) suggests that Gd/Gp is independent of the putative crack lengths if the lengths are assumed to be identical with each other (ad = ap). However, the independence is valid only when the infinitesimal crack extension and the infinite model geometry are used. For the problems in this work, the models have finite geometry dimensions. The ad and ap may affect the Gd/Gp ratio even when ad = ap is assumed. Thus, the sensitivity of the putative crack lengths on the Gd/Gp ratio need to be clarified. By taking different elastic mismatches listed in Table 2, the Gd/Gp ratios calculated by numerical method as functions of a for various putative crack lengths are plotted in Fig. 10, where the width of the unit cell model w = 250 lm is used. For the present analyses, the ad = ap = 0.03 lm is the shortest lengths that could be implemented because there exists limit on the smallest mesh size in ABAQUS (the available smallest mesh size depends on the typical geometry dimensions of the finite element model). Moreover, the theoretical solution derived by HH [28], which is based on the infinitesimal crack extensions assumption, is also presented for comparison. It is seen that the putative crack lengths have marked influences on the ERR ratios. The dependence of Gd/Gp on the lengths becomes more significant as a increases towards 1.0. The calculated Gd/Gp gradually approaches the HH’s solution as ad = ap decreases. As a increases, the singularity of the surface crack decreases (forming weak singularity) and the region of K-dominant field diminishes (almost vanishes when a > 0.5). The vanishing K-dominant field indicates that even the shortest possible putative crack extension could be influenced by the remotely applied load [27]. In such case, the putative crack deflection or penetration is no longer a fracture process confined well within the boundary of stress singularity field of the preexisting surface crack. Therefore, for the cases of finite geometries in this work, the dependence of Gd/Gp on the ad and ap is unavoidable for large a. The proper values of the ad and ap still need to be defined.

Fig. 10. The Gd/Gp ratios as functions of a for various putative crack lengths.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

15

Nevertheless, how to determine the ad and ap is an issue that has not yet been well solved in previous studies. The physical meanings of the ad and ap were believed to relate to some microstructural features of the interface and the penetrated material, e.g. typical flaw sizes, misfit dislocations or grain boundary orientation [27,53,54]. Direct experimental measurement of such microstructural geometries is difficult, especially the ones associated with the interface, and suffers from problems in selecting the ‘correct’ type of the features which may vary for different materials. Experimental evidence showed that the HH’s solution overestimates the tendency for crack deflection at the interface [27,54]. It implies that the finite ad and ap should be used for practical materials. The infinitesimal ad and ap may be reached only for perfect interfaces and material, such that the HH’s solution could be expected to hold. Due to the difficulty in exact estimation of the putative crack lengths, the shortest lengths ad = ap = 0.03 lm are chosen in this work. It should be emphasized that the chosen lengths provide the upper bound for the failure maps in Fig. 10, in which the crack deflection is preferred to be predicted. Similarly, as will be seen in Section 5.1.3, the predictions for the TBCs would slightly prefer to the formation of failure mode C under large modulus ratio. It is recommended that general engineering design should account for the slight uncertainty between the failure mode B and mode C when the bond coat is much stiffer than the top coat. Furthermore, it is relevant to note that Martin and Leguillon et al. [53,55–57] attempted to propose more advanced criteria to address the issue of surface crack deflection/penetration at an interface. The energetic incremental theory was used rather than the energy release rate which is based on the energetic differential theory. The surface crack deflection/penetration is assessed by selecting the path which maximises the additional energy. This method is appealing because the ad and ap would not explicitly appear in the criterion and their influence on the result vanishes. The criterion is derived from an asymptotic analysis of the surface crack. Further work is needed to validate the applicability of the criterion in the case of finite model width. 4.2. Putative crack length for crack initiation at weak locations Under increasing applied load, the stored energy in the system is able to release at weak locations to induce the initiation of new cracks. The energetic method is employed to determine the new crack onset, in which the lengths of putative crack extensions adi or api need to be defined for evaluating the ERRs (as illustrated in Fig. 7). By choosing EBC = 200 GPa, the normalized driving forces at the TC/BC interface crack tips, GTC=BC =E0TC hTC e2 , for various w are plotted in Fig. 11, in which GTC/BC is the ERR at the tip of the TC/BC interface crack and E0TC ¼ ETC =ð1 m2 Þ for plane strain problems. The applied strain e is defined as the ratio of imposed displacement load to the width of the model w. It is seen that the driving forces rapidly increase to high levels and then gradually decrease as the crack lengths increase. The driving force of a crack is sensitive to its crack length. When employing the energetic method for determining the crack onset, it indicates that the selection of putative crack length has a significant influence on the calculated ERR, so as to affect the prediction of crack onset. Therefore, the putative crack length needs to be carefully selected to represent the actual crack initiation in the coatings. Experimental studies [58,59] reported that the applied strain for TC/BC interface crack onset is in the range of 0.5–1%. By taking GTC/BC = CTC/BC = 50 J/m2, hTC = 300 lm, ETC = 20 GPa, and 0.5% e 1%, the normalized ERR for TC/BC interface crack onset can be estimated to be 0.08 GTC=BC =E0TC hTC e2 0.33. It corresponds to the crack length 0.95 lm adi 4.4 lm (according to the results in Fig. 11). The selection of smaller putative crack length indicates the premature crack initiation and longer ‘‘pop-in” crack propagation, and vice versa. In order to obtain the relatively conservative results, the putative crack length for crack initiation is chosen to be adi = api = 2 lm in all analyses.

Fig. 11. Normalized driving forces for various surface crack spacing w.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

16

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 12. Variations of fracture parameters for various mesh sizes along the interface: (a) normalized ERR and (b) normalized crack length.

4.3. Mesh size along the interface For the nodal release approach, crack propagation is modeled by sequentially releasing the nodal pairs along the crack path. The distance between the neighboring nodal pairs, which is also the mesh size along the path, represents the increment of crack extension in the fracture process. The effect of mesh size along the interface on the variation of fracture parameter is examined herein. The evolutions of normalized ERR and length of the TC/BC interface crack for various mesh sizes are presented in Fig. 12, in which w = 250 lm and EBC = 200 GPa are used. It is seen that the evolutions show the zigzag patterns in the process of crack propagating along the TC/BC interface. The smoother results are obtained for the more refined meshes. It implies that the finite element mesh along the interface is the reason inducing the zigzag evolutions of the fracture parameters. Nonetheless, the mesh size has insignificant influence on the coatings spallation strain. This is easily understood because the linear elastic fracture process is assumed in the analysis. The increment of crack extension in the previous fracture process would not change the driving force at a crack tip under given applied load. In order to obtain the acceptable computational cost and sufficiently refined fracture parameters, the mesh size along the interfaces (including TC/BC interface and BC/Sub interface) is chosen to be 1 lm in all analyses.

5. Results and discussion 5.1. Initially stress-free coatings with brittle interface In this section, the coatings are assumed to be initially stress-free, and both the TC/BC and BC/Sub interfaces are considered to be brittle with kw ¼ 1.

5.1.1. Variation of fracture parameters under failure mode A In order to illustrate the effect of modulus mismatch, variations of fracture parameters for cases of a = 0.04 and 0.73 under w = 250 lm are compared in Fig. 13. The ERR at TC/BC interface crack tip (GTC/BC) and ERR at the putative crack tip e TC=BC ¼ GTC=BC =CTC=BC and G e BC ¼ GBC =CTC=BC . The in the bond coat (GBC) are both normalized by CTC/BC, which is written as G a TC=BC ¼ 2aTC=BC =w. The applied strain total length of TC/BC interface crack 2aTC=BC is normalized by w, which is written as e is defined as the ratio of imposed displacement to the width of the model w. Note that spallation of the coating is defined in the case that the total length of an interface crack reaches 0.9 times the model width (0.9w), similarly hereinafter. It is seen in Fig. 13 that variations of ERRs for both elastic mismatches are qualitatively similar. TC/BC interface crack initiates when GTC/BC reaches the fracture toughness. The small crack extension always results in a drop of driving force, indicating the stable crack propagation in these cases. The interface crack propagates fast at early stage but gradually approaches to steady state. This phenomenon coincides with the conclusions made by Xu et al. [60]. For putative crack in bond coat, its ERR is higher than that of interface crack (GBC GTC/BC). There is a decrease of GBC after the delamination onset. Larger modulus mismatch, or say stiffer bond coat, leads to a faster increase of GBC. It implies that new crack tends to initiate in the stiffer bond coat when the TC/BC interface crack is extending towards the boundary. Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

17

Fig. 13. Effect of modulus mismatch on the variations of ERRs under failure mode A (w = 250 lm and eonset is the interfacial delamination onset strain): (a) a = 0.04; (b) a = 0.73.

The effect of surface crack density is shown in Fig. 14, where a is kept to be 0.83. Note that the smaller w denotes the higher crack density. It is seen that the magnitude of GBC qualitatively rises with the increase of crack density, indicating that the bond coat tends to crack under high surface crack density. Particularly, for model width greater than 500 lm, GTC/BC abruptly jump to a large value so that the ‘pop-in’ crack extension occurs. Fig. 15 shows the variations of driving force with the delamination length at the pop-in crack extension moment. As the interface crack grows, the driving forces quickly increase to the peaks and then gradually decrease. The large driving force and long unstable delamination length may be one reason responsible for the short service life of TBCs with low surface crack density. Fig. 16(a)–(c) shows the rx, ry, and sxy stress fields under short crack (2aTC/BC/w = 0.1), and (d)–(f) shows the stress fields under long crack (2aTC/BC/w = 0.5). Phase angle w is an important parameter to characterize the component ratio of crack tip driving forces. The phase angle, | w|, as a function of delamination length is plotted in Fig. 17. The w is independent of applied load for specific crack length and is observed always greater than 45°. Therefore, the shear traction plays a dominant role during the delamination propagation. It is interesting that K1 is observed to gradually decrease from positive to zero, which corresponds to the increase of | w| to 90°, and then become negative thereafter to cause the change of sign of w and decrease of |w|. It is known that negative mode I SIF for homogeneous material crack usually indicates the penetration of two crack faces. However, it is surprising that the interface crack faces are found to always open even when the K1 is negative. It could be mathematically explained from the definition of crack flank displacements behind the tip, namely Eq. (29), where the normal displacement d2 depends on both K1 and K2 for specific distance r. It is possible to obtain a positive d2 with negative K1 and appropriate K2. This phenomenon could be considered as a peculiar feature for the interfacial fracture problem. Moreover, the delamination onset and coating spallation are both delayed for the TBCs with higher surface crack density, as shown in Fig. 18. It can be concluded that the interfacial fracture resistance is enhanced by introducing dense vertical cracks. The modulus mismatch has little influence on the crack onset strains as well as the spallation strain when the modulus ratio is large enough (it is EBC/ETC 4 herein). 5.1.2. Variations of fracture parameters under failure modes B and C This section aims to study the variations of fracture parameters under failure modes B and C. Effect of BC/Sub interfacial fracture toughness, CBC/Sub, on the fracture behavior are discussed. The brittle and tough BC/Sub interfaces are assumed when the toughness ratio, CBC/Sub/CTC/BC, equals 2 and 20, respectively. To ensure the appearance of failure mode B, the fracture toughness of bond coat is chosen to be CBC = 50 J/m2, as will be explained by the failure map in Section 5.1.3. Results for cases of brittle and tough BC/Sub interfaces are compared in Fig. 19, where ERRs are normalized by their corresponding fracture toughnesses, and the interface crack length is normalized by w. Under failure mode B, the rapid release of energy leads to unstable propagation of the vertical crack through the bond coat. Then, the vertical crack deflects along the BC/Sub interface and steadily extends along the BC/Sub interface. The steady delamination along tough interface is important for failure mode B since the strain tolerance is significantly improved. After the vertical crack penetrates into the bond coat, the ERR at putative TC/BC interface crack tip keeps a low value. This phenomenon indicates that TC/BC interfacial delamination is unlikely to occur for relatively brittle BC/Sub interface after the surface crack penetrating into bond coat. The stress distributions at different moments are given in Fig. 20. To ensure the formation of failure mode C, the fracture toughness of bond coat is chosen to be a medium value of CBC = 500 J/m2. As shown in Figs. 21 and 22, ERR at TC/BC interface crack tip reaches its fracture toughness earlier than that of putative crack in bond coat, which indicates that the surface crack tends to deflect. ERR at the putative crack tip in bond coat Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

18

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 14. Effect of surface crack density on the variations of ERRs under failure mode A (a = 0.83 and eonset is the interfacial delamination onset strain): (a) w = 250 lm; (b) w = 500 lm; (c) w = 750 lm; (d) w = 1000 lm.

Fig. 15. Energy release rate as a function of delamination length when unstable crack propagation occurs.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

19

Fig. 16. Local stress distributions under failure mode A when crack lengths are: (a)–(c) 2aTC/BC/w = 0.1; (d)–(f) 2aTC/BC/w = 0.5 (w = 250 lm, a = 0.83, and deformation scale factor is 5).

Fig. 17. Phase angle w as a function of normalized TC/BC interface crack length under failure mode A.

generally increases with the propagation of TC/BC interface crack. Different from that under failure mode A, new crack initiates in bond coat before the TC/BC interface crack reaches the critical length. For the case of brittle BC/Sub interface, see Fig. 21, newly initiated crack is observed to unstably and rapidly extend throughout the bond coat and BC/Sub interface until its complete spallation. The large energy stored in bond coat leads to strong crack tip driving force and rapid release of energy along the brittle interface. However, as shown in Fig. 22, tough BC/Sub interface can make the unstable delamination become stable when the interface crack length gets long enough.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

20

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 18. The delamination onset and coating spallation strains for TC/BC interface under failure mode A, regarding different surface crack densities and modulus mismatches.

Fig. 19. Variation of fracture parameters under failure mode B with w = 250 lm and a = 0.83: (a) brittle BC/Sub interface CBC/Sub/CTC/BC = 2; (b) tough BC/ Sub interface CBC/Sub/CTC/BC = 20.

Nonetheless, this behavior also reflects that the driving force at BC/Sub interface crack tip under mode C is generally stronger than the one under mode B. Thus, larger CBC/Sub is required for mode C than mode B to ensure the occurrence of steady BC/ Sub interfacial delamination. Moreover, in both brittle and tough cases, the TC/BC interface cracks are arrested after new crack initiates in bond coat, which implies that the new cracks become the only paths to release the energy. The stress distributions under failure mode C at different moments are given in Fig. 23, where crack tips are marked in the figures. Variations of phase angle at tough BC/sub interface under modes B and C are plotted in Fig. 24. The phase angle of BC/ Sub interface crack increases with the propagation of interfacial delamination, and the one of TC/BC interface crack under mode C varies similarly with that of mode A. Shear tractions are mainly responsible for interfacial delamination at both the weak interfaces. To get a better understanding of the pop-in crack propagation behaviors in the bond coat as well as at BC/Sub interface, the driving forces are compared and illustrated in Fig. 25. It is seen that the driving forces of vertical crack in bond coat generally increase with the crack length, which implies that the vertical crack unstably propagates faster and faster. The driving force under mode C is higher than that under mode B, due to the greater CBC of mode C and thus the larger stored energy when new crack onset occurs. Conversely, the driving forces at BC/Sub interface crack tips gradually decrease with the delamination length, which indicates the interface crack could be arrested when its length is long enough. Tough BC/Sub interface is observed to benefit in its weaker driving force and the shorter unstable delamination length, resulting in the better strain tolerance. Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

21

Fig. 20. Local stress distributions under failure mode B when crack lengths are: (a)–(c) aBC/hBC = 0.5; (d)–(f) 2aBC/Sub/w = 0.1; (g)–(i) 2aBC/Sub/w = 0.5 (w = 250 lm, a = 0.83, CBC/Sub/CTC/BC = 20, and deformation scale factor is 2. aBC is the crack length in bond coat, and aBC/Sub is the half crack length along BC/Sub interface).

Fig. 21. Variation of fracture parameters under failure mode C with brittle BC/Sub interface (w = 250 lm, a = 0.83, and CBC/Sub/CTC/BC = 2).

Furthermore, the crack onset strain and coating spallation strain under failure modes A, B and C (in the same case of w = 250 lm, a = 0.83) are shown in Fig. 26. Note that the ‘crack onset’ means TC/BC interface cracks or vertical crack in bond coat starts to grow from the root of pre-existing surface crack. It is seen that the CBC/Sub has a significant influence on the coating

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

22

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 22. Variation of fracture parameters under failure mode C with tough BC/Sub interface (w = 250 lm, a = 0.83 and CBC/Sub/CTC/BC = 20): (a) variables for TC/BC interface crack and putative crack in the bond coat; (b) ERR and crack length for BC/Sub interface crack.

Fig. 23. Local stress distributions under failure mode C when crack lengths are: (a)–(c) aBC/hBC = 0.5; (d)–(f) 2aBC/Sub/w = 0.1; (g)–(i) 2aBC/Sub/w = 0.5 (w = 250 lm, a = 0.83, CBC/Sub/CTC/BC = 20, and deformation scale factor is 1).

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

23

Fig. 24. Phase angle |w| as a function of normalized interface crack length under failure modes B and C for the case of w = 250 lm, a = 0.83 and CBC/Sub/CTC/ = 20.

BC

Fig. 25. Comparison of driving force for vertical crack in bond coat and interface crack at BC/Sub interface under failure mode B and C.

spallation strain for failure modes B and C. The mode C has the largest crack onset strain (equal with that of mode A) and coating spallation strain. Above observations suggest that spallation of dense vertically cracked TBCs is essentially the competition between vertical crack deflection and penetration at TC/BC interface, as well as the crack initiation in bond coat. Controlling the formation of failure modes could be implemented by selecting proper mechanical properties of the coating components. Of the many influence factors, fracture toughness ratio of bond coat to TC/BC interface is a dominant one. Generally, tough bond coat tends to induce the TBCs spalling from TC/BC interface, the brittle one is more likely to fracture at BC/Sub interface, whereas delamination at both the interfaces forms when the bond coat fracture toughness is medium. Nonetheless, the increasing crack tip driving force and the pop-in crack extension in bond coat indicate that it fails to provide resistance to crack propagation once the vertical crack starts to grow therein. When this type of cracking happens, BC/Sub interface becomes the only path to resist crack growth, and thus it should tough enough to ensure the delamination propagation is stable. From the point of view of toughening, employing coating with high-density surface cracks in the top coat, the medium toughness of bond coat and tough BC/Sub interface can provide good fracture resistance and strain tolerance. 5.1.3. Failure maps under various surface crack densities Failure map of TBCs under various surface crack densities is numerically predicted and illustrated in Fig. 27, in which the modulus ratio, EBC/ETC, and fracture toughness ratio, CBC/CTC/BC, are considered as design variables. Previous results in this work show that new crack may initiate in bond coat after the surface crack deflects, but it is unlikely to delaminate at TC/BC interface after the surface crack penetrates into bond coat. Therefore, when constructing the failure maps, the boundary for Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

24

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 26. The crack onset strain and coating spallation strain under failure mode A, B, and C.

Fig. 27. Failure maps for various surface crack densities.

surface crack penetration could be simply determined by energy criterion given in Eq. (5), while the boundary for crack deflection is obtained by calculating the maximum ERR of putative crack in bond coat before the coating spalling off. In the maps, failure mode C vanishes when the modulus ratio is low, which is caused by the slow increasing ERR at the putative crack tip in compliance bond coat. Surface crack density is shown to have a significant influence on the formation of failure modes. As the density increases, regions of mode C enlarge and thus the failure mode tends to change from mode A to mode C. In other word, for the given material properties, it is more likely to crack in bond coat and at BC/Sub interface for coatings with high surface crack density. Fortunately, the fracture toughness of metallic-metallic BC/Sub interface is usually much higher than that of the ceramic-metallic TC/BC interface. Thus, introducing high-density surface cracks in TBCs benefits in reducing the coating stiffness and thus reliefs the interfacial mismatch stress. Also, it is easier to induce the coating spalling from the tougher interface and releasing energy from multiple paths so that the fracture resistance could be further improved.

5.2. Effect of mixed-mode interface toughness The interface toughness is usually a function of the mode mixity acting on the interface, as expressed in Eq. (39). The effect of mixed-mode interface toughness on the formation of failure mode is studied in this section. For simplicity, the parameter kw of TC/BC interface is assumed to be identical with the one of BC/Sub interface. The kw = 1.0 is used to represent the ideally brittle interface, while kw = 0.5 and 0.25 are utilized to study the mode mixity dependent interfacial fracture process. Moreover, the coatings are assumed to be initially stress-free. Fig. 28 plots the typical evolutions of fracture parameters under failure mode A for kw – 1.0 (w = 250 lm, a = 0.83, and kw = 0.5 is used herein). The ERRs at TC/BC interface crack tips are normalized by mixed-mode toughness CTC/BC(w), while the Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

25

Fig. 28. Evolutions of fracture parameters under failure mode A for kw = 0.5 (w = 250 lm and a = 0.83).

ERR at the putative crack tip in the bond coat is normalized by mode 1 TC/BC toughness C1TC=BC ¼ 50 J=m2 . During the interfacial fracture process, a pop-in crack extension occurs when the phase angle |w| reaches 90°, leading to the decrease of ERR at the putative crack tip in bond coat. The fracture processes for cases kw – 1.0 differ from the one for the ideally brittle interface (kw = 1.0) whose pop-in extension only arises at the root of the surface crack. The difference is caused by the nonuniform variations of the interface toughness for cases kw – 1.0. According to the variations of |w| along the TC/BC interface, as shown in Fig. 17, the |w| gradually increases to 90° and then decreases, which corresponds to the increase and then decrease of CTC/BC(w). The interface crack with |w| = 90° has the highest toughness and the strongest crack tip driving force. Once the interface crack extends forward with an increment, some stored strain energy releases from the coatings. However, the remained crack tip driving force still exceeds the interface toughness at the new crack tip, such that the pop-in extension occurs. The fracture processes under failure mode B for cases kw – 1.0 are qualitatively similar to the ones presented in Fig. 19. None pop-in crack extension is observed during the fracture of BC/Sub interface because the interface toughness generally increases as the crack extends along the interface (see Fig. 24). Moreover, the fracture processes under failure mode C for cases kw – 1.0 are also qualitatively similar to the ones shown in Figs. 21 and 22, except that a pop-in extension arises when the |w| of TC/BC interface crack reaches 90°. The crack onset strains and coating spallation strains for various kw are compared in Fig. 29, in which the case of w = 250 lm and a = 0.83 is used. The fracture toughness CBC = 50 J/m2 is chosen to ensure the formation of mode B, and CBC = 500 J/ m2 is employed to induce the formation of mode C. Moreover, the tough BC/Sub interface is assumed for modes B and C (CBC/ Sub/CTC/BC = 20). It is seen that the parameter kw has a significant influence on the crack onset and coatings spallation strains. Generally, the strain tolerance of the coatings is enhanced with the decrease of kw, due to the enhancement of the interface toughness for |w| – 0. The strain tolerance under failure mode C is superior to the ones under modes A and B, which confirms that the introduction of mode C is effective in improving the fracture resistance of the coatings. When the mixed-mode interface toughness is considered, the surface crack deflection criterion in Eq. (5) is still valid, but the toughness Cd should be recognized as C(w). Thus, the criterion can be rewritten as

C1 Gd < : Cp Gp f1 þ tan2 ½ð1 kw Þwg

ð43Þ

Considering the expression in Eq. (43), the failure maps for various kw are presented in Fig. 30, in which the toughness ratio CBC =C1TC=BC is chosen as a variable. The parameter kw has a significant influence on the formation of failure modes. As kw decrease, it is easier to crack in the bond coat to form failure modes B and C, and the tougher bond coat is required to induce the formation of mode A. 5.3. Effect of residual stresses TBCs for gas turbine application are generally fabricated by thermal spraying method at high-temperature. After cooling down to room temperature, residual stresses are introduced into the coatings due to the thermal expansion mismatches among different layers. Effect of residual stresses on the formation of failure modes is discussed in this section. Herein, Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

26

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 29. Crack onset strains and coating spallation strains for various kw under (a) failure mode A; (b) failure mode B; (c) failure mode C (w = 250 lm, a = 0.83, and CBC/Sub/CTC/BC = 20).

Fig. 30. Failure maps for the case of w = 250 lm under various kw (C1TC=BC is the mode 1 TC/BC interface toughness).

6 typical thermal expansion coefficients of the top coat, bond coat, and substrate are, respectively, given as ath °C1, TC = 9 10 6 1 6 1 th th aBC = 13.6 10 °C , and asub = 12 10 °C [61]. The interfaces are assumed to be ideally brittle with kw = 1.0. Fig. 31 plots typical variations of fracture parameters under failure mode A for coatings containing residual stresses (w = 250 lm and a = 0.83). The residual stresses are introduced by steadily cooling the coatings with initial temperature T0 = 1000 °C to room temperature (20 °C). Afterwards, mechanical tensile loading is imposed to induce the fracture and spallation of the coatings. The definition of applied strain is the same with that in Section 5.1, which is the ratio of imposed displacement to the original width of the model. For simplicity of the analyses, the residual stresses induced contact between the vertical crack faces is not considered. This simplicity is believed to have negligible influence on the failure map as well as the strain tolerance.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

27

Fig. 31. Evolutions of fracture parameters under failure mode A with initially residual stresses (w = 250 lm and a = 0.83). The coatings are initially cooled down from T0 = 1000 °C to room temperature (20 °C) prior to the mechanical loading.

Fig. 32. Crack onset strains and coating spallation strains for various initially stress states under (a) failure mode A; (b) failure mode B (T0 = 20 °C denotes the initially stress-free state, and T0 = 600 °C and 1000 °C denotes that the coatings are cooled down from 600 °C and 1000 °C, respectively).

In Fig. 31, it is seen that the evolution of ERR at the TC/BC interface crack tip experiences a decrease from an initial state and then an increase, whereas the ERR of putative crack in bond coat tends to increases as the TC/BC interface crack extends. Above difference in variations of ERRs result from the superposition of applied tensile load on the different residual stressed coating layers. After cooling down to room temperature, tensile residual stresses (in the loading direction) arise in bond coat and compressive stresses develop in the top coat (because aTC < aBC herein). As the imposed mechanical loading increases, the stresses in top coat gradually change from compressive to tensile states, leading to the decrease and then increase of ERR at TC/BC interface crack tip, whereas the tensile stress in bond coat is gradually enhanced, inducing the increase of ERR at putative crack tip in bond coat. Under failure mode C, the fracture of TC/BC interface is qualitatively similar to the one under mode A, and fracture of bond coat and BC/Sub interface is similar with the one of initially stress-free coatings. The crack onset and coating spallation strains under failure modes A and C are presented in Fig. 32, in which the results of initially stress-free and residual stress states are compared (the material properties for different initial stress states are the same). Under mode A, the residual stresses benefit in slightly delaying the cracking and spalling of the coatings, which attributes to the initially compressive stresses in the top coat. However, the crack onset and spallation strains of the residual stressed coatings under mode C is reduced, due to the enhanced tensile stresses in bond coat. Nonetheless, mode C still has better strain tolerance than mode A. It suggests that design of proper residual stresses for a failure mode facilitates the improvement of strain tolerance. Under failure mode C, the residual stresses should be reduced as much as possible. Failure maps of TBCs under various initial stress states are plotted in Fig. 33. Residual stresses have a significant influence on the surface crack deflection/penetration as well as the new crack onset behavior in bond coat. It seems unlikely to form

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

28

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 33. Failure maps for the case of w = 250 lm under various initially residual stresses.

the failure mode B for a stiff bond coat, e.g. for the case of EBC/ETC = 10 and T0 = 1000 °C, the toughness ratio CBC/CTC/BC is as low as 0.55 and it is not practical. The regions of failure mode C enlarge as the residual stresses enhance. Hence, it is easier to form the mode C when high residual stresses arise in the coatings. 6. Conclusions In this work, mechanisms governing the formation of failure modes in dense vertically cracked thermal barrier coatings (TBCs) was numerically investigated. An energy-based criterion was employed to deal with the vertical crack deflection/penetration behavior at a bimaterial interface, and an interaction integral method for discontinuous materials was utilized to calculate the fracture parameters. Three-layer unit cell models were developed to study the cracking behavior under different failure modes. Based on the numerical results, it is demonstrated that the top coat/bond coat interfacial fracture resistance can be enhanced by introducing high-density surface cracks. Fracture toughness ratio of bond coat to top coat/bond coat interface is a dominant factor governing the transition of failure modes. Generally, tough bond coat tends to induce the TBCs spalling from TC/BC interface, the brittle one is more likely to induce the fracture occurring at bond coat/substrate interface, and delamination at both weak interfaces forms when the fracture toughness is medium. Nonetheless, bond coat fails to provide resistant to crack propagation due to the increasing driving force and unstable crack propagation therein. With the increase of surface crack density, cracking tends to occur in bond coat, leading to delamination at bond coat/substrate interface or delamination at both weak interfaces. Introducing high-density surface cracks in TBCs benefits in reducing the coating stiffness and thus reliefs the interfacial mismatch stress. Also, it is easier to induce the coating release energy from multiple paths so that the fracture resistance could be further improved to obtain more strain-tolerant TBCs. Furthermore, the mixed-mode interface toughness have significant influence on the formation of failure modes. The stronger dependence of interface toughness on the mode mixity benefits in delaying the cracking in the coatings. Moreover, initially residual stresses have opposite influence on the strain tolerance of the coatings under different failure modes. Acknowledgements The authors would like to thank Dr. Yixiu Shu and Prof. Yazhi Li at Northwestern Polytechnical University for valuable discussions on the fracture analyses. Dr. Peng Jiang and Bowen Lv at Xi’an Jiaotong University are gratefully acknowledged for their useful comments on the manuscript. This work is supported by China 973 Program (2013CB035700) and National Natural Science Foundation of China (11472204 and 11602188). References [1] Padture NP, Gell M, Jordan EH. Thermal barrier coatings for gas-turbine engine applications. Science 2002;296:280–4. [2] Evans AG, Mumm D, Hutchinson J, Meier G, Pettit F. Mechanisms controlling the durability of thermal barrier coatings. Prog Mater Sci 2001;46:505–53. [3] Fan X, Jiang W, Li J, Suo T, Wang T, Xu R. Numerical study on interfacial delamination of thermal barrier coatings with multiple separations. Surf Coat Technol 2014;244:117–22. [4] Cao X, Vassen R, Stoever D. Ceramic materials for thermal barrier coatings. J Eur Ceram Soc 2004;24:1–10. [5] Zhou B, Kokini K. Effect of surface pre-crack morphology on the fracture of thermal barrier coatings under thermal shock. Acta Mater 2004;52:4189–97. [6] Zhang W, Fan X, Wang T. The surface cracking behavior in air plasma sprayed thermal barrier coating system incorporating interface roughness effect. Appl Surf Sci 2011;258:811–7. [7] Fan X, Xu R, Zhang W, Wang T. Effect of periodic surface cracks on the interfacial fracture of thermal barrier coating system. Appl Surf Sci 2012;258:9816–23. [8] Zhou B, Kokini K. Effect of preexisting surface cracks on the interfacial thermal fracture of thermal barrier coatings: an experimental study. Surf Coat Technol 2004;187:17–25. [9] Taylor T, Appleby D, Weatherill A, Griffiths J. Plasma-sprayed yttria-stabilized zirconia coatings: structure-property relationships. Surf Coat Technol 1990;43:470–80.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

29

[10] Gell M, Xie L, Ma X, Jordan EH, Padture NP. Highly durable thermal barrier coatings made by the solution precursor plasma spray process. Surf Coat Technol 2004;177:97–102. [11] Guo H, Vaßen R, Stöver D. Atmospheric plasma sprayed thick thermal barrier coatings with high segmentation crack density. Surf Coat Technol 2004;186:353–63. [12] Kokini K, Banerjee A, Taylor TA. Thermal fracture of interfaces in precracked thermal barrier coatings. Mater Sci Eng, A 2002;323:70–82. [13] Wu C-W, Chen G-N, Zhang K, Luo G-X, Liang N-G. The effect of periodic segmentation cracks on the interfacial debonding: study on interfacial stresses. Surf Coat Technol 2006;201:287–91. [14] Chen X, Chen G, Zhang K, Luo G, Xiao J. On the thermally induced interfacial delamination of a segmented coating. Surf Coat Technol 2008;202:2878–84. [15] Fan X, Zhang W, Wang T, Liu G, Zhang J. Investigation on periodic cracking of elastic film/substrate system by the extended finite element method. Appl Surf Sci 2011;257:6718–24. [16] Zhou B, Kokini K. Effect of pre-existing surface crack morphology on the interfacial thermal fracture of thermal barrier coatings: a numerical study. Mater Sci Eng, A 2003;348:271–9. [17] Zhu W, Yang L, Guo JW, Zhou YC, Lu C. Numerical study on interaction of surface cracking and interfacial delamination in thermal barrier coatings under tension. Appl Surf Sci 2014;315:292–8. [18] Mei H, Pang Y, Huang R. Influence of interfacial delamination on channel cracking of elastic thin films. Int J Fract 2007;148:331–42. [19] Zhu W, Yang L, Guo J, Zhou Y, Lu C. Numerical study on interaction of surface cracking and interfacial delamination in thermal barrier coatings under tension. Appl Surf Sci 2014;315:292–8. [20] Fan X, Zhang W, Wang T, Sun Q. The effect of thermally grown oxide on multiple surface cracking in air plasma sprayed thermal barrier coating system. Surf Coat Technol 2012;208:7–13. [21] Qian L, Zhu S, Kagawa Y, Kubo T. Tensile damage evolution behavior in plasma-sprayed thermal barrier coating system. Surf Coat Technol 2003;173:178–84. [22] Li X, Liang L, Xie J, Chen L, Wei Y. Thickness-dependent fracture characteristics of ceramic coatings bonded on the alloy substrates. Surf Coat Technol 2014;258:1039–47. [23] Chen ZB, Wang ZG, Zhu SJ. Tensile fracture behavior of thermal barrier coatings on superalloy. Surf Coat Technol 2011;205:3931–8. [24] Ming-Yuan H, Hutchinson JW. Crack deflection at an interface between dissimilar elastic materials. Int J Solids Struct 1989;25:1053–67. [25] Dai-Heng C. A crack normal to and terminating at a bimaterial interface. Eng Fract Mech 1994;49:517–32. [26] Romeo A, Ballarini R. A cohesive zone model for cracks terminating at a bimaterial interface. Int J Solids Struct 1997;34:1307–26. [27] Lee W, Yoo Y-H, Shin H. Reconsideration of crack deflection at planar interfaces in layered systems. Compos Sci Technol 2004;64:2415–23. [28] He MY, Evans AG, Hutchinson JW. Crack deflection at an interface between dissimilar elastic materials: role of residual stresses. Int J Solids Struct 1994;31:3443–55. [29] Martínez D, Gupta V. Energy criterion for crack deflection at an interface between two orthotropic media. J Mech Phys Solids 1994;42:1247–71. [30] Stern M, Becker E, Dunham R. A contour integral computation of mixed-mode stress intensity factors. Int J Fract 1976;12:359–68. [31] Yau J, Wang S, Corten H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. J Appl Mech 1980;47:335–41. [32] Wang S, Yau J, Corten H. A mixed-mode crack analysis of rectilinear anisotropic solids using conservation laws of elasticity. Int J Fract 1980;16:247–59. [33] Raju I, Shivakumar K. An equivalent domain integral method in the two-dimensional analysis of mixed mode crack problems. Eng Fract Mech 1990;37:707–25. [34] Johnson J, Qu J. An interaction integral method for computing mixed mode stress intensity factors for curved bimaterial interface cracks in nonuniform temperature fields. Eng Fract Mech 2007;74:2282–91. [35] Sukumar N, Huang Z, Prévost JH, Suo Z. Partition of unity enrichment for bimaterial interface cracks. Int J Numer Meth Eng 2004;59:1075–102. [36] Gosz M, Moran B. An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions. Eng Fract Mech 2002;69:299–319. [37] Okada H, Ohata S. Three-dimensional J-integral evaluation for cracks with arbitrary curvatures and kinks based on domain integral method for quadratic tetrahedral finite element. Eng Fract Mech 2013;109:58–77. [38] Dolbow J, Gosz M. On the computation of mixed-mode stress intensity factors in functionally graded materials. Int J Solids Struct 2002;39:2557–74. [39] Walters MC, Paulino GH, Dodds Jr RH. Computation of mixed-mode stress intensity factors for cracks in three-dimensional functionally graded solids. J Eng Mech 2006;132:1–15. [40] Yu H, Wu L, Guo L, Du S, He Q. Investigation of mixed-mode stress intensity factors for nonhomogeneous materials using an interaction integral method. Int J Solids Struct 2009;46:3710–24. [41] Yu H, Wu L, Guo L, He Q, Du S. Interaction integral method for the interfacial fracture problems of two nonhomogeneous materials. Mech Mater 2010;42:435–50. [42] Guo L, Guo F, Yu H, Zhang L. An interaction energy integral method for nonhomogeneous materials with interfaces under thermal loading. Int J Solids Struct 2012;49:355–65. [43] Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech 1968;35:379–86. [44] Kim J-H, Paulino GH. Consistent formulations of the interaction integral method for fracture of functionally graded materials. J Appl Mech 2005;72:351–64. [45] Hutchinson JW, Suo Z. Mixed Mode Cracking in Layered Materials. Adv Appl Mech 1992;29:63–191. [46] Banks-Sills L, Boniface V, Eliasi R. Development of a methodology for determination of interface fracture toughness of laminate composites—the 0 /90 pair. Int J Solids Struct 2005;42:663–80. [47] Ranjbar-Far M, Absi J, Shahidi S, Mariaux G. Impact of the non-homogenous temperature distribution and the coatings process modeling on the thermal barrier coatings system. Mater Des 2011;32:728–35. [48] Cao XQ, Vassen R, Stoever D. Ceramic materials for thermal barrier coatings. J Eur Ceram Soc 2004;24:1–10. [49] Evans AG, Mumm DR, Hutchinson JW, Meier GH, Pettit FS. Mechanisms controlling the durability of thermal barrier coatings. Prog Mater Sci 2001;46:505–53. [50] Arai M, Okajima Y, Kishimoto K. Mixed-mode interfacial fracture toughness for thermal barrier coating. Eng Fract Mech 2007;74:2055–69. [51] Kim SS, Liu YF, Kagawa Y. Evaluation of interfacial mechanical properties under shear loading in EB-PVD TBCs by the pushout method. Acta Mater 2007;55:3771–81. [52] Xu ZH, Yang Y, Huang P, Li X. Determination of interfacial properties of thermal barrier coatings by shear test and inverse finite element method. Acta Mater 2010;58:5972–9. [53] Leguillon D, Lacroix C, Martin E. Interface debonding ahead of a primary crack. J Mech Phys Solids 2000;48:2137–61. [54] Ahn B, Curtin W, Parthasarathy T, Dutton R. Criteria for crack deflection/penetration criteria for fiber-reinforced ceramic matrix composites. Compos Sci Technol 1998;58:1775–84. [55] Martin E, Leguillon D, Lacroix C. A revisited criterion for crack deflection at an interface in a brittle bimaterial. Compos Sci Technol 2001;61:1671–9. [56] Leguillon D, Lacroix C, Martin E. Crack deflection by an interface–asymptotics of the residual thermal stresses. Int J Solids Struct 2001;38:7423–45. [57] Martin E, Poitou B, Leguillon D, Gatt J. Competition between deflection and penetration at an interface in the vicinity of a main crack. Int J Fract 2008;151:247–68. [58] Yang L, Zhong ZC, You J, Zhang QM, Zhou YC, Tang WZ. Acoustic emission evaluation of fracture characteristics in thermal barrier coatings under bending. Surf Coat Technol 2013;232:710–8.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037

30

B. Li et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

[59] Yang L, Zhong Z, Zhou Y, Zhu W, Zhang Z, Cai C, et al. Acoustic emission assessment of interface cracking in thermal barrier coatings. Acta Mech Sin 2016;32:342–8. [60] Xu R, Fan X, Zhang W, Song Y, Wang T. Effects of geometrical and material parameters of top and bond coats on the interfacial fracture in thermal barrier coating system. Mater Des 2013;47:566–74. [61] Li B, Fan X, Zhou K, Wang T. Effect of oxide growth on the stress development in double-ceramic-layer thermal barrier coatings. Ceram Int 2017;43:14763–74.

Please cite this article in press as: Li B et al. Mechanisms governing the failure modes of dense vertically cracked thermal barrier coatings. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.11.037