Audric Saillard

Mohammed Cherkaoui
email: mcherkaoui@me.gatech.edu

Laurent Capolungo

Esteban P Busso

G W Woodruff

Stress influence on high temperature oxide scale growth: modeling and investigation on a thermal barrier coating system

Keywords: modelling, oxide films, morphology, numerical simulation, thermal barrier coatings modeling, interfacial thermodynamics, high-temperature oxidation, numerical simulation, thermal barrier coatings

come

Stress influence on high temperature oxide scale growth: modeling and investigation on a thermal barrier coating system.

Introduction

Thermal barrier coatings are used in aircraft and industrial gas-turbine engine to allow a high functioning temperature while maintaining sustainable conditions for the metallic structure.

However, at high temperature, an oxide scale (commonly named thermally grown oxide, TGO) develops between the coating, usually made of Yttria-stabilized Zirconia (YSZ), and the metallic bond coat (BC), limiting the lifetime of the coating and as a consequence the long term reliability of the structure. The oxide scale, mainly constituted of Alumina oxide, presents significantly different thermal expansion and creep properties than the surrounding YSZ and BC layers. This induces high stresses and deformations upon cooling at room temperature, leading to crack nucleation and eventually to the spallation of the top coat. Thus, the growth of the TGO during high temperature operations is the most important phenomenon responsible for the TBC intrinsic failure [1,2]. However, if the global mechanisms driving the oxide scale development are well known, the factors for varying kinetics and the development of a rough and undulated scale are poorly understood.

Several studies have investigated the development of undulations at the TGO/BC interface, pointing out the key roles of oxide growth stresses and interface geometry [2][3][4][5]. These works focus on plastic and creep deformation during cooling and heating cycles, or stress-driven material diffusion within the BC [6,7], but all consider a uniform growth of the TGO. To our knowledge, no study has been performed considering a local effect of stress on the oxidation kinetics at the BC/TGO interface, which could result in the non-uniform growth of the oxide scale, a factor which has been identified as strongly detrimental to the TBC lifetime [2]. It would result in additional stress concentrations upon cooling and might lead to large local scale deformations of the oxide layer, accelerating the development of cracks and spallations. Even limited oxide scale morphology variations might become significant under the effect of thermal cycling. Therefore, it is of primary interest to understand the mechanisms leading to a nonuniform oxide growth, quantify the contributions and identify the most influential physical parameters. This should allow for the optimization of the TBC material system design in terms of composition and fabrication process, and also help predicting the location of failure.

A mechanism is proposed here relying on an influence of the stress/strain state on each side of the sharp oxide/metal interface on the local oxidation kinetics. A specific formulation is derived for this effect considering the thermodynamic force driving the sharp oxide/metal phase boundary propagation, and specially including the contributions of the large volumetric strain and the mass consumption associated with metal oxidation. A continuum mechanics/diffusion framework completes the model, accounting for the diffusion processes through the material layers and along the interface, globally driving the oxidation, as well as the large phase transformation eigenstrain development and its accommodation. The numerical simulation of this complex formulation required the development of an original numerical scheme coupling a finite element method for the calculation of the stress and diffusion field, with an external routine for the resolution of the phase boundary local composition and propagation. The non-uniform oxide scale growth and its coupling with stress development are simulated and investigated in the case of an Alumina layer growing between an electron-beam physical vapor deposited (EB-PVD) Yttria-stabilized Zirconia thermal barrier coating and a Ni-Pt-Al bond coat at 1100°C.

Formulation of the stress influence on the local oxidation kinetics

The growth of the Alumina scale in the thermal barrier coating environment occurs mainly at the oxide/metal interface through the direct oxidation of the metallic bond coat [8]. The oxygen consumption during the oxidation reaction engenders the creation of vacancies at the oxide boundary, which diffuse counter to oxygen ions towards the free surface exposed to an oxygenrich atmosphere. The oxidation chemical reaction being much faster than the species diffusion, the inner oxidation process is limited by the local availability of oxygen at the interface.

Consequently, coupled quasi-equilibrium concentrations of oxygen and vacancies rapidly establish at the interface and are maintained during the phase propagation. The mechanism proposed here for a stress influence on the local oxidation kinetics relies on a deviation from the stress-free quasi-equilibrium concentration induced by the mechanical state at the propagating interface. Since the gradient of concentration drives species diffusion, the local incoming flux of oxygen at the phase boundary is in turn affected and eventually the oxidation kinetics (oxide phase propagation) is locally modified.

The principle of a thermodynamic effect induced by the stress development accompanying metal oxidation at high temperature on the process kinetics has been reviewed by Evans [9]. He demonstrated through a simple direct influence on interface vacancy concentration that the accommodation work in response of the metal oxidation-associated volume expansion might significantly alter and even stop the uniform growth of a Zirconia film. Here this principle is extended and a complete local treatment is proposed, allowing for the simulation of a complex non-uniform phenomenon. A continuum mechanics description of a dissipative sharp interface propagation is derived to include the influence of the stress state and accommodation work on the phase boundary propagation kinetics. This thermodynamic approach, extensively developed over the last decades (see [10] for general treatment and references), allows for a direct description of the stress and composition influences on the local interface motion kinetics. Very few works have considered a mechanical coupling with phase boundary propagation during high-temperature effect of epitaxial stress on the growth kinetics of a ZrO 2 layer on a zirconium substrate at 600°C [11]. Creton et al. used a similar approach to investigate stress development associated with the oxidation of an uranium dioxide substrate at 300°C [12]. On top of the classical free energy change under small strains, the derivation presented here considers the mass diffusion fluxes, the mass consumption and the deformation work during the phase transformation process, and it incorporates the contribution of the large volumetric strain associated with the metal oxidation.

A control volume is considered, undergoing a phase change such that at time t , there is a surface, ( ) t Γ , which defines the boundary between the untransformed ( ) + Ω and transformed ( ) -Ω volumes of the material (see Figure 1). During this process, ( ) t Γ propagates through the material with a velocity Γ v , with a normal component defined by:

n v Γ = ⋅ v n (1)
where n is the normal vector to the interface.

Mechanical description

The volume is modeled with two continuum media separated by a sharp interface, as shown in Figure 1. The location vectors for each particle in the reference configuration are given by x and, at time t , by y , respectively. Then, the displacement of the particle is ( , ) ( , ) t t = u x y x x and the corresponding deformation gradient and material velocity at time t are = ∂ x F y and t = ∂ v u . The quantities in the undeformed or reference configuration will be identified by a subscript zero, see The interface is assumed coherent at the continuum mechanics description scale, implying that ( , ) t y x is continuous across the interface, but that the material velocities and the deformation gradients are discontinuous. Using subscripts + and -for the limiting values of a quantity α when the interface is approached from the + Ω and - Ω domains, respectively, its jump is [ ] [ ] on ( )

n v t = - Γ v F n . ( 2
)
For a local and updated description of the volume change associated with the phase oxidation, an effective measure is employed, I Π , obtained through: det( ) det( )

eq I V V - - + + Π = = F F (3)
where eq V -stands for the volume of the transformed and accommodated phase, equivalent to the volume of the original phase V + before its oxidation. Finally, the equilibrium condition at the interface requires that, if any surface stress phenomena are neglected, the traction, t , must be continuous across the boundary:

[ ] [ ] on ( ) t σ = = Γ t n 0 . (4)

Chemical description

The considered volume, Ω , is composed of a different phase on each side of the sharp boundary ( ) t Γ . In the considered materials at high temperature, both ionic species and defects diffuse within each phase bulk and through the interface. The two phases are assumed to present distinct chemical compositions, such that the propagation of the child phase over the parent one locally induces mass "trapping" within the formed material and "release" from the dissolved one.

The global mass conservation for a species i over a volume embedding the phase boundary is given by:

( ) [ ] { } ( ) 0 0 Ω det V i i i i I i n dC d C C v dS dt + - Γ   ∇ ⋅ + = -⋅ + -Π     ∫ ∫ J F J n (5)
where i J and i C represents the species flux vector and concentration, respectively. The terms involving the concentration in this expression denote the two contributions in mass storage during the phase boundary propagation: phase transformation and phases composition change.

The local mass balance associated with the phase boundary propagation is given by:

[ ] { } i i I i n C C v + - ⋅ = -Π J n (6)
The model of a sharp interface requires conditions relating the chemical potentials on each side of the discontinuity. At this point, there is a need to differentiate interstitial species from those occupying regular lattice sites, also known as substitutional species [10,13]. Due to the lattice constraint, i.e. the conservation of the finite number of lattice sites available for a given substitutional species, these species can only diffuse through a vacancy mechanism. Thus the flux of a substitutional species i is coupled with an opposite flux of vacancies v i : v for any substitutional species

i i i = - J J . ( 7
)
Although a specific notation for the vacancy is then introduced for each substitutional species i, it is important to note that a kind of vacancy corresponds to a given sublattice and not to an atomic species. Therefore in this notation several indices might refer to the same type of defect, which would have to be treated altogether as one species.

The diffusion potential of substitutional species, also called relative or virtual potential, is defined by [14,15]:

F
v i i i µ µ µ = - % . ( 8
)
This potential can be used instead of the true species one in energy balances to avoid dealing directly with defects. The diffusion potential of interstitial species is assumed to reduce to the chemical potential. The local chemical equilibrium at the phase boundary requires that the diffusion potentials must be continuous across the interface, for both interstitial and substitutional species:

0 i µ   =   , (9)
in which the relative potential must be used for substitutional species [10].

Driving force for the phase boundary propagation

The propagation of the phase boundary is accompanied by an energy dissipation, whose value per unit length of propagation constitutes the driving force for the process. It is derived in this section. The control volume described in Figure 1

e G P = -Ψ & . ( 10
)
Note that this formulation assumes isothermal conditions and constant kinetic energy for the body.

Power of external forces

In the absence of body forces, the power is that developed by the external forces on the

δ δ µ = ⋅ - ⋅ ∫ ∫ t v J n % , ( 11
)
where t is the traction, v is the material velocity, and i µ % and i J are the diffusion potential and flux, respectively, of any diffusing species through the external boundary of outward normal vector Ω n . Furthermore, the Einstein notation is used so that the double indices indicate a summation over the species. Then,

( ) ( ) Ω Ω σ σ ⋅ = ⋅ = ⋅ t v n v v n .
Application of the Gauss theorem on the domain Ω as described in Figure 2, embedding the surface discontinuity Γ and subjected to the action of ( )

i i σ µ = - f v J % yields: ( ) (
)

Ω e i i i i P dV dS σ µ σ µ Γ = ∇ ⋅ - + - ⋅     ∫ ∫ v J v J n % % . ( 12
) Furthermore, [ ] [ ] [ ] σ σ σ ⋅ = ⋅ = ⋅ v n n v n v , since σ σ σ + - = = n
n n over Γ by traction continuity (Equation ( 4)). This term can be rewritten in terms of the Hadamard compatibility conditions (Equation ( 2)) as follows,

[ ] [ ] n v σ σ ⋅ = - ⋅ n v n F n . ( 13
)
An admissible discontinuity in the deformation gradients across the interface Γ requires that [ ] ∝ F n . Thus, let r be an arbitrary proportionality vector so that [ ] = ⊗ F r n , where the tensor product is used:

( ) i j ij a b ⊗ = a b
. Substituting into Equation (13) leaves:

[ ] ( ) [ ] n n v v σ σ σ ⋅ = - ⋅ ⊗ = - n v n r n n : F . ( 14
)
where : represents the double contraction operator (

ij ij A B ≡ A : B
). Moreover, the continuity of tractions implies that

{ } i i i i I i n C C v µ µ + - ⋅ = -Π J n % . ( 15
)
Introducing the developments into Equation ( 12), the surface integral is:

( ) [ ] { } 1 2 e i i I i n P C C v dS σ σ µ Γ + - + - Γ   = - + - -Π     ∫ : F . ( 16
)

Time evolution of the Helmholtz free energy

The time evolution of the Helmholtz free energy, Ψ & , can be written in terms of the total free energy density W:

Ω d WdV dt Ψ = ∫ & . ( 17
)
As illustrated in Figure 3, the change of the total free energy with time includes three components: one associated with the time variation of the total free energy density W within Ω - and Ω + , i.e. Ω overall; and two others associated with the change in volume of the domains Ω -

{ } { } ( ) , , e e i i I i i I I n C C W W v dS µ µ κ γ Γ + + -- + - Γ Ψ = - -Π + -Π + ∫ & . ( 18
)
where κ represents the local interface curvature and γ I is the phase boundary energy.

Driving force for phase boundary propagation

The global rate of dissipation induced by the phase boundary propagation can now be obtained by substitution of Equations ( 16) and ( 18) into (10):

{ } { } ( ) [ ] { } , , 1
2 e i i I i i i I i i I e e I n G P C C C C W W v dS µ µ µ κ γ σ σ Γ Γ Γ + - + + -- Γ + - + -  = -Ψ = - -Π + -Π +    - + + -Π   ∫ : F & % . ( 19
)
The term within brackets represents the energy dissipated by the phase boundary local propagation for given compositions as well as stress and strain situation on each side of the sharp interface. The three first components account for the diffusion of species at the interface, the chemical reaction taking place at the phase boundary and the morphology of the discontinuity, respectively, while the two last include a contribution of the local mechanical stress state. Overall they constitute the driving force for the phase boundary propagation, f I . It can be deduced from this formulation that interstitial species do not contribute to the driving force ( i i µ µ = % so the diffusion term cancels with the chemical free energy change), but that only substitional ones do through the vacancy potential (

i i i V µ µ µ -= - %
). Assuming that oxygen is the only species diffusing through a vacancy mechanism at the interface, we obtain:

{ } ( ) [ ] { } + , , V O V O 1 2 e e I I I I f C C W W µ µ κ γ σ σ + - - + - + - = -Π + - + + -Π : F (20)
The driving force for the phase boundary local propagation then includes two stress contributions: the work required to complete the deformation associated with the phase

Phase boundary propagation kinetics

Under purely elastic and isothermal conditions, the second law of thermodynamics requires that on Γ , 0

I n f v ≥ , and implies that, when the driving force is different from zero, the motion of Γ induces dissipation. We can note that this statement is probably too restrictive for most solids undergoing inelastic deformation. In such cases, the motion of the interface Γ is expected

to start when I R f f ≥ , where R
f is a threshold value related to the nature of the interaction between the moving boundary and different types of defects. The actual driving force for the phase boundary propagation must then be defined as the difference between the two terms. Since we are interested in the stress influence only, the driving force corresponding to a flat stress-free phase boundary is denoted 0 I f in order to highlight deviation terms. Assuming furthermore that induced local oxygen concentration deviations at the interface are negligible in front of the reference value, we obtain:

( ) ( ) { } ( ) [ ] { } 0 0 0 V V O V V O , , 1
2

I I I e e I I f f C C W W µ µ µ µ κ γ σ σ + + + - - - + - + - = + - - - Π + - + + -Π : F (21)
Under the non-equilibrium thermodynamic conditions which prevail here, the kinetics of phase transformation needs to be defined as a kinetic relation between the rate n v at which the weak discontinuity moves and the driving force.

v M f = (22)
where M I is the phase boundary mobility and constitutes a temperature-dependent parameter for the phase boundary propagation process. The considered phase transformation, metal oxidation, induces a composition change. Therefore mass conservation, Equation (5), must be resolved along with the phase boundary propagation. Thus, a stress-induced deviation in the propagation kinetics is locally coupled with a composition variation, counter-balanced by mass diffusion.

Mass diffusion framework

The Alumina scale inner growth at high-temperature is controlled by the diffusion through the oxide scale of oxygen ions, assumed to operate through a vacancy mechanism. Therefore the oxide scale development, at the macroscopic level, is described by a diffusion formulation associated to mass conservation at the propagating interface. A classical formulation for composition driven species diffusion through thick film under isothermal conditions considers the chemical potential µ as the driving force. The associated defect current density is given by [17]:

i i i i C D R ϕ µ θ = - ∇ J ( 23
)
where i D ϕ is the diffusivity of species i in phase φ, R is the molar gas constant and θ is the absolute temperature in Kelvin. The chemical potential for the defect species is expressed by: ( )

ο ln i i i R c µ µ θ = + (24)
where ο i µ is the standard chemical potential of species i in a reference state of phase and temperature, and i c is the species mole fraction. Alumina scales developed at high temperature on metallic alloys exhibit a polycrystalline nature, such that two diffusion paths exist through the oxide scale: the grain bulks and the grain boundaries. Furthermore, the metal/oxide interface also constitutes a specific diffusion path.

While grain boundaries as well as the phase boundary occupy very low volume fractions, they present much higher diffusivities such that their contribution can generally not be neglected [18].

The grain boundary contribution is included classically by considering a homogenized formulation for the diffusion through the material layer. The validity of such a formulation is only direct for one-dimensional models. It is justified here for a 2D model by the fact that no particular roughness of the oxide/metal phase boundary has been reported in direct relation with grain boundary locations, despite significant grain size and difference in diffusivities between the two paths. This proves that in typical cases the fast diffusion along the interface operates a large redistribution of the incoming species flux along the boundary, resulting in a quite uniform oxide growth. The fast diffusion along the interface is included through a 1D path following the boundary between the metal and oxide phases domains. In order to obtain the diffusion flow, the path is attributed a section, w I , which stands for the phase boundary thickness, and a specific diffusivity, I i D , is applied:

( ) I I I i i i I i w D C J R µ θ = - ∇ (25)
where I i J is the flux of species i along the interface (per unit thickness of the 2D domain) and I ∇ is the tangent gradient along the interface.

Continuum mechanics framework

The possible origins of stresses associated with an oxide scale growth have been reviewed by several authors [9,19,20]. Here we focus on the large volume eigenstrain induced by the direct alloy oxidation, since it is most likely to constitute the main contribution to the local stress development at the interface. Furthermore, it induces a direct coupling between the local stress situation and the phase boundary propagation likely to significantly influence local developments. We consider here an effective measure of the volume ratio, denoted Π, possibly including additionally to the Pilling-Bedworth ratio other chemical contributions in volume change, as for instance from interdiffusion of species at the interface accompanying the oxidation reaction [21,22]. The free energy released by the oxidation reaction is such that a strain accommodation work able to suppress the process would require a stress level about two orders of magnitude higher than typical yield or fracture stresses [9]. Consequently, the volume eigenstrain associated with phase oxidation is spontaneously accommodated. This generally results in the development of an anisotropic strain, exhibiting a main component in the direction of "easiest" expansion. Such a phenomenon has been identified in experiments in which the strain developed essentially along the global direction of free expansion, i.e. towards the free surface [23]. In this situation the oxidation strain can simply be modeled as transversely isotropic with a main component on the surface normal vector. However in a more general case and in order to better describe local developments at the interface, no arbitrary anisotropy should be enforced.

Therefore we consider here the development of an isotropic volume eigenstrain, which is simultaneously accommodated by a plastic flow. The enforced true volumetric strain is:

( ) In order to assess the stress influence on the oxidation process, it is essential to provide a realistic mechanical framework. Here an elasto-visco-plastic behavior is used for each material response. The model retained combines uncoupled linear elasticity, rate-independent plasticity and creep. Small elastic strains and isotropic behaviors are assumed and the rate of deformation is taken as a strain measure. Hook's law relates the true (Cauchy) stress tensor, σ, to the elastic (logarithmic) strain component, el ε , through the fourth-order elastic stiffness tensor, C:

1 ln 3 vol ε Π = Π (26)
el = σ Cε .
A strain additive decomposition is used:

el pl vol ε Π = + + ε ε I ε ( 27
)
in which pl ε stands for the plastic accommodation strain component. Incremental plasticity formulations are used, and both rate-independent plasticity and creep are modeled through the Mises stress potential (equivalent deviatoric stress, temperature [24]. The corresponding equivalent strain rate, cr ε & % , is then given by:

cr N cr cr d A ε σ = & % % (28)
in which N cr and A cr are parameters dependent on material's composition and microstructure. and processes related to the interface propagation significantly influence the fields' evolutions and when the local propagation kinetics in turn depends on these affected values at the discontinuity. In this study aimed at investigating long-term high temperature oxidation, the interface displacement is gradual so that no fast morphological changes take place. However, the phase boundary propagation locally induces a large volume eigenstrain which considerably affects the strain and stress fields at the interface, which in turn influences the local propagation kinetics. Furthermore, the oxidation kinetics is coupled with mass diffusion, including fast transport along the interface. The resolution challenge then comes from the multi-physics treatment of a discontinuity propagation coupled with its environment. General problems of continuum mechanics are traditionally solved using the finite element method (FEM) in which a domain is modeled by a mesh that is locally refined to capture high spatial variations while keeping a smooth evolution over the elements. The definition and propagation of a sharp discontinuity in this environment are no longer straight-forward. This explains why most numerical predictions of complex phase boundary propagation / stress development couplings have been historically performed using two limited approaches. The first of these addresses simple interfacial geometries (linear, cylindrical or spherical) which are conserved in time, and allows for the use of a finite difference method over a grid evolving with the discontinuity motion [25,26]. The second approach addresses more complex geometries by employing FEM but limiting the analysis to a predefined interface propagation; this allows to design a mesh which remains adapted at all calculation time steps [8,27]. Specific numerical methods have been developed and are the subjects of extensive research for the treatment of phase boundary propagation within FEM, such as phase field [28] or level set [29]. However limitations remain on the allowed propagation kinetics formulations (diffuse interface or not field-dependent evolution) and the inclusion of local processes (such as interface diffusion). Thus the numerical treatment of the propagation of an evolving sharp interface dependent on field values in a continuum framework still constitutes a major challenge.

Numerical simulation

In order to implement the presented framework for oxide-metal phase boundary propagation at high temperature, a specific numerical tool has been developed allowing for the propagation of an evolving interface following a complex law dependent on field values. It is based on a classical FEM for the fields' resolution, performed sequentially with a specially developed external routine run in pre/post-calculations at each time step to calculate the phase boundary composition and propagation. It essentially couples the concept of phase field [30] (diffuse interface of given width described by a field variable) for the smooth treatment of a moving sharp interface over a fixed mesh, and the front tracking method [31] for the phase boundary propagation (the phase field evolution is assigned by the external routine from the mechanical and diffusion fieldsdependent propagation kinetics). A chart flow presenting the resolution scheme is provided in Figure 4. Two steps are sequentially performed for each time increment.

First, a prediction of interface displacement over the time increment is computed using diffusion fluxes along with stress and strain tensors at the sharp phase boundary extracted from the previous FEM analysis. Vacancy concentration and displacement at tracking points are resolved through time integration of the system of equations formed by: mass conservation, Equation ( 5), phase boundary propagation, Equations ( 22) and ( 21 integration, using finite differences over the interface tracking points, of the concentrationdependent mass conservation statement:

[ ] ( ) , , 2 , , 1
I i t t i t i I i t t i i t t i I w C C tD C A C B t +∆ +∆ +∆ - = ∆ ∇ - + + ⋅ ∆ J n (29) ( ) ( ) { } ( ) ( ) { } 0 1 2 1 / / 0 1 / 1 / , , 1 2 1 /
ln

avg ox avg ox ox avg ox i ox avg ox I I C C C C i ox ox i i I I I I I C C C C i i i i e i e i I A C M R C C B C M f C R f f W W σ σ θ θ κ γ σ σ - - - - + - + - = Π - = Π -Π + + +   = - + + -Π   : F
Here the upperscript i denotes a given tracking point along the interface. Equation ( 24) has been included in this development and linearized providing small deviations from the average concentration, C avg , along the phase boundary. Furthermore, C ox stands for the oxide stoichiometric value of the concentration, and C 0 for the initial value corresponding to 0 I f . Once the concentrations at the phase boundary tracking points are obtained, the displacements are calculated through the explicit time integration of the normal propagation velocity:

( ) , , 1
ox I i t t i t i i i C t A C B +∆ Π = + ∆ + x x n (30)
Then the curvatures at each tracking point are updated. The algorithm developed for time integration combines the previously described loop with a predictor-corrector scheme using adaptive time increment setting. This scheme is required to solve for the strong coupling of propagation and then concentration with interface geometry, resulting from the curvature term in the driving force and fast diffusion along the interface.

Then, data for the next FEM analysis are built. It includes the updated concentration along the interface and the evolution of the phase field, which defines the local material properties. An updated value of the field, φ i , is defined for each element of the FEM mesh from the new position of the phase boundary, as shown in Figure 5. For elements lying within the diffuse layer, In a second step, a transient coupled stress/diffusion FEM analysis is then performed, considering the displaced phase boundary as being fixed, to resolve for the time evolution of the diffusion and mechanical fields. The oxygen concentration is obtained using Fick's first law for the diffusion fluxes and a classical mass conservation statement. The material displacement, and the associated stress and strain fields, are solved through continuum mechanics, as described in section 4. The local material diffusivity only depends on the element location with respect to the node-based interface (metal or oxide side), while mechanical properties vary linearly, as defined by the phase field, from the untransformed to the transformed materials values. Therefore the interface is only diffuse for the mechanical fields' resolution.

An iterative predictor/corrector method is employed for the global resolution scheme as shown in Figure 4, with a convergence criterion on the maximum local difference in driving force for the phase boundary propagation, in order to ensure simulation convergence and accuracy despite sequentially evolving but highly coupled mechanical/diffusion and phase fields. While the developed numerical framework requires a very fine mesh over the propagation zone and is

Simulation description

The framework developed is applied to a practical case concerning a thermal barrier coating (TBC) formed by electron-beam physical vapor deposition (EB-PVD) of Yttria stabilized zirconia (YSZ) over a Pt-aluminized NiAl bond coat (BC) adapted from the model by Busso et al. [8,32]. The simulation describes the long-term isothermal growth of the oxide layer (named thermally grown oxide, TGO) at the operating temperature, 1100°C, by oxidation of the metallic bond coat at the metal/oxide interface. The simulation domain describes a cross-section perpendicular to the YSZ/BC interface plane and is meshed with generalized plane strain elements. The YSZ and BC layers are modeled, and the initial YSZ/BC interface is idealized as a wavy-type segment, as introduced by Busso et al. from scanning electron microscopy micrograph observations. Here the initial interface is described by a generic quarter sin period, as presented in Figure 6, of dimensions representative of a typical local roughness after coating deposition (8µmwidth and 3.2µm-depth) [8,33]. Symmetry boundary conditions are enforced for displacement at the lateral sides and for the thickness of the domain, which forces the system to uniformly expand in the in-plane directions (x 1 and x 3 ). Additionally, in order to account for the in-plane [32]. Finally the simulation domain is fixed at the lower boundary but free to expand at the top surface. The materials parameters are given in Table 1 for the bond coat [34], the thermally grown oxide [35] and the YSZ coating [8,35]. For the oxidation volume eigenstrain, The Pilling-Bedworth ratio is taken from the primary oxidation reaction as 1.28 Π =

= = & &
, which gives an isotropic volume strain of 0.08

vol ε Π ≈
. Since from the numerical scheme implemented the propagating metal/oxide phase boundary description within the finite elements domain relies on the mesh, the oxidation zone is finely meshed using 36000 1st order elements of dimensions 0.02 0.04 µm.

The long term inner TGO growth at high temperature is limited by the oxygen inward diffusion through the developed scale. The initial configuration is here obtained from the uniform growth of a thin oxide layer (0.5 µm) from a stress-free system in a time obtained from the empirical kinetics law in [8] (about 30 hr). A fixed unit concentration is considered at the free surface for oxygen, providing that the absorption process is infinitely fast compared to the ionic transport. The oxygen solubility in the bond coat is assumed negligible and constant diffusivities are employed. Here effective diffusion parameters are formed by multiplying the diffusivities with the stress-free relative vacancy concentration at the phase boundary,

V V O / C C δ =
, which constitutes the scaling factor for the diffusion gradient and is enforced at the initial step. This parameter is fitted for the TGO so that the scale growth kinetics matches globally the value of the experimental data over the simulated time. A value of ≈ is thus assigned. The diffusion parameter for the YSZ layer is obtained through the ratio of its diffusivity to the one of the Alumina scale, which gives a ratio around 7 10 from [36] and [37]. Such a high value makes

n I I v t M f =
, when the stress and curvature contributions are neglected.

Since the inner oxidation is globally driven by the diffusion of oxygen ions through the oxide scale, the mass balance must be achieved at the phase boundary, so that the mass rate per unit surface consumed in the oxidation reaction matches the incoming diffusion flux. This gives from Equation ( 6):

[ ] 0 O 0 O O at I I I M f t t C C + - ⋅ = = -Π J n (31)
Thus this product is calculated from the initial oxide scale thickness and chemical concentrations, along with the oxide diffusivity data, and only one of the two parameters is to estimate.

Observation of the oxide/metal interface morphology has been used to set 0 I f . The curvaturerelated thermodynamic force, I κ γ , can easily be quantitatively estimated since the interfacial energy is known and the curvature range can be identified from cross-section micrographs. Since

Results and discussion

The results from the numerical simulations are presented in this section. They include the stress and curvature effects on the local propagation kinetics of the oxide-metal phase boundary, the volume eigenstrain generation from the direct metal oxidation at the inner interface and the time-dependent and time-independent behaviors of the different materials, and the specific mass diffusion along the propagating phase boundary.

The oxidized bond coat layer in the TBC environment as predicted after 100 and 335 hours of isothermal oxidation is presented on In order to investigate the stress situation which lead to this non-uniform oxidation, the two contributions of the stress influence on the process kinetics, namely the work required to complete the deformation associated with the metal phase oxidation at the stressed interface, although it can be seen that it increases from the peak to the valley region. The higher value at the valley compared to the peak is due to the fact that the TGO and the YSZ coating are less strain compliant than the metallic bond coat, inducing a more limited in-plane accommodation of the oxidation volume eigenstrain at the valley in comparison with the peak. On top of being rather uniform, the data are also quite constant in time, which backs up the validity of the model of a phase transformation eigenstrain exhibiting a constant anisotropy with the main direction along the local normal to the phase boundary [8]. However, the variations of the main transformation strain component along the interface, even if limited, greatly influence the local traction. From this analysis, the traction is clearly the main component of the local oxidation kinetics variations.

( ) [ ]
It is essentially compressive in the valley region, inducing the lower oxidation kinetics, slightly tensile near the peak resulting in the faster than average oxidation and vanishes at the mid-slope.

It is interesting to notice that the stress contribution gradient is significant over the 0-100 hours This highlights the progressive development of the oxide scale, but stress effects at long term are also reinforced by the role of diffusion along the interface. Indeed, it becomes more and more influent with time since the oxide thickness increase leads to a decrease of the arriving oxygen flux, thus allowing for a larger mass redistribution along the phase boundary.

The non-uniform morphology of the oxide scale has been identified as a strongly detrimental factor for the TBC lifetime [2]. Figure 9 shows the evolution with oxidation time of the ratio of the maximum local scale thickness to the minimum value. The respective influences of the oxidemetal interfacial energy and diffusivity are also presented. It can be seen that very significant thickness non-uniformity develops with time in the reference case, with a maximum local thickness getting close to twice the minimum one after 500 hours. The development is quite linear with time, with a rate increase at around 200 hours. As noticed previously, the nonuniformity is limited after 100 hours of oxidation, with a ratio value around 1.1 (10% thickness difference), but more significant after 335 hours where it reaches almost 1.5 (50%). These values are consistent with local thickness ratios than can be measured on SEM micrographs, see for instance [33]. The fast mass diffusion along the inner interface is a key process contributing to the significant variations in local oxidation kinetics. Therefore the phase boundary diffusivity influences the rate of non-uniform morphological development through the kinetics of mass redistribution. However, it can be seen on Figure 9 that the effect of an interface diffusivity variation by an order of magnitude from the reference case considered is limited on the nonuniformity development. This indicates that the diffusion process is not limiting the interface evolution in this case. Furthermore, no clear influence can be deduced since a decrease of the diffusivity eventually leads to an increase in non-uniformity. This situation is due to the strong development, since the equilibrium radius is more than twenty times lower than the one of the initial interface. The phenomenon is only limited by the process kinetics. The effect of a curvature increase is also presented, showing again the very beneficial contribution. This graph provides two guidelines for the TBC material system design in order to limit non-uniformities development: to increase the interfacial energy, as stated earlier, and to increase the stressinduced flow rate at the phase boundary in order to limit the local traction magnitude. These properties might be tailored through the use of additives, notably reactive elements, which are known to segregate at the oxide/metal interface and improve its adhesion [41]. Globally, the maximum traction magnitude at high temperature is essentially set by the weak bond coat strength, which should then be as low as possible in this environment. Lowering the initial interface roughness would also be beneficial since it would result in a decrease of the initial maximum traction [33], thus slowing down the non-uniformity development kinetics.

Conclusions

A study of oxide scale growth developing a non-uniform morphology in an EB-PVD thermal Figure 10. Curvature radius vs traction amplitude.

Formulation of the stress influence on the local oxidation kinetics

The growth of the Alumina scale in the thermal barrier coating environment occurs mainly at the oxide/metal interface through the direct oxidation of the metallic bond coat [8]. The oxygen effect of epitaxial stress on the growth kinetics of a ZrO 2 layer on a zirconium substrate at 600°C [11]. Creton et al. used a similar approach to investigate stress development associated with the oxidation of an uranium dioxide substrate at 300°C [12]. On top of the classical free energy change under small strains, the derivation presented here considers the mass diffusion fluxes, the mass consumption and the deformation work during the phase transformation process, and it incorporates the contribution of the large volumetric strain associated with the metal oxidation.

A control volume is considered, undergoing a phase change such that at time t , there is a surface, ( ) t Γ , which defines the boundary between the untransformed ( ) + Ω and transformed ( ) -Ω volumes of the material (see Figure 1). During this process, ( ) t Γ propagates through the material with a velocity Γ v , with a normal component defined by:

n v Γ = ⋅ v n (1)
where n is the normal vector to the interface.

Mechanical description

The volume is modeled with two continuum media separated by a sharp interface, as shown in Figure 1. The location vectors for each particle in the reference configuration are given by x and, at time t , by y , respectively. Then, the displacement of the particle is ( , )

6/34 Authors

The interface is assumed coherent at the continuum mechanics description scale, implying that ( , ) t y x is continuous across the interface, but that the material velocities and the deformation gradients are discontinuous. Using subscripts + and -for the limiting values of a quantity α when the interface is approached from the + Ω and - Ω domains, respectively, its jump is

n v t = - Γ v F n . ( 2
)
For a local and updated description of the volume change associated with the phase oxidation, an effective measure is employed, I Π , obtained through: det( ) det( )

eq I V V - - + + Π = = F F (3)
where eq Vstands for the volume of the transformed and accommodated phase, equivalent to the volume of the original phase V + before its oxidation. Finally, the equilibrium condition at the interface requires that, if any surface stress phenomena are neglected, the traction, t , must be continuous across the boundary:

[ ] [ ] on ( ) t σ = = Γ t n 0 . ( 4
)

Chemical description

The considered volume, Ω , is composed of a different phase on each side of the sharp boundary ( ) t Γ . In the considered materials at high temperature, both ionic species and defects diffuse within each phase bulk and through the interface. The two phases are assumed to present distinct chemical compositions, such that the propagation of the child phase over the parent one The global mass conservation for a species i over a volume embedding the phase boundary is given by:

( ) [ ] { } ( ) 0 0 Ω det V i i i i I i n dC d C C v dS dt + - Γ   ∇ ⋅ + = -⋅ + -Π     ∫ ∫ J F J n (5)
where i J and i C represents the species flux vector and concentration, respectively. The terms involving the concentration in this expression denote the two contributions in mass storage during the phase boundary propagation: phase transformation and phases composition change.

The local mass balance associated with the phase boundary propagation is given by:

[ ] { } i i I i n C C v + - ⋅ = -Π J n (6)
The model of a sharp interface requires conditions relating the chemical potentials on each side of the discontinuity. At this point, there is a need to differentiate interstitial species from those occupying regular lattice sites, also known as substitutional species [10,13]. Due to the lattice constraint, i.e. the conservation of the finite number of lattice sites available for a given substitutional species, these species can only diffuse through a vacancy mechanism. Thus the flux of a substitutional species i is coupled with an opposite flux of vacancies v i : v for any substitutional species

i i i = - J J . ( 7
)
Although a specific notation for the vacancy is then introduced for each substitutional species i, it is important to note that a kind of vacancy corresponds to a given sublattice and not to an atomic species. Therefore in this notation several indices might refer to the same type of defect, which would have to be treated altogether as one species.

The diffusion potential of substitutional species, also called relative or virtual potential, is defined by [14,15]:

v i i i µ µ µ = - ɶ . ( 8
)
This potential can be used instead of the true species one in energy balances to avoid dealing directly with defects. The diffusion potential of interstitial species is assumed to reduce to the chemical potential. The local chemical equilibrium at the phase boundary requires that the diffusion potentials must be continuous across the interface, for both interstitial and substitutional species:

0 i µ   =   , ( 9
)
in which the relative potential must be used for substitutional species [10].

Driving force for the phase boundary propagation

The propagation of the phase boundary is accompanied by an energy dissipation, whose value per unit length of propagation constitutes the driving force for the process. It is derived in this section. The control volume described in Figure 1 is considered to undergo external mechanical forces and mass fluxes. The global rate of dissipation, G, is given by the difference between the rate of external work (in a general thermodynamic sense), e P , and the rate of change of the total Helmholtz free energy, Ψ ɺ :

e G P = -Ψ ɺ . ( 10
)
Note that this formulation assumes isothermal conditions and constant kinetic energy for the body.

Power of external forces

In the absence of body forces, the power is that developed by the external forces on the boundary Ω δ and by the flow of diffusing species across Ω δ :

F
δ δ µ = ⋅ - ⋅ ∫ ∫ t v J n ɶ , ( 11
)
where t is the traction, v is the material velocity, and i µ ɶ and i J are the diffusion potential and flux, respectively, of any diffusing species through the external boundary of outward normal vector Ω n . Furthermore, the Einstein notation is used so that the double indices indicate a summation over the species. Then,

( ) ( ) Ω Ω σ σ ⋅ = ⋅ = ⋅ t v n v v n .
Application of the Gauss theorem on the domain Ω as described in Figure 2, embedding the surface discontinuity Γ and subjected to the action of ( )

i i σ µ = - f v J ɶ yields: ( ) (
)

Ω e i i i i P dV dS σ µ σ µ Γ = ∇ ⋅ - + - ⋅     ∫ ∫ v J v J n ɶ ɶ . ( 12
) Furthermore, [ ] [ ] [ ] σ σ σ ⋅ = ⋅ = ⋅ v n n v n v , since σ σ σ + - = = n
n n over Γ by traction continuity (Equation ( 4)). This term can be rewritten in terms of the Hadamard compatibility conditions (Equation ( 2)) as follows,

[ ] [ ] n v σ σ ⋅ = - ⋅ n v n F n . ( 13
)
An admissible discontinuity in the deformation gradients across the interface Γ requires that [ ] ∝ F n . Thus, let r be an arbitrary proportionality vector so that [ ] = ⊗ F r n , where the tensor product is used:

( ) i j ij a b ⊗ = a b
. Substituting into Equation ( 13) leaves:

[ ] ( ) [ ] n n v v σ σ σ ⋅ = - ⋅ ⊗ = - n v n r n n : F . ( 14
)
where : represents the double contraction operator (

ij ij A B ≡ A : B
). Moreover, the continuity of tractions implies that

{ } i i i i I i n C C v µ µ + - ⋅ = -Π J n ɶ . ( 15
)
Introducing the developments into Equation ( 12), the surface integral is:

( ) [ ] { } 1 2 e i i I i n P C C v dS σ σ µ Γ + - + - Γ   = - + - -Π     ∫ : F . ( 16
)

Time evolution of the Helmholtz free energy

The time evolution of the Helmholtz free energy, Ψ ɺ , can be written in terms of the total free energy density W:

Ω d WdV dt Ψ = ∫ ɺ . ( 17
)
As illustrated in Figure 3, the change of the total free energy with time includes three components: one associated with the time variation of the total free energy density W within Ω - and Ω + , i.e. Ω overall; and two others associated with the change in volume of the domains Ω -

{ } { } ( ) , , e e i i I i i I I n C C W W v dS µ µ κ γ Γ + + -- + - Γ Ψ = - -Π + -Π + ∫ ɺ . ( 18
)
where κ represents the local interface curvature and γ I is the phase boundary energy.

Driving force for phase boundary propagation

The global rate of dissipation induced by the phase boundary propagation can now be obtained by substitution of Equations ( 16) and ( 18) into ( 10):

{ } { } ( ) [ ] { } , , 1
2 e i i I i i i I i i I e e I n G P C C C C W W v dS µ µ µ κ γ σ σ Γ Γ Γ + - + + -- Γ + - + -  = -Ψ = - -Π + -Π +    - + + -Π   ∫ : F ɺ ɶ . ( 19
)
The term within brackets represents the energy dissipated by the phase boundary local propagation for given compositions as well as stress and strain situation on each side of the sharp interface. The three first components account for the diffusion of species at the interface, the chemical reaction taking place at the phase boundary and the morphology of the discontinuity, respectively, while the two last include a contribution of the local mechanical stress state. Overall they constitute the driving force for the phase boundary propagation, f I . It can be deduced from this formulation that interstitial species do not contribute to the driving force ( i i µ µ = ɶ so the diffusion term cancels with the chemical free energy change), but that only substitional ones do through the vacancy potential (

i i i V µ µ µ -= - ɶ
). Assuming that oxygen is the only species diffusing through a vacancy mechanism at the interface, we obtain:

{ } ( ) [ ] { } + , , V O V O 1 2 e e I I I I f C C W W µ µ κ γ σ σ + - - + - + - = -Π + - + + -Π : F (20)
The driving force for the phase boundary local propagation then includes two stress contributions: the work required to complete the deformation associated with the phase

Phase boundary propagation kinetics

Under purely elastic and isothermal conditions, the second law of thermodynamics requires that on Γ , 0

I n f v ≥ , and implies that, when the driving force is different from zero, the motion of Γ induces dissipation. We can note that this statement is probably too restrictive for most solids undergoing inelastic deformation. In such cases, the motion of the interface Γ is expected

to start when I R f f ≥ , where R
f is a threshold value related to the nature of the interaction between the moving boundary and different types of defects. The actual driving force for the phase boundary propagation must then be defined as the difference between the two terms. Since we are interested in the stress influence only, the driving force corresponding to a flat stress-free phase boundary is denoted 0 I f in order to highlight deviation terms. Assuming furthermore that induced local oxygen concentration deviations at the interface are negligible in front of the reference value, we obtain:

( ) ( ) { } ( ) [ ] { } 0 0 0 V V O V V O , , 1 2
I I I e e I I f f C C W W µ µ µ µ κ γ σ σ + + + - - - + - + - = + - - - Π + - + + -Π : F (21)
Under the non-equilibrium thermodynamic conditions which prevail here, the kinetics of phase transformation needs to be defined as a kinetic relation between the rate n v at which the weak discontinuity moves and the driving force. Considering a quasi-static propagation process,

0 I f ≈ ,
v M f = (22)
where M I is the phase boundary mobility and constitutes a temperature-dependent parameter for the phase boundary propagation process. The considered phase transformation, metal oxidation, induces a composition change. Therefore mass conservation, Equation ( 5), must be resolved along with the phase boundary propagation. Thus, a stress-induced deviation in the propagation kinetics is locally coupled with a composition variation, counter-balanced by mass diffusion.

Mass diffusion framework

The Alumina scale inner growth at high-temperature is controlled by the diffusion through the oxide scale of oxygen ions, assumed to operate through a vacancy mechanism. Therefore the oxide scale development, at the macroscopic level, is described by a diffusion formulation associated to mass conservation at the propagating interface. A classical formulation for composition driven species diffusion through thick film under isothermal conditions considers the chemical potential µ as the driving force. The associated defect current density is given by [17]:

i i i i C D R ϕ µ θ = - ∇ J (23)
where i D ϕ is the diffusivity of species i in phase φ, R is the molar gas constant and θ is the absolute temperature in Kelvin. The chemical potential for the defect species is expressed by: ( )

ο ln i i i R c µ µ θ = + (24)
where ο i µ is the standard chemical potential of species i in a reference state of phase and temperature, and i c is the species mole fraction. Alumina scales developed at high temperature on metallic alloys exhibit a polycrystalline nature, such that two diffusion paths exist through the oxide scale: the grain bulks and the grain boundaries. Furthermore, the metal/oxide interface also constitutes a specific diffusion path.

While grain boundaries as well as the phase boundary occupy very low volume fractions, they present much higher diffusivities such that their contribution can generally not be neglected [18].

The grain boundary contribution is included classically by considering a homogenized formulation for the diffusion through the material layer. The validity of such a formulation is only direct for one-dimensional models. It is justified here for a 2D model by the fact that no particular roughness of the oxide/metal phase boundary has been reported in direct relation with grain boundary locations, despite significant grain size and difference in diffusivities between the two paths. This proves that in typical cases the fast diffusion along the interface operates a large redistribution of the incoming species flux along the boundary, resulting in a quite uniform oxide growth. The fast diffusion along the interface is included through a 1D path following the boundary between the metal and oxide phases domains. In order to obtain the diffusion flow, the path is attributed a section, w I , which stands for the phase boundary thickness, and a specific diffusivity, I i D , is applied:

( ) I I I i i i I i w D C J R µ θ = - ∇ (25)
where I i J is the flux of species i along the interface (per unit thickness of the 2D domain) and I ∇ is the tangent gradient along the interface.

Continuum mechanics framework

The possible origins of stresses associated with an oxide scale growth have been reviewed by several authors [9,19,20]. Here we focus on the large volume eigenstrain induced by the direct alloy oxidation, since it is most likely to constitute the main contribution to the local stress development at the interface. Furthermore, it induces a direct coupling between the local stress situation and the phase boundary propagation likely to significantly influence local developments. We consider here an effective measure of the volume ratio, denoted Π, possibly including additionally to the Pilling-Bedworth ratio other chemical contributions in volume change, as for instance from interdiffusion of species at the interface accompanying the oxidation reaction [21,22]. The free energy released by the oxidation reaction is such that a strain accommodation work able to suppress the process would require a stress level about two orders of magnitude higher than typical yield or fracture stresses [9]. Consequently, the volume eigenstrain associated with phase oxidation is spontaneously accommodated. This generally results in the development of an anisotropic strain, exhibiting a main component in the direction of "easiest" expansion. Such a phenomenon has been identified in experiments in which the strain developed essentially along the global direction of free expansion, i.e. towards the free surface [23]. In this situation the oxidation strain can simply be modeled as transversely isotropic with a main component on the surface normal vector. However in a more general case and in order to better describe local developments at the interface, no arbitrary anisotropy should be enforced.

Therefore we consider here the development of an isotropic volume eigenstrain, which is simultaneously accommodated by a plastic flow. The enforced true volumetric strain is:

( ) In order to assess the stress influence on the oxidation process, it is essential to provide a realistic mechanical framework. Here an elasto-visco-plastic behavior is used for each material response. The model retained combines uncoupled linear elasticity, rate-independent plasticity and creep. Small elastic strains and isotropic behaviors are assumed and the rate of deformation is taken as a strain measure. Hook's law relates the true (Cauchy) stress tensor, σ, to the elastic (logarithmic) strain component, el ε , through the fourth-order elastic stiffness tensor, C:

1 ln 3 vol ε Π = Π (26)
el = σ Cε .
A strain additive decomposition is used:

el pl vol ε Π = + + ε ε I ε (27)
in which pl ε stands for the plastic accommodation strain component. Incremental plasticity formulations are used, and both rate-independent plasticity and creep are modeled through the Mises stress potential (equivalent deviatoric stress,

3 2 d d d σ = σ : σ ɶ
) and associated flow. The creep relaxation of the considered oxide and metal phases follows a power law at high temperature [24]. The corresponding equivalent strain rate, cr ε ɺ ɶ , is then given by:

cr N cr cr d A ε σ = ɺ ɶ ɶ (28)
in which N cr and A cr are parameters dependent on material's composition and microstructure. tasks: calculating the interface displacement following a given kinetics and solving for the continuous fields on each side of the discontinuity. Strong couplings appear when conservations and processes related to the interface propagation significantly influence the fields' evolutions and when the local propagation kinetics in turn depends on these affected values at the discontinuity. In this study aimed at investigating long-term high temperature oxidation, the interface displacement is gradual so that no fast morphological changes take place. However, the phase boundary propagation locally induces a large volume eigenstrain which considerably affects the strain and stress fields at the interface, which in turn influences the local propagation kinetics. Furthermore, the oxidation kinetics is coupled with mass diffusion, including fast transport along the interface. The resolution challenge then comes from the multi-physics treatment of a discontinuity propagation coupled with its environment. General problems of continuum mechanics are traditionally solved using the finite element method (FEM) in which a domain is modeled by a mesh that is locally refined to capture high spatial variations while keeping a smooth evolution over the elements. The definition and propagation of a sharp discontinuity in this environment are no longer straight-forward. This explains why most numerical predictions of complex phase boundary propagation / stress development couplings have been historically performed using two limited approaches. The first of these addresses simple interfacial geometries (linear, cylindrical or spherical) which are conserved in time, and allows for the use of a finite difference method over a grid evolving with the discontinuity motion [25,26]. The second approach addresses more complex geometries by employing FEM but limiting the analysis to a predefined interface propagation; this allows to design a mesh which remains adapted at all calculation time steps [8,27]. Specific numerical methods have been propagation within FEM, such as phase field [28] or level set [29]. However limitations remain on the allowed propagation kinetics formulations (diffuse interface or not field-dependent evolution) and the inclusion of local processes (such as interface diffusion). Thus the numerical treatment of the propagation of an evolving sharp interface dependent on field values in a continuum framework still constitutes a major challenge.

In order to implement the presented framework for oxide-metal phase boundary propagation at high temperature, a specific numerical tool has been developed allowing for the propagation of an evolving interface following a complex law dependent on field values. It is based on a classical FEM for the fields' resolution, performed sequentially with a specially developed external routine run in pre/post-calculations at each time step to calculate the phase boundary composition and propagation. It essentially couples the concept of phase field [30] (diffuse interface of given width described by a field variable) for the smooth treatment of a moving sharp interface over a fixed mesh, and the front tracking method [31] for the phase boundary propagation (the phase field evolution is assigned by the external routine from the mechanical and diffusion fieldsdependent propagation kinetics). A chart flow presenting the resolution scheme is provided in First, a prediction of interface displacement over the time increment is computed using diffusion fluxes along with stress and strain tensors at the sharp phase boundary extracted from the previous FEM analysis. Vacancy concentration and displacement at tracking points are resolved through time integration of the system of equations formed by: mass conservation, Equation ( 5), phase boundary propagation, Equations ( 22) and ( 21), and along-interface diffusion, Equation (25)

( ) ( ) { } ( ) ( ) { } 0 1 2 1 / / 0 1 / 1 / , , 1 2 1 /
A C M R C C B C M f C R f f W W σ σ θ θ κ γ σ σ - - - - + - + - = Π - = Π -Π + + +   = - + + -Π   : F
Here the upperscript i denotes a given tracking point along the interface. Equation ( 24) has been included in this development and linearized providing small deviations from the average concentration, C avg , along the phase boundary. Furthermore, C ox stands for the oxide stoichiometric value of the concentration, and C 0 for the initial value corresponding to 0 I f . Once the concentrations at the phase boundary tracking points are obtained, the displacements are calculated through the explicit time integration of the normal propagation velocity:

( ) , , 1
ox I i t t i t i i i C t A C B +∆ Π = + ∆ + x x n (30)
Then the curvatures at each tracking point are updated. The algorithm developed for time integration combines the previously described loop with a predictor-corrector scheme using adaptive time increment setting. This scheme is required to solve for the strong coupling of propagation and then concentration with interface geometry, resulting from the curvature term in the driving force and fast diffusion along the interface.

Then, data for the next FEM analysis are built. It includes the updated concentration along the interface and the evolution of the phase field, which defines the local material properties. An updated value of the field, φ i , is defined for each element of the FEM mesh from the new position of the phase boundary, as shown in Figure 5. For elements lying within the diffuse layer, An iterative predictor/corrector method is employed for the global resolution scheme as shown in Figure 4, with a convergence criterion on the maximum local difference in driving force for the phase boundary propagation, in order to ensure simulation convergence and accuracy despite sequentially evolving but highly coupled mechanical/diffusion and phase fields. While the developed numerical framework requires a very fine mesh over the propagation zone and is

Simulation description

The framework developed is applied to a practical case concerning a thermal barrier coating (TBC) formed by electron-beam physical vapor deposition (EB-PVD) of Yttria stabilized zirconia (YSZ) over a Pt-aluminized NiAl bond coat (BC) adapted from the model by Busso et al. [8,32]. The simulation describes the long-term isothermal growth of the oxide layer (named thermally grown oxide, TGO) at the operating temperature, 1100°C, by oxidation of the metallic bond coat at the metal/oxide interface. The simulation domain describes a cross-section perpendicular to the YSZ/BC interface plane and is meshed with generalized plane strain elements. The YSZ and BC layers are modeled, and the initial YSZ/BC interface is idealized as a wavy-type segment, as introduced by Busso et al. from scanning electron microscopy micrograph observations. Here the initial interface is described by a generic quarter sin period, as presented in Figure 6, of dimensions representative of a typical local roughness after coating deposition (8µmwidth and 3.2µm-depth) [8,33]. Symmetry boundary conditions are enforced for displacement at the lateral sides and for the thickness of the domain, which forces the system to uniformly expand in the in-plane directions (x 1 and x 3 ). Additionally, in order to account for the in-plane 1 for the bond coat [34], the thermally grown oxide [35] and the YSZ coating [8,35]. For the oxidation volume eigenstrain, The Pilling-Bedworth ratio is taken from the primary oxidation reaction as 1.28 Π =

, which gives an isotropic volume strain of 0.08 vol ε Π ≈ . Since from the numerical scheme implemented the propagating metal/oxide phase boundary description within the finite elements domain relies on the mesh, the oxidation zone is finely meshed using 36000 1st order elements of dimensions 0.02 0.04 µm.

The long term inner TGO growth at high temperature is limited by the oxygen inward diffusion through the developed scale. The initial configuration is here obtained from the uniform growth of a thin oxide layer (0.5 µm) from a stress-free system in a time obtained from the empirical kinetics law in [8] (about 30 hr). A fixed unit concentration is considered at the free surface for oxygen, providing that the absorption process is infinitely fast compared to the ionic transport. The oxygen solubility in the bond coat is assumed negligible and constant diffusivities are employed. Here effective diffusion parameters are formed by multiplying the diffusivities with the stress-free relative vacancy concentration at the phase boundary,

V V O / C C δ =
, which constitutes the scaling factor for the diffusion gradient and is enforced at the initial step. This parameter is fitted for the TGO so that the scale growth kinetics matches globally the value of the experimental data over the simulated time. A value of ≈ is thus assigned. The diffusion parameter for the YSZ layer is obtained through the ratio of its diffusivity to the one of the Alumina scale, which gives a ratio around 7 10 from [36] and [37]. Such a high value makes Since the inner oxidation is globally driven by the diffusion of oxygen ions through the oxide scale, the mass balance must be achieved at the phase boundary, so that the mass rate per unit surface consumed in the oxidation reaction matches the incoming diffusion flux. This gives from Equation ( 6):

[ ] 0 O 0 O O at I I I M f t t C C + - ⋅ = = -Π J n (31)
Thus this product is calculated from the initial oxide scale thickness and chemical concentrations, along with the oxide diffusivity data, and only one of the two parameters is to estimate.

Observation of the oxide/metal interface morphology has been used to set 0 I f . The curvaturerelated thermodynamic force, I κ γ , can easily be quantitatively estimated since the interfacial energy is known and the curvature range can be identified from cross-section micrographs. Since

Results and discussion

The results from the numerical simulations are presented in this section. They include the stress and curvature effects on the local propagation kinetics of the oxide-metal phase boundary, the volume eigenstrain generation from the direct metal oxidation at the inner interface and the time-dependent and time-independent behaviors of the different materials, and the specific mass diffusion along the propagating phase boundary.

The oxidized bond coat layer in the TBC environment as predicted after 100 and 335 hours of isothermal oxidation is presented on In order to investigate the stress situation which lead to this non-uniform oxidation, the two contributions of the stress influence on the process kinetics, namely the work required to complete the deformation associated with the metal phase oxidation at the stressed interface, although it can be seen that it increases from the peak to the valley region. The higher value at the valley compared to the peak is due to the fact that the TGO and the YSZ coating are less strain compliant than the metallic bond coat, inducing a more limited in-plane accommodation of the oxidation volume eigenstrain at the valley in comparison with the peak. On top of being rather uniform, the data are also quite constant in time, which backs up the validity of the model of a phase transformation eigenstrain exhibiting a constant anisotropy with the main direction along the local normal to the phase boundary [8]. However, the variations of the main transformation strain component along the interface, even if limited, greatly influence the local traction. From this analysis, the traction is clearly the main component of the local oxidation kinetics variations.

It is essentially compressive in the valley region, inducing the lower oxidation kinetics, slightly tensile near the peak resulting in the faster than average oxidation and vanishes at the mid-slope.

It is interesting to notice that the stress contribution gradient is significant over the 0-100 hours This highlights the progressive development of the oxide scale, but stress effects at long term are also reinforced by the role of diffusion along the interface. Indeed, it becomes more and more influent with time since the oxide thickness increase leads to a decrease of the arriving oxygen flux, thus allowing for a larger mass redistribution along the phase boundary.

The non-uniform morphology of the oxide scale has been identified as a strongly detrimental factor for the TBC lifetime [2]. Figure 9 shows the evolution with oxidation time of the ratio of the maximum local scale thickness to the minimum value. The respective influences of the oxidemetal interfacial energy and diffusivity are also presented. It can be seen that very significant thickness non-uniformity develops with time in the reference case, with a maximum local thickness getting close to twice the minimum one after 500 hours. The development is quite linear with time, with a rate increase at around 200 hours. As noticed previously, the nonuniformity is limited after 100 hours of oxidation, with a ratio value around 1.1 (10% thickness difference), but more significant after 335 hours where it reaches almost 1.5 (50%). These values are consistent with local thickness ratios than can be measured on SEM micrographs, see for instance [33]. The fast mass diffusion along the inner interface is a key process contributing to the significant variations in local oxidation kinetics. Therefore the phase boundary diffusivity influences the rate of non-uniform morphological development through the kinetics of mass redistribution. However, it can be seen on Figure 9 that the effect of an interface diffusivity variation by an order of magnitude from the reference case considered is limited on the nonuniformity development. This indicates that the diffusion process is not limiting the interface evolution in this case. Furthermore, no clear influence can be deduced since a decrease of the diffusivity eventually leads to an increase in non-uniformity. This situation is due to the strong development, since the equilibrium radius is more than twenty times lower than the one of the initial interface. The phenomenon is only limited by the process kinetics. The effect of a curvature increase is also presented, showing again the very beneficial contribution. This graph provides two guidelines for the TBC material system design in order to limit non-uniformities development: to increase the interfacial energy, as stated earlier, and to increase the stressinduced flow rate at the phase boundary in order to limit the local traction magnitude. These properties might be tailored through the use of additives, notably reactive elements, which are known to segregate at the oxide/metal interface and improve its adhesion [41]. Globally, the maximum traction magnitude at high temperature is essentially set by the weak bond coat strength, which should then be as low as possible in this environment. Lowering the initial interface roughness would also be beneficial since it would result in a decrease of the initial maximum traction [33], thus slowing down the non-uniformity development kinetics.

Conclusions

A study of oxide scale growth developing a non-uniform morphology in an EB-PVD thermal

  . Garruchet et al. derived an interfacial thermodynamic model to investigate an

Figure 1 .

 1 Figure 1. Since continuum media are considered, any volume change is the result of a mechanical

  of the velocity and the deformation gradient can be related through the Hadamard compatibility condition:

  is considered to undergo external mechanical forces and mass fluxes. The global rate of dissipation, G, is given by the difference between the rate of external work (in a general thermodynamic sense), e P , and the rate of change of the total Helmholtz free energy, Ψ & :

  continuity of diffusion potentials at the interface, Equation (9), along with the local mass balance, Equation (6), yields:

C

  and Ω + , due to the motion of the external boundary on one hand, and of the interface Γ on the other hand. Only this last term is related to the phase boundary propagation and is then of interest here. The total free energy density W includes contributions from the stored elastic energy density, e W , and the chemical free energy density, stand for the chemical potential and concentration, respectively, of any species i present in the phase considered. Note that the Einstein notation is again used so that the double indices indicate a summation over the species. At the interface, a specific thermodynamic force induced by the gradient of local curvature along the interface applies, tending to flatten the boundary in order to reduce its surface and thus the global interface energy[16]. This effect can be directly introduced in the time evolution of the free energy induced by the phase boundary propagation, which then overall can be expressed by:

  and the change in elastic energy of the transformed volume over the process,

  flow. The creep relaxation of the considered oxide and metal phases follows a power law at high

  ), and along-interface diffusion, Equation(25). Assuming constant stresses and strains, normal fluxes and phase boundary curvature, the concentration transient evolution is resolved by direct implicit time

  the normal distance to the phase boundary line allows building a field evolving linearly over the finite thickness, as defined in the figure. The updated field is not directly enforced, but will result from a progressive time-dependent variation in order to ensure a smooth phase transformation, as defined in Figure5. Times for the onset of the local phase field evolution as well as the corresponding transformation kinetics are constructed based on phase boundary propagation increments. Additionnaly, a node-based interface is defined through identification of a path formed by FEM mesh nodes following the true phase boundary. It is used to enforce the oxygen concentration boundary conditions for the next diffusion analysis, interpolated from the values computed at the interface tracking points.

  interface displacements (due to the sequential evolutions of the interface and the domain fields), it allows for the treatment of complex phase boundary composition and propagation evolutions dependent on discontinuous field values and coupled with specific interface processes such as fast mass diffusion and volume eigenstrain generation. It thus meets the specific and stringent requirements of the oxidation formulation. The FEM computations have been performed using the commercial code ABAQUS, while specific programs developed under MATLAB have been used for the pre/post-calculations.

  , as it is known, transparent for oxygen diffusion. Symmetry boundary conditions are applied on the lateral sides of the domain for diffusion (zero flux normal to the boundary). To the author's knowledge, no values are available in the literature for the vacancy diffusivity along the BC/TGO phase boundary. It is here estimated through the surface selfdiffusivity of Ni[38] and gives a diffusion parameter of experimental data (1.38 J/m 2[39]) and values used in other studies (1.1J/m 2[40], 1J/m 2[6]). Finally, in order to set the oxidation kinetics, the driving force for a stress-free flat metal/oxide phase boundary, 0 I f , as well as the interface mobility, I M , have to be defined. No bibliographic or experimental data of these parameters have been found in the literature, and they have therefore been set from physical considerations. First, the product of the two values sets the initial metal/oxide boundary

If

  for the interface evolution, it thus provides references for the setting of this parameter. The first considered situation is that the effect of the curvature-related thermodynamic force alone on the phase-boundary propagation should be negligible on the initial interface geometry corresponding to the experimentally observed roughness. This allows for the uniform growth of the convoluted oxide scale when stresses are negligible. The initial waviness exhibiting a curvature radius of the order of 10µm, it requires: situation considers the minimum bound of the roughness. From observation, it is estimated to 0.1µm. Consequently, smaller morphologies should be naturally flattened, giving: has been chosen, which yields an interface mobility of

  Figure 7 a). Following the study by Busso and Qian [8], these times have been selected because they represent a lower and an upper bound, respectively, of the oxidation time leading to the critical conditions resulting in the TBC spallation upon cooling at rest temperature. Also shown on Figure 7 a) for comparison are the oxide/metal different times which would correspond to a uniform oxidation. Note that the undeformed configuration is presented in order to point out the oxidation kinetics contribution to the oxide scale development, eliminating the effect of the large phase transformation volume eigenstrain in the true oxide thickness growth. During the first 100 hours of oxidation, the phase propagation has remained quite uniform. But during the following 235 hours, local oxidation kinetics have significantly differed leading to the development of significant non-uniformities in the oxidation layer thickness. The local history of oxidation layer thickness at different locations along the phase boundary is presented on Figure 7 b), along with the average value for comparison. The points A, B and C correspond to the peak, the mid-slope and the valley, respectively, of the oxide-metal phase boundary as illustrated on Figure 7 a). The average oxidized zone thickness follows the empirical kinetics law, giving 0.8 µm after 100 hours and 1.8µm after 335 hours, but dissimilar local developments have been induced by the varying stress field along the inner interface. The highest oxidation kinetics is found to be at the peak and the slowest at valley, the mid-slope presenting a slightly above-average value. The main deviation is clearly located at the valley, where the oxidation is almost stopped after 300 hours. The corresponding relative oxide zone thickness variations with respect to the average value at the three locations after 335 hours are +8% (A), +2%(B) and -29% (C), respectively.

  and the change in elastic energy of the transformed volume over the process, , are presented on Figure 8 a). The data are averages of the simulation steps contributions weighted by the average step propagation increments over the time Page 25 of 80 http://mc.manuscriptcentral.com/pm-Authors ranges of interest in order to reflect the influences histories. It can be inferred that the local oxidation kinetics variations are essentially driven by the phase transformation work, which exhibits a much higher variation amplitude along the interface than the elastic energy change. This situation is due to the large phase transformation eigenstrain which leads to plastic accommodation of the components in the constrained directions, resulting in a similar dominant in-plane compressive state along the interface. However, the undulated geometry induces a varying situation along the free direction, which is the one of main development for the volume eigenstrain, resulting in high accommodation work variations along the phase boundary. The traction and the normal component of the deformation gradient jump along the phase boundary are presented on Figure 8 b), making explicit the main contributions to the oxidation work variations. The deformation normal to the phase boundary is quite uniform, around 1 0.28 I Π -= ,

  if little non-uniformity in the development can be observed at the term of the period.

  Figure9show that an increase by an order of magnitude of the interfacial energy results in a

( 2 )( 3 )

 23 barrier coating system is presented. It is based on a mechanism of stress-affected oxidation kinetics at the metal/oxide sharp phase boundary, for which a specific formulation based on continuum thermodynamics has been derived. A continuum mechanics/diffusion framework is presented to complete the model, including the volume eigenstrain generation from the main oxidation at the inner interface, the time-dependent and time-independent accommodation behaviors of the different materials, as well as the specific mass diffusion along the propagating phase boundary. A specific numerical simulation tool has been developed in order to solve for the Authors complex formulation while allowing the morphological evolution of the oxide-metal phase boundary. The study reveals that:(1) A stress-affected local oxidation kinetics assisted by a large mass redistribution through diffusion along the phase boundary might induce the development of an oxide scale of significantly non-uniform morphology, a situation likely to dramatically reduce the TBC system lifetime. The main influence on the oxidation kinetics local variation comes from the phase transformation work accompanying the propagation of the phase boundary under a traction profile related to the geometry. This work is induced by the accommodation of the large volume eigenstrain associated with the metallic phase oxidation, which develops mainly along the normal direction to the interface. The interfacial energy is a key property for the attenuation of the development of a rough interface, competing with geometry-induced traction variations. Therefore increasing the interfacial energy and/or the stress-induced flow rate near the phase boundary at the high functioning temperature would significantly reduce the development of non-uniformities in the TGO thickness.Assessment of the presented mechanism of stress-affected local oxide growth, its main influences and the provided recommendations should be achieved next to validate the model and the optimization directions. Experimental testing of the interfacial energy influence would be difficult due to the lack of data for different bond coat/Alumina interfaces and the likely variations induced by minor alloying elements, additives and impurities not only on this property but on the whole oxide growth process. However, the influence of the traction magnitude should be experimentally accessible, at least qualitatively, through comparison of the developed non-thickness for different specific experiments. A direct measurement of the traction along the metal/oxide phase boundary is not possible with today's techniques which detect an average stress level within the resolution zone, and therefore only allow recording the highest local stress which is the in-plane compression in thermally grown oxide scales. But the traction magnitude can be affected by experimental means while keeping the same material system. For instance it would be interesting to perform two similar isothermal oxidation experiments, one with the YSZ coating and another one with the coating removed, since the relaxation of the constraint imposed by the top layer should greatly affect the out-of-plane stress magnitude. Note that for a meaningful comparison the oxide growth kinetics should be similar in magnitude and development direction. Therefore the temperature gradient within the YSZ coating should be accounted for such that the oxide scales grow at similar temperatures in the two experiments. Another solution could be to directly apply a lateral stress to the thermal barrier coating system during an isothermal oxidation experiment. However, the imposed stress should be of the order of the in-plane growth stress within the TGO to be significant. At this value, high plastic deformation of the soft metallic bond coat and/or damage of the YSZ coating are likely, which might prevent from meaningful comparisons between experiments at different stress levels. Finally, a last track would be to perform isothermal oxidation experiments with bond coats of same composition but different grain sizes. This should affect the flow rate near the BC/TGO interface and as a consequence the local traction magnitude. In order to assess the influence, the initial TGO waviness should be similar for the different experiments, and the stress relaxation variation with grain size should be previously quantified and be significant.

F

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. Thermodynamic forces on the domain boundary and the interface.

Figure 3 .

 3 Figure 3. Infinitesimal time evolution of the system.

Figure 4 .

 4 Figure 4. Chart flow of the resolution scheme.

Figure 5 .

 5 Figure 5. Schematic of the phase field definition.

Figure 6 .

 6 Figure 6. Schematic of the TBC system model (not to scale).

Figure 7 .

 7 Figure 7. a) Oxidized bond coat layer after 100 and 335 hours of isothermal oxidation. The dashed lines show the interfaces corresponding to a uniform oxidation. b) Oxidation kinetics at different locations.

Figure 8 .

 8 Figure 8. Average a) stress contributions and b) tractions and normal deformation gradients jump over the indicated time period vs normalized coordinate along the oxide-metal boundary.

Figure 9 .

 9 Figure 9. Evolution with oxidation time of the ratio of the maximum local scale thickness to the minimum value. The results are presented for different oxide-metal interface energies and diffusivities.

  In thermal barrier coating (TBC) systems, an oxide layer develops at high temperature below the ceramic coating, leading at long term to the mechanical failure of the structure upon cooling. This study investigates a mechanism of stress-affected oxidation likely to induce the growth of a non-uniform oxide scale detrimental to the TBC lifetime. A continuum thermodynamics formulation is derived accounting for the influence of the stress and strain situation at the sharp metal/oxide phase boundary on the local oxidation kinetics. It specially includes the contributions of the large volumetric strain and the mass consumption associated with metal oxidation. A continuum mechanics/mass diffusion framework is used along with the developed formulation for the interface evolution to study the growth of an oxide layer coupled with local stress development. The implementation of the model has required the development of a specific simulation tool, based on a finite element method completed with an external routine for the phase boundary propagation. Results on an electron-beam physical vapor deposited (EB-PVD) TBC case are presented. The processes resulting in a non-uniform oxide scale growth are analyzed and the main influences are discussed.Keywords: modeling; interfacial thermodynamics; high-temperature oxidation; numerical simulation; thermal barrier coatings Thermal barrier coatings are used in aircraft and industrial gas-turbine engine to allow a high functioning temperature while maintaining sustainable conditions for the metallic structure.However, at high temperature, an oxide scale (commonly named thermally grown oxide, TGO) develops between the coating, usually made of Yttria-stabilized Zirconia (YSZ), and the metallic bond coat (BC), limiting the lifetime of the coating and as a consequence the long term reliability of the structure. The oxide scale, mainly constituted of Alumina oxide, presents significantly different thermal expansion and creep properties than the surrounding YSZ and BC layers. This induces high stresses and deformations upon cooling at room temperature, leading to crack nucleation and eventually to the spallation of the top coat. Thus, the growth of the TGO during high temperature operations is the most important phenomenon responsible for the TBC intrinsic failure[1, 2]. However, if the global mechanisms driving the oxide scale development are well known, the factors for varying kinetics and the development of a rough and undulated scale are poorly understood.Several studies have investigated the development of undulations at the TGO/BC interface, pointing out the key roles of oxide growth stresses and interface geometry[2][3][4][5]. These works focus on plastic and creep deformation during cooling and heating cycles, or stress-driven material diffusion within the BC[6, 7], but all consider a uniform growth of the TGO. To our knowledge, no study has been performed considering a local effect of stress on the oxidation kinetics at the BC/TGO interface, which could result in the non-uniform growth of the oxide scale, a factor which has been identified as strongly detrimental to the TBC lifetime[2]. It would result in additional stress concentrations upon cooling and might lead to large local scale deformations of the oxide layer, accelerating the development of cracks and spallations. Even

  morphology variations might become significant under the effect of thermal cycling. Therefore, it is of primary interest to understand the mechanisms leading to a nonuniform oxide growth, quantify the contributions and identify the most influential physical parameters. This should allow for the optimization of the TBC material system design in terms of composition and fabrication process, and also help predicting the location of failure.A mechanism is proposed here relying on an influence of the stress/strain state on each side of the sharp oxide/metal interface on the local oxidation kinetics. A specific formulation is derived for this effect considering the thermodynamic force driving the sharp oxide/metal phase boundary propagation, and specially including the contributions of the large volumetric strain and the mass consumption associated with metal oxidation. A continuum mechanics/diffusion framework completes the model, accounting for the diffusion processes through the material layers and along the interface, globally driving the oxidation, as well as the large phase transformation eigenstrain development and its accommodation. The numerical simulation of this complex formulation required the development of an original numerical scheme coupling a finite element method for the calculation of the stress and diffusion field, with an external routine for the resolution of the phase boundary local composition and propagation. The non-uniform oxide scale growth and its coupling with stress development are simulated and investigated in the case of an Alumina layer growing between an electron-beam physical vapor deposited (EB-PVD) Yttria-stabilized Zirconia thermal barrier coating and a Ni-Pt-Al bond coat at 1100°C.

  Authorsconsumption during the oxidation reaction engenders the creation of vacancies at the oxide boundary, which diffuse counter to oxygen ions towards the free surface exposed to an oxygenrich atmosphere. The oxidation chemical reaction being much faster than the species diffusion, the inner oxidation process is limited by the local availability of oxygen at the interface.Consequently, coupled quasi-equilibrium concentrations of oxygen and vacancies rapidly establish at the interface and are maintained during the phase propagation. The mechanism proposed here for a stress influence on the local oxidation kinetics relies on a deviation from the stress-free quasi-equilibrium concentration induced by the mechanical state at the propagating interface. Since the gradient of concentration drives species diffusion, the local incoming flux of oxygen at the phase boundary is in turn affected and eventually the oxidation kinetics (oxide phase propagation) is locally modified.The principle of a thermodynamic effect induced by the stress development accompanying metal oxidation at high temperature on the process kinetics has been reviewed by Evans[9]. He demonstrated through a simple direct influence on interface vacancy concentration that the accommodation work in response of the metal oxidation-associated volume expansion might significantly alter and even stop the uniform growth of a Zirconia film. Here this principle is extended and a complete local treatment is proposed, allowing for the simulation of a complex non-uniform phenomenon. A continuum mechanics description of a dissipative sharp interface propagation is derived to include the influence of the stress state and accommodation work on the phase boundary propagation kinetics. This thermodynamic approach, extensively developed over the last decades (see[10] for general treatment and references), allows for a direct description of the stress and composition influences on the local interface motion kinetics. Very few works have considered a mechanical coupling with phase boundary propagation during high-temperature . Garruchet et al. derived an interfacial thermodynamic model to investigate an

Figure 1 .

 1 Figure 1. Since continuum media are considered, any volume change is the result of a mechanical

  of the velocity and the deformation gradient can be related through the Hadamard compatibility condition:

  "trapping" within the formed material and "release" from the dissolved one.

  continuity of diffusion potentials at the interface, Equation (9), along with the local mass balance, Equation (6), yields:

C

  and Ω + , due to the motion of the external boundary on one hand, and of the interface Γ on the other hand. Only this last term is related to the phase boundary propagation and is then of interest here. The total free energy density W includes contributions from the stored elastic energy density, e W , and the chemical free energy density, stand for the chemical potential and concentration, respectively, of any species i present in the phase considered. Note that the Einstein notation is again used so that the double indices indicate a summation over the species. At the interface, a specific thermodynamic force induced by the gradient of local curvature along the interface applies, tending to flatten the boundary in order to reduce its surface and thus the global interface energy[16]. This effect can be directly introduced in the time evolution of the free energy induced by the phase boundary propagation, which then overall can be expressed by:

  and the change in elastic energy of the transformed volume over the process,

5. Numerical simulation 5 . 1 .

 51 Numerical frameworkImplementation of a general numerical model for a sharp propagating phase boundary coupled with an evolving domain is very challenging. The solution can be determined by coupling two

  developed and are the subjects of extensive research for the treatment of phase boundary

Figure 4 .

 4 Figure 4. Two steps are sequentially performed for each time increment.

  the normal distance to the phase boundary line allows building a field evolving linearly over the finite thickness, as defined in the figure. The updated field is not directly enforced, but will result from a progressive time-dependent variation in order to ensure a smooth phase transformation, as defined in Figure5. Times for the onset of the local phase field evolution as well as the corresponding transformation kinetics are constructed based on phase boundary propagation increments. Additionnaly, a node-based interface is defined through identification of a path formed by FEM mesh nodes following the true phase boundary. It is used to enforce the oxygen concentration boundary conditions for the next diffusion analysis, interpolated from the values computed at the interface tracking points.In a second step, a transient coupled stress/diffusion FEM analysis is then performed, considering the displaced phase boundary as being fixed, to resolve for the time evolution of the diffusion and mechanical fields. The oxygen concentration is obtained using Fick's first law for the diffusion fluxes and a classical mass conservation statement. The material displacement, and the associated stress and strain fields, are solved through continuum mechanics, as described in section 4. The local material diffusivity only depends on the element location with respect to the node-based interface (metal or oxide side), while mechanical properties vary linearly, as defined by the phase field, from the untransformed to the transformed materials values. Therefore the interface is only diffuse for the mechanical fields' resolution.

  interface displacements (due to the sequential evolutions of the interface and the domain fields), it allows for the treatment of complex phase boundary composition and propagation evolutions dependent on discontinuous field values and coupled with specific interface processes such as fast mass diffusion and volume eigenstrain generation. It thus meets the specific and stringent requirements of the oxidation formulation. The FEM computations have been performed using the commercial code ABAQUS, while specific programs developed under MATLAB have been used for the pre/post-calculations.

  , as it is known, transparent for oxygen diffusion. Symmetry boundary conditions are applied on the lateral sides of the domain for diffusion (zero flux normal to the boundary). To the author's knowledge, no values are available in the literature for the vacancy diffusivity along the BC/TGO phase boundary. It is here estimated through the surface selfdiffusivity of Ni[38] and gives a diffusion parameter of experimental data (1.38 J/m 2[39]) and values used in other studies (1.1J/m 2[40], 1J/m 2[6]). Finally, in order to set the oxidation kinetics, the driving force for a stress-free flat metal/oxide phase boundary, 0 I f , as well as the interface mobility, I M , have to be defined. No bibliographic or experimental data of these parameters have been found in the literature, and they have therefore been set from physical considerations. First, the product of the two values sets the initial metal/oxide boundary stress and curvature contributions are neglected.

If

  for the interface evolution, it thus provides references for the setting of this parameter. The first considered situation is that the effect of the curvature-related thermodynamic force alone on the phase-boundary propagation should be negligible on the initial interface geometry corresponding to the experimentally observed roughness. This allows for the uniform growth of the convoluted oxide scale when stresses are negligible. The initial waviness exhibiting a curvature radius of the order of 10µm, it requires: bound of the roughness. From observation, it is estimated to 0.1µm. Consequently, smaller morphologies should be naturally flattened, giving:

Figure 7 a

 7 ). Following the study by Busso and Qian[8], these times have been selected because they represent a lower and an upper bound, respectively, of the oxidation time leading to the critical conditions resulting in the TBC spallation upon cooling at rest temperature. Also shown on Figure7a) for comparison are the oxide/metal different times which would correspond to a uniform oxidation. Note that the undeformed configuration is presented in order to point out the oxidation kinetics contribution to the oxide scale development, eliminating the effect of the large phase transformation volume eigenstrain in the true oxide thickness growth. During the first 100 hours of oxidation, the phase propagation has remained quite uniform. But during the following 235 hours, local oxidation kinetics have significantly differed leading to the development of significant non-uniformities in the oxidation layer thickness. The local history of oxidation layer thickness at different locations along the phase boundary is presented on Figure 7 b), along with the average value for comparison. The points A, B and C correspond to the peak, the mid-slope and the valley, respectively, of the oxide-metal phase boundary as illustrated on Figure 7 a). The average oxidized zone thickness follows the empirical kinetics law, giving 0.8 µm after 100 hours and 1.8µm after 335 hours, but dissimilar local developments have been induced by the varying stress field along the inner interface. The highest oxidation kinetics is found to be at the peak and the slowest at valley, the mid-slope presenting a slightly above-average value. The main deviation is clearly located at the valley, where the oxidation is almost stopped after 300 hours. The corresponding relative oxide zone thickness variations with respect to the average value at the three locations after 335 hours are +8% (A), +2%(B) and -29% (C), respectively.

  and the change in elastic energy of the transformed volume over the process, on Figure 8 a). The data are averages of the simulation steps contributions weighted by the average step propagation increments over the time Authors ranges of interest in order to reflect the influences histories. It can be inferred that the local oxidation kinetics variations are essentially driven by the phase transformation work, which exhibits a much higher variation amplitude along the interface than the elastic energy change. This situation is due to the large phase transformation eigenstrain which leads to plastic accommodation of the components in the constrained directions, resulting in a similar dominant in-plane compressive state along the interface. However, the undulated geometry induces a varying situation along the free direction, which is the one of main development for the volume eigenstrain, resulting in high accommodation work variations along the phase boundary. The traction and the normal component of the deformation gradient jump along the phase boundary are presented on Figure 8 b), making explicit the main contributions to the oxidation work variations. The deformation normal to the phase boundary is quite uniform, around 1

  if little non-uniformity in the development can be observed at the term of the period.

Figure 9

 9 Figure9show that an increase by an order of magnitude of the interfacial energy results in a

( 2 )( 3 )

 23 barrier coating system is presented. It is based on a mechanism of stress-affected oxidation kinetics at the metal/oxide sharp phase boundary, for which a specific formulation based on continuum thermodynamics has been derived. A continuum mechanics/diffusion framework is presented to complete the model, including the volume eigenstrain generation from the main oxidation at the inner interface, the time-dependent and time-independent accommodation behaviors of the different materials, as well as the specific mass diffusion along the propagating phase boundary. A specific numerical simulation tool has been developed in order to solve for the Authors complex formulation while allowing the morphological evolution of the oxide-metal phase boundary. The study reveals that:(1) A stress-affected local oxidation kinetics assisted by a large mass redistribution through diffusion along the phase boundary might induce the development of an oxide scale of significantly non-uniform morphology, a situation likely to dramatically reduce the TBC system lifetime. The main influence on the oxidation kinetics local variation comes from the phase transformation work accompanying the propagation of the phase boundary under a traction profile related to the geometry. This work is induced by the accommodation of the large volume eigenstrain associated with the metallic phase oxidation, which develops mainly along the normal direction to the interface. The interfacial energy is a key property for the attenuation of the development of a rough interface, competing with geometry-induced traction variations. Therefore increasing the interfacial energy and/or the stress-induced flow rate near the phase boundary at the high functioning temperature would significantly reduce the development of non-uniformities in the TGO thickness.Assessment of the presented mechanism of stress-affected local oxide growth, its main influences and the provided recommendations should be achieved next to validate the model and the optimization directions. Experimental testing of the interfacial energy influence would be difficult due to the lack of data for different bond coat/Alumina interfaces and the likely variations induced by minor alloying elements, additives and impurities not only on this property but on the whole oxide growth process. However, the influence of the traction magnitude should be experimentally accessible, at least qualitatively, through comparison of the developed non-thickness for different specific experiments. A direct measurement of the traction along the metal/oxide phase boundary is not possible with today's techniques which detect an average stress level within the resolution zone, and therefore only allow recording the highest local stress which is the in-plane compression in thermally grown oxide scales. But the traction magnitude can be affected by experimental means while keeping the same material system. For instance it would be interesting to perform two similar isothermal oxidation experiments, one with the YSZ coating and another one with the coating removed, since the relaxation of the constraint imposed by the top layer should greatly affect the out-of-plane stress magnitude. Note that for a meaningful comparison the oxide growth kinetics should be similar in magnitude and development direction. Therefore the temperature gradient within the YSZ coating should be accounted for such that the oxide scales grow at similar temperatures in the two experiments. Another solution could be to directly apply a lateral stress to the thermal barrier coating system during an isothermal oxidation experiment. However, the imposed stress should be of the order of the in-plane growth stress within the TGO to be significant. At this value, high plastic deformation of the soft metallic bond coat and/or damage of the YSZ coating are likely, which might prevent from meaningful comparisons between experiments at different stress levels. Finally, a last track would be to perform isothermal oxidation experiments with bond coats of same composition but different grain sizes. This should affect the flow rate near the BC/TGO interface and as a consequence the local traction magnitude. In order to assess the influence, the initial TGO waviness should be similar for the different experiments, and the stress relaxation variation with grain size should be previously quantified and be significant.

F

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. Thermodynamic forces on the domain boundary and the interface.

Figure 3 .

 3 Figure 3. Infinitesimal time evolution of the system.

Figure 4 .

 4 Figure 4. Chart flow of the resolution scheme.

Figure 5 .

 5 Figure 5. Schematic of the phase field definition.

Figure 6 .

 6 Figure 6. Schematic of the TBC system model (not to scale).

Figure 7 .

 7 Figure 7. a) Oxidized bond coat layer after 100 and 335 hours of isothermal oxidation. The dashed lines show the interfaces corresponding to a uniform oxidation. b) Oxidation kinetics at different locations.

Figure 8 .

 8 Figure 8. Average a) stress contributions and b) tractions and normal deformation gradients jump over the indicated time period vs normalized coordinate along the oxide-metal boundary.

Figure 9 .

 9 Figure 9. Evolution with oxidation time of the ratio of the maximum local scale thickness to the minimum value. The results are presented for different oxide-metal interface energies and diffusivities.

Figure 10 .

 10 Figure 10. Curvature radius vs traction amplitude.

5.1. Numerical framework

	Philosophical Magazine	17/34
	tasks: calculating the interface displacement following a given kinetics and solving for the
	continuous fields on each side of the discontinuity. Strong couplings appear when conservations
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	Implementation of a general numerical model for a sharp propagating phase boundary coupled
	with an evolving domain is very challenging. The solution can be determined by coupling two

  Authorsdeformations practically imposed by the thick substrate, displacement boundary conditions are effectively applied in these directions. At the oxidation temperature, a value typical of creep rates exhibited by single crystal superalloys in service is used,

	22/34
	ε	11	ε	33	9 1 10 s --
	F o
	r
	P
	e e r
	R e v i e w
			O n
					l y

Table 1 .

 1 Mechanical properties of the bond coat, the TGO and YSZ at 1100ºC.

	32/34	Authors
	o
	r
	P
	e e r
		R e v i e w
		O n
		l y

Page 34 of 80 http

	Philosophical Magazine & Philosophical Magazine Letters
	Philosophical Magazine,
	Vol. X, No. X, Month 200X, 000-000
	F o r P e e r R e v i e w O n l y 277x160mm (96 x 96 DPI) F o r P e e r R e v i e w O n l y 120x69mm (96 x 96 DPI) F o r P e e r R e v i e w O n l y 155x70mm (96 x 96 DPI) F o r P e e r R e v i e w O n l y 180x187mm (96 x 96 DPI) F o r P e e r R e v i e w O n l y 257x180mm (72 x 72 DPI) F o r P e e r R e v i e w O n l y F o r P e e r R e v i e w O n l y 79x79mm (150 x 150 DPI) F o r P e e r R e v i e w O n l y 79x79mm (150 x 150 DPI) F o r P e e r R e v i e w O n l y 79x79mm (150 x 150 DPI) F o r P e e r R e v i e w O n l y 79x79mm (150 x 150 DPI) F o r P e e r R e v i e w O n l y 79x79mm (150 x 150 DPI) F o r F o r P e P e e r e r R e R v e i e v i e w 79x79mm (150 x 150 DPI) w O n l y O n l
	91x132mm (120 x 120 DPI)	y

://mc.manuscriptcentral.com/pm-pml *Corresponding author. Email: mcherkaoui@me.gatech.edu

Stress

influence on high temperature oxide scale growth: modeling and investigation on a thermal barrier coating system.

	A. Saillard a , M. Cherkaoui a *, L. Capolungo a and E.P. Busso b
	a UMI 2958, Georgia Tech -CNRS, Metz, France
	G. W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta,
	GA 30332-0405, USA
	LPMM-CNRS Université Paul Verlaine, Ile du Saulcy, 57045 Metz, Cedex 1, France
	b Centre des Matériaux, Mines ParisTech, -UMR CNRS 7633, Evry, France

  a linear kinetic relation is assumed:

	Philosophical Magazine	13/34
	n	I I
	F o r
	P
	e
	e r
	R
	e
		v i e
		w
		O n l
		y

Table 1 .

 1 Mechanical properties of the bond coat, the TGO and YSZ at 1100ºC.

	32/34	Authors
	o r
	P
	e
	e r
		R
		e
		v i e
		w
		O n l
		y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters