Oxide growth and damage evolution in thermal barrier coatings

Engineering Fracture Mechanics 78 (2011) 2139–2152

Contents lists available at ScienceDirect

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

Oxide growth and damage evolution in thermal barrier coatings T.S. Hille a,b,1, S. Turteltaub a, A.S.J. Suiker a,⇑ a b

Delft University of Technology, Faculty of Aerospace Engineering, P.O. Box 5058, 2600 GB Delft, The Netherlands Materials Innovation Institute (M2i), Mekelweg 2, 2628 CD Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 29 June 2010 Received in revised form 5 April 2011 Accepted 6 April 2011 Available online 13 April 2011 Keywords: TBC systems Multi-physics Oxide growth Partition-of-unity method Cohesive law

a b s t r a c t Cracking in thermal barrier coatings (TBC) is triggered by the development of a thermallygrown oxide (TGO) layer that develops during thermal cycling from the oxidation of aluminum present in the bond coat (BC). In the present communication a numerical model is presented that describes the interactive development of the TGO morphology and the fracture processes in TBC systems in a mesh-independent way. The evolution of the TGO–BC mixture zone is described by an oxygen diffusion–reaction model. The partition-of-unity method is employed for the simulation of discrete cracking, where cracks can nucleate and propagate across finite elements at arbitrary locations and orientations. The validity of the model is demonstrated through the analysis of a representative TBC system subjected to a specific thermal cycling process. A parametric analysis demonstrates the sensitivity of the response of the TBC system to the fracture strength of the top coat. The simulation results indicate that cracks appear primarily at the current location of the BC–TGO interface and may nucleate at early stages of thermal cycling. These results are in good agreement with recent experimental observations reported in the literature. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Thermal barrier coatings (TBCs) are applied to the metallic components of high-temperature gas turbines used in aircrafts and power stations [1,2]. A TBC maintains the integrity of the structural core of a turbine by protecting it against oxidation during service. Moreover, a TBC provides a significant temperature difference across its thickness, which makes it possible to operate the turbine more efficiently at higher temperatures. In order to adequately perform these two tasks, a TBC system is typically composed of three layers: (i) the bond coat (BC), (ii) the thermally-grown oxide (TGO) and (iii) the top coat (TC). The BC, which in this study is represented by a NiCoCrAlY alloy, adheres the coating system to the substrate (e.g., a Ni-based superalloy). The TGO, a relatively thin layer of alumina (Al2O3), appears on top of the BC through oxidation under service conditions due to oxygen diffusion through the TC. Consequently, further oxygen diffusion to the substrate is reduced significantly, which avoids the oxidation of the substrate. The TC, an yttria-stabilized zirconia, provides a substantial temperature difference across the TBC thickness due to its relatively low heat conductivity. The TC considered in the present study is taken to be applied by electron beam physical vapor deposition (EB-PVD) [3], which results in a relatively flat contact area with the structure below. In aircraft engines, the thermo-mechanical loading applied to TBC systems is characterized by the specific flight cycles. At elevated operation temperatures, the thickness of the TGO increases as aluminum from the BC is continuously oxidized. When the engine is shut down, large in-plane compressive stresses arise in the TGO [4]. The compressive stresses during

⇑ Corresponding author. 1

E-mail address: [email protected] (A.S.J. Suiker). Present address: Siemens AG Energy Sector, Mellinghofer Str. 55, 45473 Mülheim an der Ruhr, Germany.

0013-7944/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2011.04.003

2140

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

Nomenclature Latin symbols c oxygen concentration C elastic stiffness tensor d length parameter in the cohesive model d1, d2 depths to characterize the undulation D oxygen diffusion coefficient D⁄ oxygen diffusion coefficient of undamaged material D elastic compliance tensor E Young’s modulus E energy dissipated during crack opening fd damage loading function fp rate-independent yield function effective fracture strength ft g rational function that characterizes the traction–separation law Gt effective fracture toughness h thickness of a layer H height of the computational domain I second-order identity tensor K elastic stiffness at crack face contact l1, l2 lengths to characterize the undulation L length of the computational domain m damage exponent of diffusivity M oxygen molarity in alumina n volume fraction n normal outward unit vector pO2 partial oxygen pressure R reaction rate in oxidation process S oxygen sink t time t tractions on a cohesive surface T effective traction on a cohesive surface u displacement x position in Euclidean space X geometrical interpolation variable for characterizing the anisotropy of the top coat Greek symbols a coefficient of thermal expansion b weighing factor for mode I and II cracking contributions C gamma function c shape parameter for the cohesive law d (auxiliary) crack opening  d actual crack opening D equivalent crack opening e strain g viscosity parameter h temperature j history variable l shear modulus m Poisson’s ratio n hardening rate parameter p ratio of the circumference of a circle to its diameter r stress-like parameter in the cohesive model r stress tensor ry yield strength rVM von Mises equivalent stress Subscripts/superscripts 0 related to the initial state of a process BC related to the BC

2141

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

d e g iso n vp s S th trans TGO TC u

related to the damage process in the cohesive zone related to elastic deformation related to chemical growth related to isotropic behaviour normal component related to viscoplastic deformation shear (tangential) component related to the substrate related to thermal deformation related to transversely isotropic behaviour related to the TGO related to the TC related to the ultimate state of a process

cooling are induced because the thermal expansion coefficient of the TGO is lower than that of the BC, and the interface between the two components is coherent in the elastic regime. Nevertheless, local tensile stresses can arise at geometrical imperfections of the TGO layer. These tensile stresses may cause the development of micro-cracks, which, after coalescence, ultimately can lead to spallation failure of the TC [5,6]. TGO growth and damage evolution in TBCs has been studied extensively over the past decade [7–12], where the importance of the TGO morphology on the lifetime of TBC systems is clearly pointed out. However, modeling approaches have been limited to predefined cracking paths and a prescribed evolution of the TGO geometry. In the present contribution, the TGO growth and damage evolution are modeled without these restrictions. To this end, the numerical implementation of the models is carried out using a mesh-independent formulation. In Section 2, the growth of the TGO layer is modeled using a diffusion–reaction model that connects the TGO volume fraction to the concentration of free oxygen. Inherent to this approach is the appearance of a BC–TGO mixture zone. The constitutive models of the layer materials and the mixture zone are presented in Section 3; this section also describes the cohesive zone model used for the simulation of discrete cracks. The material parameters of the models are discussed in Section 4. Section 5 presents a numerical analysis of a representative TBC system subjected to a thermal cycling process. A parametric study is performed for different fracture strengths of the TC to determine its influence on the durability of a TBC system. In Section 6, the main conclusions of the present study are summarized and validated against experimental observations reported in the literature.

2. Growth model for the TGO layer The TGO layer is composed of alumina (a-Al2O3) that results from the oxidation of aluminum present in the BC, as shown schematically in Fig. 1. The oxygen diffuses through the TC towards the BC, while the aluminum is contained in the BC. In the TGO layer, the diffusivity of oxygen is orders of magnitude larger than the diffusivity of aluminum [13,14]. Hence, the development of the TGO occurs primarily at the actual TGO–BC interface [15]. Moreover, it is assumed that the oxidation x2

c=c

_

hTC

~ ~

hTGO

TC

Oxygen diffusion

TGO TGO-BC mixture

Aluminum in BC hBC

~ ~

n=1

2Al + 3O → Al2O3

0
BC

n=0 c .n = 0

Time Dwelling period

Δ

Fig. 1. Schematization of the growth process of the TGO layer during the dwell phase of a thermal cycling process. The cross-section left shows the TBC components (TC, TGO and BC) across the height at the beginning of the dwell phase. The increase in the TGO volume fraction n by oxidation leads to an increasing thickness hTGO of the region containing pure alumina.

2142

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

reaction takes place only during the so-called dwell phase, which is the high-temperature plateau of a typical thermal cycle [2,16]. 2.1. Oxygen diffusion–reaction model The formation of alumina is assumed to occur gradually in regions previously occupied by the BC material. During the oxidation process, zones of newly formed Al2O3 are mixed with the BC material until eventually the TGO material fully occupies the region. Correspondingly, the growth of the TGO layer is described by the volume fractions of the TGO material, nTGO, and the BC material, nBC. These volume fractions take values between 0 and 1 and satisfy the necessary constraint nBC + nTGO = 1. For convenience, the TGO volume fraction is denoted as n: = nTGO, hence nBC = 1  n. During the oxidation process the volume fraction varies from n = 0 (pure BC material) to n = 1 (pure TGO material). The formation of the TGO layer is governed by the availability of the reactants, i.e., the free (mobile) oxygen and aluminum. The distribution of concentration of free oxygen c is governed by the following diffusion–reaction equation:

@c  r  ðDrcÞ ¼ S; @t

ð1Þ

where D is the oxygen diffusion coefficient and S is a sink term representing the spatial confinement of mobile oxygen due to TGO formation, i.e., the mobile oxygen does not further diffuse when it reacts with aluminum. A similar equation governs the concentration of aluminum available for oxidation. However, since the diffusion of this aluminum is negligible, its time rate of change is directly proportional to the oxidation rate, which in turn is proportional to the term S in (1). An expression for the sink term S follows from a reaction model. Here, it is assumed that the rate of TGO formation (i.e., @n/@t) is proportional to the product of the concentrations of mobile oxygen and aluminum. In turn, the concentration of Al is assumed to scale with the volume fraction 1  n of the BC (where the available aluminum is concentrated). Correspondingly, the rate of TGO formation can be expressed as

@n ¼ Rð1  nÞc; @t

ð2Þ

where R is a constant of proportionality that specifies the reaction rate. The rate of formation of TGO can be related to the sink term S in (1) as follows:

S¼M

@n ; @t

ð3Þ

where M is the ratio of confined oxygen atoms to formed TGO volume, which is equal to the molarity of oxygen in alumina. In view of (1)–(3), the concentration of mobile oxygen c and the volume fraction of TGO n can be determined by solving the following (coupled) system of equations:

@c  r  ðDrcÞ ¼ MRð1  nÞc; @t @n ¼ Rð1  nÞc: @t

ð4Þ

The boundary conditions for the diffusion problem are reflected by the oxygen concentration c and normal oxygen flux Drc  n prescribed at the top and bottom of the TBC system, respectively, with n the normal outward unit vector (see also Fig. 1). The oxygen diffusion coefficient of the TC is typically much larger than that of the TGO (DTC  DTGO) and it may assumed that the oxygen concentration at the TC–TGO interface is approximately equal to the prescribed value c at the top of the TBC. In addition, the normal oxygen flux at the bottom may be assumed to be negligible (i.e., rc  n = 0). The evolving boundaries that separate the TGO–BC mixture zone from the (pure) TGO and the (pure) BC layers can be determined from the conditions n = 1 and n = 0 respectively (see Fig. 1). TGO For the first thermal cycle, the thickness of the TGO is given as h0 , which reflects the thickness due to a pre-oxidation performed prior to deposition of the TC layer, and the concentration of oxygen is taken as zero below the TGO–BC interface (i.e., the mixture zone has zero thickness initially). 2.2. Growth strain The formation of the TGO is typically associated with the volumetric expansion of the material [10,17–19,5,20–23,13]. In the BC–TGO mixture zone, where the additional TGO is formed, the volumetric expansion is accounted for through a growth strain tensor eg of the form

eg ¼ neg I; g

ð5Þ

where e is the growth strain associated to the newly formed TGO and I is the second-order identity tensor. In general, the total strain is expressed as e = ee + evp + eth + eg, where ee,evp and eth are the elastic, viscoplastic and thermal strains, respectively. In view of this strain decomposition, the growth strain provides an important coupling between the oxidation process and the associated state of stress. The constitutive models describing the elastic and viscoplastic behaviors of the TBC components can be found in the subsequent section.

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

2143

3. Constitutive models This section contains the mechanical constitutive models, a fracture model and a description for enhanced diffusion due to cracking in the materials present in a TBC system. A rule of mixtures based on Voigt’s assumption is used to derive the models that describe the material response in the mixture zone. 3.1. Constitutive models for the TBC components The TC material is a brittle ceramic whose mechanical response prior to fracture can be simulated using an elastic model. Due to the EB–PVD deposition process, the microstructure of the TC is not homogeneous [18]. Close to the TGO interface, the TC consists of small grains that are randomly oriented. However, further away from the interface, the morphology of the TC is characterized by larger grains forming a columnar structure that is oriented perpendicular to the interface and extends across the thickness of the TC. Consequently, the behavior of the TC can be described using an isotropic model in the region close to the TGO interface whereas a transversely isotropic model is more appropriate for the region further away from the interface. The transition from the isotropic to the transversely isotropic elastic properties can be captured by means of a position-dependent interpolation function for the effective elastic stiffness tensor CTC of the TC, i.e.,

with

b TC ðXÞ ¼ ð1  XÞCTC þ XCTC CTC ¼ C iso trans

ð6Þ

8 0 > < x x 2 iso b X ¼ X ðx2 Þ ¼ x  xiso > trans : 1

ð7Þ

if

x2 6 xiso ;

if

xiso < x2 < xtrans ;

if

x2 P xtrans :

TC In (6), CTC iso and Ctrans represent the isotropic and transversely isotropic elastic stiffness tensors, respectively. As indicated in (7), the interpolation variable X depends linearly on the vertical coordinate x2 shown in Fig. 5. In the region x2 6 xiso the mechanical response is fully isotropic whereas in the region x2 P xtrans it is fully transversely isotropic. The tensor CTC iso is characterized by a Young’s modulus E and a Poisson’s ratio m. The tensor CTC trans is given by the transversely isotropic material properties E1, E2, l2, m12, and m13, where the subindices 1, 2 and 3 refer to a Cartesian coordinate system such that, as shown in Fig. 5, the x2-axis coincides with the direction perpendicular to the plane of isotropy. The compliance tensor DTC trans is the inverse of the transversely isotropic elastic tensor CTC trans and is defined as

2

 TC 1 DTC trans ¼ Ctrans

1

6 6 m12 6 6 m 13 16 6 ¼ 6 0 E1 6 6 6 0 6 4 0

3

m12 E1 E2 m12

m13

0

0

0

m12

0

0

1

0

0

0

0 E1

7 0 7 7 0 7 7 7 ; 0 7 7 7 0 7 7 E1 5

0

0

0

0

l2

0

0

2ð1 þ m13 Þ

0

0

ð8Þ

l2

where Voigt’s notation has been used such that the stress and strain tensor components ‘‘11, 22, 33, 23, 13, and 12’’ correspond to the matrix rows (and columns) 1–6, respectively. The TGO is modeled as isotropically linear–elastic. The BC is modeled as isotropically elasto-viscoplastic, with the viscoplastic deformation behavior characterized by the Perzyna model, see e.g., [24,25]. The Perzyna viscoplastic model accounts for rate-effects that may occur at elevated temperatures, such as plastic creep and relaxation. The yield strength ry thereby evolves in agreement with an exponentially-saturating hardening law, i.e.,





ry ¼ r^ y ðjvp Þ ¼ ry0 þ ryu  ry0 ð1  exp ðnjvp ÞÞ; y 0

ð9Þ

y u

where r is the initial yield strength, r is the ultimate yield strength, and n is a parameter that sets the hardening rate. The viscoplastic history variable jvp appearing in (9) evolves according to

j_ vp ¼

hf p i

g

with f p ¼

rVM  ry ; ry

ð10Þ

where j_ vp represents the time derivative of jvp, rVM is the von Mises stress, and fp is the rate-independent (plastic) yield surface normalized by the actual yield strength ry. The Macaulay brackets hxi ¼ 12 ðx þ jxjÞ ensure that only positive values of x are taken into account. Consequently, viscoplastic deformations only develop when overstress is generated, as reflected by positive values of the rate-independent yield surface fp, where the parameter g quantifies the rate-dependency of the response. Finally, in the case of von Mises plasticity the viscoplastic strain rate is computed in the standard way, where the flow direction is determined from the stress derivative of the rate-independent yield surface, i.e.,

e_ vp ¼ j_ vp

@f p ; @ rBC

with rBC the stress in the bond coat.

ð11Þ

2144

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

The effective constitutive behavior in the BC–TGO mixture zone is modeled according to Voigt’s assumption. Hence, the total strain e in a material point of the mixture zone is assumed to be the same for both constituents, i.e.,

e ¼ eBC ¼ eTGO :

ð12Þ

Correspondingly, the effective stress state r in the mixture region is obtained as the volume average of the stresses of the constituents, i.e.,

r ¼ ð1  nÞrBC þ nrTGO ; where r

BC

and r

TGO

ð13Þ

represent the stress tensors in the BC and TGO phases, respectively.

3.2. Cohesive zone model A cohesive zone model, which specifies the traction–separation relationship across a developing crack, is used to simulate the fracture process under thermo-mechanical loading. When the fracture strength of a material is reached at some location, a surface of displacement discontinuity is introduced. The traction across this surface of discontinuity, as described by the cohesive zone model, gradually reduces from the fracture strength to zero as the displacement jump across the surface increases. Correspondingly, the surface of discontinuity is interpreted as a crack. The cohesive zone model for a single phase material is described in detail in [26], and is extended in the present study to describe the failure response in the BC–TGO mixture zone. For completeness, the cohesive zone model is briefly summarized below. The traction t = (tn, ts) acting on the crack faces is related to the crack opening d = (dn, ds), where the subscripts ‘‘n’’ and ‘‘s’’ refer to the directions normal and tangential to the crack face, respectively. An equivalent crack opening D is defined as

D¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hdn i2 þ b2 d2s ;

ð14Þ

where b is a weighting factor for the mode I and mode II contributions. The equivalent crack opening is used to compute the equivalent traction T as

T ¼ Tb ðD; jd Þ ¼

8 < g^ðDÞ : g^ðjd Þ

if

D

jd

f d ¼ 0 and j_d > 0;

otherwise;

ð15Þ

where jd is a damage history variable, fd is a loading function, and g is the effective traction–separation law. The upper and lower expressions in (15) provide the equivalent traction during, respectively, crack growth and unloading/reloading. The specific form of the effective traction–separation law used in (15) is 3

g ¼ g^ðDÞ ¼

rd2 Dc 3

ðd þ DÞ2þc

ð16Þ

:

b equals the fracIn the above expression, the model parameters r and d are chosen such that the maximum of the function T b equals the material fracture toughness Gt, see [26] for more details on the calibration ture strength ft and the integral of T procedure. Furthermore, the parameter c determines the overall shape of the traction–separation law. Using a loading function

f d ¼ ^f d ðD; jd Þ ¼ D  jd ;

ð17Þ

the Karush–Kuhn–Tucker relations reflecting the loading and unloading conditions are as follows:

f d j_ d ¼ 0

f d 6 0;

j_ d P 0:

ð18Þ

After evaluating (15), the normal and shear tractions can be computed as

8 < dn T if tn ¼ D : Kdn if 2 ds T: ts ¼ b D

dn > 0; dn < 0;

ð19Þ

Special care needs to be taken with respect to the continuity of tractions at the onset of cracking. In particular, the traction t = rn computed prior to the introduction of a surface of displacement discontinuity with normal n needs to be the same as the traction obtained from the traction–separation law at zero actual crack opening. For the model given by (14)–(19), the traction continuity can be satisfied upon using a shift d0 in the crack opening vector, see [26] for more details. In the BC–TGO mixture zone, in view of the different stresses rBC and rTGO in the BC and TGO constituents as reflected by (13), this vector needs to be computed separately for each constituent, leading to dBC and dTGO . In addition, the 0 0 traction–separation response in the mixture zone is computed applying Voigt’s assumption to the BC and TGO constituents,

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

2145

Fig. 2. The cohesive law of (1) a BC–TGO mixture (illustrated here for a TGO volume fraction of n = 0.5) is the volume–weighted average of the cohesive laws of (2) the pure BC constituent and (3) the pure TGO constituent.

i.e., the actual crack opening  d is considered to be the same for both constituents. Accordingly, the crack openings for the BC and TGO constituents, as used in (14), are

di ¼ di0 þ d with i ¼ BC; TGO:

ð20Þ

In correspondence with (20), the description of fracture in the BC and TGO phases requires separate history variables jd,BC,jd,TGO that satisfy (18) individually. Furthermore, the cohesive law is evaluated separately for each phase, yielding the individual tractions t BC ¼ ^t BC ðdBC ; jd;BC Þ and t TGO ¼ ^t TGO ðdTGO ; jd;TGO Þ. The effective cohesive traction t follows as the volume average of the tractions of the constituents, i.e.,

t ¼ ð1  nÞt BC þ nt TGO :

ð21Þ

Fig. 2 illustrates the equivalent traction–separation (TD) law in the BC–TGO mixture zone, together with the corresponding cohesive laws in the BC and TGO constituents, for a mixture with equal volume fractions (n = 0.5). Note that the TGO is assumed to be more brittle and less strong than the BC. Upon oxidation, the TGO volume fraction n increases, as a result of which the effective traction that can be transferred across the crack in the mixture zone reduces. Consequently, oxidation may force the crack to open further and thereby reduces the stress level in the surrounding bulk material. 3.3. Oxygen diffusion coefficient Increased oxidation along cracks was observed in the experimental study reported in [27]. In order to account for this behavior, the oxygen diffusion coefficient D that appears in (1) is augmented along crack trajectories as follows:

D ¼ D



m Gt Gt  E

ð22Þ

where D⁄ is the diffusion coefficient of the undamaged bulk material, Gt is the fracture toughness, m is a calibration parameter and E is the energy dissipated during crack opening, which can be computed as

Eðjd Þ ¼

Z 0

jd

1 Tb ðD; DÞ dD  jd Tb ðjd ; jd Þ: 2

ð23Þ

Note that the diffusivity becomes infinitely large when the crack has completely opened, i.e., D ? 1 as jd ? 1. In the mixture zone, the effective diffusion coefficient is determined in accordance with a volume average of the diffusion coefficients of the constituents, i.e.,

D ¼ ð1  nÞDBC þ nDTGO :

ð24Þ

4. Model parameters 4.1. Material parameters for diffusion–reaction model The values of the parameters for the TGO growth model are given in Table 1. The growth of the TGO layer is determined by the diffusion-reaction model given by (4). The parameter R appearing in this model is numerically calibrated from the experimentally observed evolution of the TGO thickness reported in previous work [28], see Fig. 3. The TBC system used

2146

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152 Table 1 Material parameters of the growth model. Parameter

Value

Equation

TC diffusivity (for oxygen) TGO diffusivity BC diffusivity Reaction rate parameter

DTC ? 1 DTGO = 2.0  1014 m2 s1 DBC = 2.0  1014 m2 s1

(1) (1) (1) (3)

Oxygen molarity in Al2O3 Growth strain

A

1

R ¼ 104 m3 mol s1 M = 1.11  105 mol m3 eg = 3.0  103

(3) (5)

B

_ Fig. 3. (A) Spatial profiles of the oxygen concentration c, the TGO volume fraction n, and the rate of change of the TGO volume fraction nð¼ @n=@tÞ after 1000 cycles, as obtained from the growth model for the TGO layer. The thickness of the mixture zone is approximately 1 lm. (B) The evolution of the TGO thickness hTGO, obtained with the calibrated diffusion–reaction model (solid line), compared to experimental measurements (crosses with error bars).

in the calibration procedure is assumed to have straight interfaces between the constituents. Because the diffusivity of the TC is orders of magnitude larger than that of the TGO and the BC [29], and the oxidation reaction does not take place in the TC, the calibration may be performed on a reduced TBC system where the TC is omitted. Here, the oxygen diffusion coefficients of the TGO and the BC are assumed to be equal (DTGO  DBC), and are adopted from [13]. The TGO and the BC have initial thicknesses of 1 lm and 139 lm, respectively, in accordance with the values measured for the TBC test samples studied in [28]. The TBC specimen is subjected to 1000 thermal cycles. Each thermal cycle lasts 1 h, and starts with a 45 min dwell phase at 1095 °C, followed by a cooling and reheating phase where the temperature is brought down at constant thermal rate of 2.2 °C/s to 105 °C and back to 1095 °C. The oxygen concentration c at the top of the TBC system (=top of the TGO layer) is assumed to be in correspondence with a partial oxygen pressure of pO2 ¼ 105 Pa (i.e., equal to the atmospheric pressure) at a dwell temperature of 1095 °C. The value of the oxygen concentration at this specific temperature could not be found in the literature, and therefore has been derived from a relation between the solubility of hydrogen in alumina and the temperature, as experimentally established by Roy and Coble [30]. Based on this relation and the working assumption that the solubility is inversely proportional to the atomic volume of the solute, the oxygen concentration at the top of the TGO layer could be estimated as c ¼ 1:55 mol m3 . In addition, at the bottom of the TBC system (=bottom of the BC layer) the oxygen flux is prescribed to be zero. 4.2. Material parameters for mechanical models As discussed in [5], the volumetric increase of the TGO layer may be assumed to result from (i) the addition of material (oxygen) and (ii) the development of a growth strain. This assumption is motivated from the fact that the oxidation process takes place in an open system, with a net influx of oxygen occurring over the upper boundary of the TBC system. Accordingly, the value of the growth strain is substantially smaller than in studies where it is supposed to fully represent the volumetric increase following from the oxidation process [13,17,18]. The actual value of the growth strain is estimated as eg = 3.0  103, which is of the same order of magnitude as the value reported in [21]. From preliminary simulations not reported here, the in-plane compressive stresses induced solely by the growth of the oxide layer then do not surpass 1 GPa, which is in accordance with the stress values reported in [5,22,23]. For the TC layer, the spatial transition from an isotropic model to a transversely isotropic model, as reflected by (7), is determined by the values xiso = 135 lm and xtrans = 200 lm, which are estimated from typical TC microstructures. As can be observed from Fig. 5, the origin from which the above values are measured is located at the bottom of the bond coating. This means that at the valley of the TGO irregularity, which is maximal along the axis of symmetry at a vertical distance of 127 lm from the origin, the TC starts as fully isotropic over a relatively small vertical distance of 135  127 = 8 lm. Above the isotropic region, the TC gradually changes from isotropic to transversely isotropic over a vertical distance of 200  135 = 65 lm, and beyond this transition zone becomes fully transversely isotropic.

2147

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152 Table 2 Thermo-mechanical bulk properties of the components of the TBC system. Parameter

Value

Equation

TC layer Coeff. of therm. expansion Young’s modulus Poisson’s ratio Young’s modulus Young’s modulus Shear modulus Poisson’s ratio Poisson’s ratio

a = 12  106 K1 E = 60 GPa m = 0.15 E1 = 68 GPa E2 = 170 GPa l2 = 77 GPa m12 = 0.04 m13 = 0.1

(8) (8) (8) (8) (8)

TGO layer Coeff. of therm. expansion Young’s modulus Poisson’s ratio

a = 8.8  106 K1 E = 393 GPa m = 0.25

BC layer Coeff. of therm. expansion Young’s modulus Poisson’s ratio Initial yield strength

a = 14  106 K1 E = 200 GPa m = 0.3 ry0 ¼ 435 MPa

ryu ¼ 565 MPa

(9)

Ultimate yield strength Hardening rate parameter Viscosity parameter

n = 77 (–) g = 2000 s

(9) (9) (10)

Substrate Coeff. of therm. expansion

a = 15  106 K1

(25)

Table 3 The fracture parameters of the components of a TBC system. The system response has been assessed with respect to three different values for the TC fracture strength. Parameter

TC layer

TGO layer

BC layer

Equation

Fracture strength, ft (MPa) Fracture toughness, Gt (N mm1) Weighting factor, b (–) Shape parameter, c (–) Elastic stiffness, K (N m3) Damage exponent for diffusivity, m (–)

200, 300, 600 0.01 1.4 0.35 1014 6

285 0.02 1.4 0.35 1014 6

500 0.3 1.4 0.35 1014 6

(7) in Ref. [26] (7) in Ref. [26] (14) (15) (19) (22)

The constitutive parameters of the various models are summarized in Tables 2 and 3. The values of the parameters are equal to those reported in [26], except for the viscous parameter g of the viscoplastic model for the BC and the various values of the TC fracture strength ftTC . The parameter g has been calibrated from experimental data for a Ni-based superalloy reported in [31], since the literature lacks such data for the BC material. The literature provides a wide range of values for the fracture strength of yttria-stabilized zirconia—the material that constitutes the TC. Depending on the measuring technique and the type of TC microstructure, values of 140 MPa and higher are reported in [32]. Other references report values in the range of 200–400 MPa [33–35]. Further, relatively high fracture strengths, ranging from 450 to 750 MPa, can be found in [36–38]. Due to this spread in reported values, a sensitivity analysis of the TBC system with respect to the TC fracture strength is carried out by considering three different values, namely 200, 300, and 600 MPa. These values will be referred to as the low, medium, and high strength cases, respectively.

5. Thermal cycling simulations The response of a TBC system exposed to 1000 thermal cycles is now analyzed, and its lifetime sensitivity with respect to the fracture strength of the TC is assessed. To this end, the models described in Sections 2 and 3 are applied to a representative domain that is subjected to the cyclic thermal loading sketched in Fig. 4. The loading conditions are representative of those applied in the experiments reported in [28]. The ‘‘partition-of-unity method’’ (also known as the ‘‘generalized finite element method’’ [39], or as the ‘‘extended finite element method’’ [40]) is employed for the kinematical description of displacement discontinuities—i.e., cracks—within the finite elements, see also [41,42]. The elements intersected by a crack are adaptively enhanced with additional degrees of freedom that are connected to a discontinuous basis function. Accordingly, the displacement jump across the crack is

2148

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

Fig. 4. Cyclic thermal loading imposed on the TBC system during the lifetime simulation. Each thermal cycle consists of a 45 min dwell phase, and a 15 min cooling and heating phase from a maximum temperature of 1095 °C to a minimum temperature of 105 °C and back to the maximum temperature.

A

B

Fig. 5. Model of the TBC system. Values of the geometrical parameters depicted in A are given in the table at the left-hand side. The mechanical boundary conditions are visualized in B. Due to symmetry only half of the domain is modeled.

accounted for using the enhanced degrees of freedom. This method allows cracks to be inserted during the simulation at the required locations in the domain and in the required orientations according to the fracture criterion. 5.1. Geometry and boundary conditions The sample modeled is representative of a TBC cross-section with a geometrical irregularity at the BC–TC interface, which was used in [28] to experimentally study the effect of damage evolution at interfacial irregularities (see Fig. 4). The BC and the TC have initial nominal thicknesses of 139 lm and 125 lm, respectively. Locally, these values vary due to the presence of the irregularity, as characterized by the geometrical parameters d1 and d2 shown in Fig. 5. The initial TGO thickness of 1 lm is in line with a pre-oxidation treatment applied to the TBC system before depositing the TC [11]. The domain has been discretized by nearly 4000 six-noded, isoparametric triangular finite elements that are equipped with a three-point Gauss quadrature. The displacement boundary conditions make use of the fact that the thickness of the substrate and the out-of-plane extension of the domain are significantly larger than the thickness of the TBC system. Thus, it is assumed that the thermal deformation of the substrate is not influenced by the presence of the TBC system. Correspondingly, the uniform thermal deformation at the substrate’s upper surface is straightforwardly imposed on the TBC system through displacement boundary conditions. Assuming an initial stress-free reference state at the (maximum) temperature h0 = 1095 °C, the horizontal displacement u1 at the substrate–BC interface follows as

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

^ 1 ðx1 ; x2 ¼ 0Þ ¼ as ðh0  hÞðL  x1 Þ; u1 ¼ u

2149

ð25Þ

where as is the coefficient of thermal expansion of the substrate and h is the current temperature. For consistency, the horizontal displacements at the two bottom corners are also imposed along the vertical boundaries of the domain, i.e.,

^ 1 ðx1 ¼ 0; x2 Þ ¼ as ðh0  hÞL and u1 ¼ u ^ 1 ðx1 ¼ L; x2 Þ ¼ 0: u1 ¼ u

ð26Þ

The tangential traction is assumed to be zero at the vertical boundaries x1 = 0 and x1 = 2L and at the vertical centerline x1 = L. Hence, the problem is considered to be symmetric with respect to the vertical centerline. In addition, the upper surface is taken as traction-free and the rigid support of the TBC by the substrate is accounted for through the condition

^ 2 ðx1 ; x2 ¼ 0Þ ¼ 0: u2 ¼ u

ð27Þ

In the direction perpendicular to the domain modeled (i.e., in x3-direction), the substrate’s deformation is imposed on the TBC system by means of a generalized plane-strain condition, i.e.,

e33 ¼ as ðh0  hÞ:

ð28Þ

The evolution of crack patterns will be studied first for a TC layer with a medium fracture strength, which is named the ‘‘reference simulation’’. Subsequently, the results of this simulation are compared to those obtained using the low and high fracture strengths reported in Table 3. 5.2. Reference simulation with medium fracture strength of top coat Fig. 6 shows the evolution of the cracks and the growth of the TGO for an increasing number of thermal cycles in a TBC system with a medium TC fracture strength of 300 MPa. The thickness of the TGO develops relatively fast during the earlier stage of the loading process, which is in accordance with a commonly observed parabolic growth rate [5]. Inherent to a diffusion process, the initial irregularity at the BC–TGO interface flattens somewhat during the growth of the TGO. During the first cycles several small cracks (represented by surfaces of discontinuity) nucleate at the BC–TGO interface at locations ranging from the crest to the valley of the irregularity. Upon thermal cycling, the BC–TGO interface progresses into the BC, and the initial cracks become fully embedded in the TGO layer. From cycle 160 to cycle 260, more cracks nucleate at the current position of the BC–TGO interface. The TGO increases in thickness and provides sufficient strain energy to propagate these cracks over a short distance along the interface. As the thermal cycling further continues, more cracks emerge, where each crack propagates as long as its location roughly coincides with the BC–TGO interface. Once the interface has moved sufficiently away from these cracks, they become arrested. However, cracks that nucleate at the BC–TGO interface after about 400 cycles remain propagating over a significant distance until the end of the cyclic loading process is reached, which occurs at 1000 thermal cycles.

Fig. 6. The evolution of crack patterns (thick lines) and the TGO interfaces (thin lines), as observed after 250, 500, 750, and 1000 thermal cycles for a TBC system with a medium fracture strength of 300 MPa for the TC.

2150

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

5.3. Comparison study of the fracture strength of the top coat To assess the influence of the TC fracture strength, the reference case with medium TC strength, ftTC ¼ 300 MPa, is compared to the low and high strength cases, ftTC ¼ 200 MPa and ftTC ¼ 600 MPa, respectively. Fig. 7 shows the crack pattern and the TGO geometry after 1000 thermal cycles for these three cases. In the simulation with the highest TC fracture strength, crack development in the TC is completely prevented, and this case does not differ significantly from the case with the medium fracture strength. The opposite holds for the simulation with the low value for the TC fracture strength. After 790 cycles, a large crack propagates through the TC, spanning the irregularity by coalescing with its mirror crack at the centerline of the sample. This induces a redistribution of stress, which forces the cracks located at the crest of the irregularity to propagate further in the direction parallel to the BC–TGO interface. Within 10 more cycles, large cracks form at the left side of the crest in the TC. This is again followed by crack propagation in the BC–TGO interface region. At this stage, the TC loses its strength in the vertical direction, which may eventually lead to spallation of (parts of) the TC. Observe that, as a consequence of the redistribution of stress due to cracking, vertical cracks nucleate in the BC layer, followed by a local bulging of the BC– TGO interface at the crest. This bulging can be ascribed to an increase in oxygen diffusivity along the crack paths, known as crack-enhanced diffusion. Fig. 8 shows the cyclic evolution of the vertical normal stress r22 in material point P (indicated in Fig. 7) at the minimum temperature during a thermal cycle, i.e., 105 °C. For the cases where the TC has a medium and a high strength, the stress response is virtually the same and increases monotonically with the number of thermal cycles. The monotonic increase in stress can be ascribed to the (approximately parabolic) growth of the TGO layer; since the TGO layer has a mismatch in elastic and thermal properties with the TC and BC layers, during its growth it continuously changes the stress pattern in its direct neighborhood.

Fig. 7. Crack pattern (thick lines) and TGO configuration (thin lines) as observed after 1000 thermal cycles in the three TBC systems characterized by a low (200 MPa), medium (300 MPa), and high (600 MPa) TC fracture strength ftTC .

Fig. 8. The evolution of the vertical normal stress r22 in a TC material point P shown in Fig. 7. The responses are evaluated at the minimum temperature in the cyclic loading process of 105 °C. The stress response is plotted for three different TC fracture strengths, ftTC ¼ 200, 300, and 600 MPa.

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

2151

Fig. 9. The evolution of the total crack length in three TBC systems characterized by a different value of the TC fracture strengths, ftTC ¼ 200, 300, and 600 MPa.

For the system with the low TC fracture strength, the stress response is qualitatively similar to that for the cases with the medium and high strength values, but the effect of cracking is more pronounced. Indeed, at 450 cycles and at 700 cycles sudden increases in stress occur. These stress increases are due to stress redistributions induced by cracks nucleating close to the crest of the irregularity and propagating towards the centerline of the TBC specimen. At 790 cycles, a large crack starts to develops across the entire width of the TGO irregularity, as a result of which the vertical normal stress in reference point P drops abruptly. Fig. 9 visualizes the evolution of the total crack length for the three investigated fracture strengths of the TC. For the case of a low TC fracture strength of 200 MPa, the total length of all cracks accumulates to 376 lm at 1000 cycles, which is about six times larger than for the other two cases. The figure further illustrates that for the case with the low TC fracture strength, at 450, 700, and 790 cycles the cracks propagate noticeably. These thermal cycles correspond to the above-mentioned values at which the stress locally changes abruptly. As already observed in Figs. 7 and 8, the cases with a TC fracture strength of 300 MPa and 600 MPa exhibit virtually the same response. 6. Conclusions The present simulations show that cracks may develop in TBC systems already at early stages of a cyclic thermal loading process. Cracks nucleate primarily at the actual position of the BC–TGO interface, due to the mismatch in elastic and thermal properties. As the TGO grows and this interface translates, an array of virtually parallel cracks emerges. However, as long as their development remains limited to the neighborhood of an irregularity, these cracks are not detrimental to the coating’s lifetime. Furthermore, catastrophic failure requires these cracks to also develop across the TC layer. A comparison study of the TC fracture strength ftTC has illustrated that crack development in the TC is about to happen at TC ft ¼ 300 MPa, see Figs. 6 and 7. Hence, for TBC systems comparable to the system studied in this communication, it may be concluded that cracking in the TC can be avoided by manufacturing TCs with a fracture strength higher than 300 MPa. The geometry analyzed in the present simulations has also been assessed experimentally, as reported in [28]. Since cracking of the TC was not observed in the experiments, it may be concluded that the TC fracture strength in the experimental samples is (well) above 300 MPa, since the present simulations show that for ftTC > 300 MPa no significant cracking occurs in the TC. Furthermore, for the medium and high TC strength cases, the size, orientation and location of the TGO cracks corresponds well with those observed experimentally, see [28]. This indicates that the present modeling approach is appropriate for realistically simulating the complex mechanisms of oxide growth and damage evolution in TBC systems subjected to thermal cycling. Acknowledgments This research was carried out under the Project Number MC7.04186 in the framework of the Research Program of the Materials innovation institute M2i (www.m2i.nl), the former Netherlands Institute for Metals Research. The authors appreciate the helpful discussions with Dr. Erik-Jan Lingen (Habanera numerical software) on numerical implementation issues within the finite element framework JEM/JIVE. In addition, the authors acknowledge the useful discussions with Dr. Wim Sloof (Delft University of Technology) and Dr. Thijs Nijdam (National Aerospace Laboratories of the Netherlands, NLR) on the thermo-mechanical behavior of TBC systems during the preparation of this article. References [1] Kear BH, Thompson ER. Aircraft gas turbine materials and processes. Science 1980;208:847–56. [2] Nicholls JR. Advances in coating design for high-performance gas turbines. MRS Bull 2003;28:659–70.

2152

T.S. Hille et al. / Engineering Fracture Mechanics 78 (2011) 2139–2152

[3] Tzimas E, Müllejans H, Peteves SD, Bressers J, Stamm W. Failure of thermal barrier coating systems under cyclic thermomechanical loading. Acta Mater 2000;48:4699–707. [4] Pindera M-J, Aboudi J, Arnold SM. The effect of interface roughness and oxide film thickness on the inelastic response of thermal barrier coatings to thermal cycling. Mater Sci Engng A 2000;284:158–75. [5] 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. [6] Chen X, Hutchinson JW, He MY, Evans AG. On the propagation and coalescence of delamination cracks in compressed coatings: with application to thermal barrier systems. Acta Mater 2003;51:2017–30. [7] Mumm DR, Evans AG. On the role of imperfections in the failure of a thermal barrier coating made by electron beam deposition. Acta Mater 2000;48:1815–27. [8] Smialek JL. A deterministic interfacial cyclic oxidation. Acta Mater 2003;51:469–83. [9] Caliez M, Chaboche J-L, Feyel F, Kruch S. Numerical simulation of EBPVD thermal barrier coatings spallation. Acta Mater 2003;51:1133–41. [10] Aktaa J, Sfar K, Munz D. Assessment of TBC systems failure mechanisms using a fracture mechanics approach. Acta Mater 2005;53:4399–413. [11] Nijdam TJ, Marijnissen GH, Vergeldt E, Kloosterman AB, Sloof WG. Development of a pre-oxidation treatment to improve the adhesion between thermal barrier coatings and NiCoCrAlY bond coatings. Oxid Met 2006;66:269–94. [12] Guo SQ, Tanaka Y, Kagawa Y. Effect of interface roughness and coating thickness on interfacial shear mechanical properties of EB–PVD yttria-partially stabilized zirconia thermal barrier coating systems. J Eur Ceram Soc 2007;27:3425–31. [13] Rösler J, Bäker M, Volgmann M. Stress state and failure mechanisms of thermal barrier coatings: role of creep in thermally grown oxide. Acta Mater 2001;49:3659–70. [14] Fox AC, Clyne TW. Oxygen transport by gas permeation through the zirconia layer in plasma sprayed thermal barrier coatings. Surf Coat Technol 2004;184:311–21. [15] Karlsson AM, Hutchinson JW, Evans AG. The displacement of the thermally grown oxide in thermal barrier systems upon thermal cycling. Mater Sci Engng A 2003;351:244–57. [16] Zhu HX, Fleck NA, Cocks ACF, Evans AG. Numerical simulations of crack formation from pegs in thermal barrier systems with NiCoCrAlY bond coats. Mater Sci Engng A 2005;404:26–32. [17] Busso EP, Lin J, Sakurai S, Nakayama M. A mechanistic study of the oxidation-induced degradation in a plasma-sprayed thermal barrier coating system. Part I: Model formulation. Acta Mater 2001;49:1515–28. [18] Busso EP, Qian ZQ. A mechanistic study of microcracking in transversely isotropic ceramic–metal systems. Acta Mater 2006;54:325–38. [19] Busso EP, Wright L, Evans HE, McCartney LN, Saunders SRJ, Osgerby S, et al. A physics-based life prediction methodology for thermal barrier coating systems. Acta Mater 2007;55:1491–503. [20] Evans HE. Stress effects in high-temperature oxidation of metals. Int Mater Rev 1995;40:1–40. [21] He MY, Hutchinson JW, Evans AG. Simulation of stresses and delamination in a plasma-sprayed thermal barrier system upon thermal cycling. Mater Sci Engng A 2003;345:172–8. [22] Karlsson AM, Evans AG. A numerical model for the cyclic instability of thermally grown oxides in thermal barrier systems. Acta Mater 2001;49:1793–804. [23] Padture NP, Gell M, Jordan EH. Materials science – thermal barrier coatings for gas-turbine engine applications. Science 2002;296:280–4. [24] Heeres OM, Suiker ASJ, de Borst R. A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model. Eur J Mech A – Solids 2002;21:1–12. [25] Armero F. Elastoplastic and viscoplastic deformations in solids and structures. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics.: Solids and structures, vol. II. John Wiley & Sons; 2004. p. 227–66 [Chapter 7]. [26] Hille TS, Suiker ASJ, Turteltaub S. Microcrack nucleation in thermal barrier coating systems. Engng Fract Mech 2009;76:813–25. [27] Yanar NM, Pettit FS, Meier GH. Failure characteristics during cyclic oxidation of yttria stabilized zirconia thermal barrier coatings deposited via electron beam physical vapor deposition on platinum aluminide an on NiCoCrAlY bond coats with processing modifications for improved performance. Metall Mater Trans A 2006;37A:1563–80. [28] Hille TS, J Nijdam T, Suiker ASJ, Turteltaub S, Sloof WG. Damage growth triggered by interface irregularities in thermal barrier coatings. Acta Mater 2009;57:2624–30. [29] Siry CW, Wanzek H, Dau C-P. Aspects of TBC service experience in aero engines. Materialwiss Werkstofftech 2001;32:650–3. [30] Roy SK, Colbe RL. Solubility of hydrogen in porous polycrystalline aluminum oxide. J Am Ceram Soc 1967;50:435–6. [31] Sinha NK. Constant-load tertiary creep in nickel-base single crystal superalloys. Mater Sci Engng A 2006;432:129–41. [32] An K, Halverson HG, Reifsnider KL, Case SW, McCord MH. Comparison of methodologies for determination of fracture strength of 8 mol% yttriastabilized zirconia electrolyte materials. J Fuel Cell Sci Technol 2005;2:99–103. [33] Cruse TA, Johnsen BP, Nagy A. Mechanical properties testing and results for thermal barrier coatings. J Therm Spray Technol 1997;6:57–66. [34] Kondoh J. Aging strengthening of 8 mol% yttria-fully-stabilized zirconia. J Alloys Compd 2004;370:285–90. [35] Choi SR, Bansal NP. Mechanical behavior of zirconia/alumina composites. Ceram Int 2005;31:39–46. [36] Johnson BP, Cruse TA, Miller RA, Brindley WJ. Compressive fatigue of a plasma sprayed ZrO2–8 wt%Y2O3 and ZrO2–10 wt%NiCoCrAlY TBC. J Engng Mater Technol 1995;117:305–10. [37] Matsuzawa M, Abe M, Horibe S. Strain rate dependence of tensile behavior and environmental effect in zirconia ceramics. ISIJ Int 2003;43:555–63. [38] Noguchi K, Fujita M, Masaki T, Mizushina M. Tensile strength of yttria-stabilized tetragonal zirconia polycrystals. J Am Ceram Soc 1989;72:1305–7. [39] Simone A, Sluys LJ. A generalized finite element method for polycrystals with discontinuous grain boundaries. Int J Numer Meth Engng 2006;67:1122–45. [40] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engng Fract Mech 2002;69:813–33. [41] Babus˘ka I, Melenk JM. The partition of unity method. Int J Numer Meth Engng 1997;40:727–58. [42] Melenk JM, Babus˘ka I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Engng 1996;139:289–314.