An Oxide Growth-Coupled Viscoplasticity Model and Its Application to Interfacial Stress Analysis near an Air Hole within a Thermal Barrier Coating

: Strength assessment for thermal barrier coatings (TBCs) is vital in the safety design of hot-section components in engines. However, several crucial factors, including thermally grown oxide (TGO) growth and creep–plasticity interaction, have been less considered in thermo-mechanical analyses for TBCs near air holes. In this study, a unified viscoplastic constitutive model incorporating TGO growth is developed and integrated into a finite element framework. The model considers multiple factors, including TGO growth, creep–plasticity interaction, interface undulation, and temperature gradient. Additionally, an analytical solution for the non-uniform temperature field of a TBC is derived. The model is then applied to calculate interfacial stresses and accumulated strain energies in the TBC near an air hole, which promote interface debonding. The obtained results can be utilized to investigate the mechanisms of hole edge delamination in TBCs, considering the combined effects of multiple complex factors. A competition for the potential failure initiation location is revealed between the first oxide layer and the evolving TGO/bond coat interface. The developed viscoplasticity model demonstrates effective capability in modelling a range of dynamic behaviors that collectively contribute to hole edge delamination failure.


Introduction
Thermal barrier coatings (TBCs) are widely utilized in hot-section components, such as combustors, turbine blades, and vanes, due to their effective thermal insulation properties [1][2][3][4].The TBC system consists of a multilayer structure, typically including the ceramic top coat (TC), thermally grown oxides (TGOs), metallic bond coat (BC), and superalloy substrate.Various types of air holes are strategically engineered in both hot-section components and applied TBCs to ensure stable combustion and efficient cooling in gas turbines.Secondary air holes, comprising primary holes and dilution holes, are machined into engine combustor liners, while film-cooling holes are present in both combustor liners and turbine blades.These holes have diameters ranging from hundreds of microns to tens of millimeters.
Near the air hole, complex stress and strain fields emerge due to thermal expansion misfit and the free-edge effect.Experimental observations have shown that under cyclic stresses, interfacial cracks can initiate from the hole edge and propagate towards the inboard region in multilayer structures, exacerbating the structural integrity of both TBCs and the underlying components [5].Therefore, conducting strength assessments for the TBC near an air hole is crucial in the safety design of TBC systems and hot-section components.Hsueh et al. [6] identified that peeling stress perpendicular to the interface and shear Coatings 2024, 14, 362 2 of 19 stress contribute to hole edge delamination and spallation in TBCs near the air holes.Liu et al. [7] developed a finite element (FE) model for a guide vane coated with TBCs and explored the impact of film hole geometries on stress distributions.Chiu et al. [8] examined temperature and stress distributions in TBCs with film-cooling holes and calciummagnesium-aluminum-silicate (CMAS) infiltration using a thermo-fluid-solid coupling FE model.However, the above studies simplified or overlooked the multilayered geometry of TBCs, particularly neglecting the presence of TGO, consequently failing to investigate crucial peeling and shear stresses near the TGO and hole edge.Jiang et al. [9,10] performed cyclic thermal tests and finite element (FE) simulations on TBCs with air holes.They observed a rise in interfacial stresses towards the hole, resulting in a notable decrease in the lifespan of TBCs.Cai et al. [11] explored the impact of TGO interface morphology on stress distribution near an inclined film-cooling hole using a two-dimensional FE model.The results indicated that the maximum tensile stress in the TC initially increased and then decreased with increasing amplitude.Importantly, the cyclic thermal load experienced by TBCs encompasses a wide temperature range, fluctuating from over 1000 • C to room temperature.The combination of severe load conditions with stress concentrations near the hole edge and TGO leads to plastic and creep deformations, which consequently interact with each other.However, in most studies, the materials of TBCs are assumed to be either elastic or ideally elastoplastic.This assumption neglects the nonlinearity of the stress-strain relationship and the strain rate effect caused by creep-plasticity interaction, potentially leading to remarkable errors in stress evaluation for TBCs with air holes under cyclic thermo-mechanical loading [12].
To accurately calculate stress and investigate the creep (viscous) and plastic behaviors of materials at elevated temperatures, two groups of viscoplasticity modelling frameworks have been developed over the past few decades: non-unified and unified constitutive modellings.In the non-unified viscoplastic model, creep and plastic strains are separately quantified as rate-dependent and -independent variables [13].Therefore, predictions concerning ratcheting and the creep-plasticity interaction may not align with experimental results [14].In contrast, the unified viscoplasticity model considers the viscoplastic strain as a unique rate-dependent quantity for the creep-plasticity interaction [15].The present authors [16][17][18][19] developed a unified Chaboche-Lemaitre constitutive model integrating a power flow rule with nonlinear anisothermal evolution of isotropic and kinematic hardenings.Subsequently, in study [20], they implemented this model into a FE modelling framework for TBCs near air holes subjected to cyclic thermo-mechanical loading.Their findings revealed that the maximum interfacial peeling and shear stresses near the hole edge increase with the hole radius, raising the likelihood of edge delamination.However, it is noteworthy that the constitutive model and numerical modelling framework developed in [20] did not account for dynamic TGO growth and interface undulations, which play crucial roles in stress and strain energy distributions and evolutions in TBCs [1].
The bond coat in TBCs, typically composed of MCrAlY, is rich in aluminum.During prolonged exposure to high-temperature environments, oxygen anions diffuse inward through the TC and react with aluminum cations from the BC, resulting in the formation of α-Al 2 O 3 , the primary component of TGO [21].Previous research has extensively simulated TGO growth, but few studies consider air holes and the creep-plasticity interaction [22][23][24][25][26][27][28].Two predominant approaches for the mechanical equivalent method to model TGO growth are discussed: element volumetric swelling [24][25][26] and element conversion with the imposition of growth strains [23,27,28].The former approach involves imposing swelling strains on a layer of elements in the initial TGO, thus limiting the growth thickness to usually not exceed the initial TGO thickness [24].The latter approach allows for largescale TGO thickening and provides a better description of the actual behavior of TGO growth [23].The present work develops a unified viscoplastic constitutive model that couples the dynamic growth of TGO and derives an analytical solution for the non-uniform temperature field of TBCs.Integrating the constitutive model into a FE framework enables the effective investigation of interfacial stresses, strain energies, and hole edge delamination mechanisms in the TBC near an air hole.This study considers the combined effects of TGO growth, creep-plasticity interaction, interface undulation, temperature gradient, sintering, etc.This paper is structured as follows: In Section 2, we elaborate on the development of an oxide growth-coupled unified viscoplastic constitutive model.Section 3 presents the derivation of an analytical solution for the temperature field of the TBC during the high-temperature holding phase.In Section 4, we explore the application of the developed constitutive model in TBCs near air holes, establishing a comprehensive viscoplasticity FE modelling framework.Here, we analyze the interfacial stresses and strain energies and examine the hole edge delamination mechanism under various factors.Finally, the paper concludes with a summary of key findings and conclusions in Section 5.

The Unified Viscoplastic Constitutive Model Incorporating TGO Growth
A unified viscoplastic constitutive model is presented to incorporate TGO growth and the creep-plasticity interaction in TBCs.This model is developed based on the advanced unified Chaboche-Lemaitre constitutive model for TBC systems proposed by the present authors in [20].

The Unified Viscoplasticity Model
In the unified viscoplasticity model, the viscoplastic strain is considered a single ratedependent quantity, which incorporates the interaction between creep and plasticity.The strain resulting from TGO growth can be linearly superimposed onto the total strain [24].Hence, the total strain ε ij is expressed by where ε e ij , ε vp ij , and ε TGO ij are the elastic, viscoplastic, and TGO growth strain tensors, respectively.The quantity ε e ij is obtained by where D ijkl and σ kl are the stiffness and corresponding stress tensors, respectively.The viscoplastic strain rate .ε vp ij is expressed by a power flow law, that is, .
p is the viscoplastic multiplier rate and defined as .
. ε vp ij , and N ij indicates the viscoplastic flow direction.Z(T) and n(T) are temperature-dependent material parameters representing the viscoplastic resistance and viscous exponent functions, respectively.S ij is the deviatoric stress tensor of σ ij , which is given by where δ ij is the Kronecker delta.In addition, X ij in Equation (3) denotes the kinematic hardening, i.e., back stress tensor.In this constitutive model, the anisothermal evolution of X ij using the Chaboche decomposition [15] is given as Coatings 2024, 14, 362 and the rate of X ij is given as .
T is the rate of temperature.Furthermore, the von Mises yield function f in Equation ( 3) is given as where R(T, p(T)) is the isotropic hardening parameter, p(T) is the viscoplastic multiplier, and σ iy (T) is the yield stress.The anisothermal evolution of R(T, p(T)) is expressed by .

R(T, p(
where H(T) represents the linear drag stress evolution, υ(T) and µ(T) are temperaturedependent material parameters.The terms H ,T (T), υ ,T (T) and µ ,T (T) are the derivatives of H(T), υ(T) and µ(T) concerning temperature.The second invariant J 2 in Equation ( 8) also depends on X ij and σ ij , which is expressed by

TGO Growth Modelling
In high-temperature environments, oxygen anions diffuse inward through the TC and react with Al cations from the BC at the TGO/BC interface and at the grain boundaries within TGO.This process leads to through-thickness and lateral growth strains, respectively, while the thickness of TGO increases [1,3,24].Experimental observations indicate that the growth of TGO thickness, denoted as h, generally conforms to the parabolic law [1,25]: where h 0 is the initial TGO thickness, A is the TGO growth coefficient, E A is the TGO growth activation energy, k B is the Boltzmann constant, t oxid is the oxidation time, and n oxid is the oxidation exponent.As mentioned above, the TGO growth strain ε TGO ij in Equation ( 1) can be divided into two components, ε ⊥ and ε ∥ , representing the strains normal and parallel to the TGO/BC interface, respectively.The corresponding growth strain rates can be calculated using Equations ( 12) and (13) [23,29]: and where C ox and D ox are the through-thickness and lateral growth strain rate constants, respectively.
During the high-temperature holding phase, TGO undergoes growth.In the constitutive model, the material properties of the layer adjacent to the TGO transition gradually from those of the BC to TGO properties throughout the holding phase.This transition occurs layer by layer, with one layer transforming in each cycle.Simultaneously, both normal and parallel TGO growth strains are applied to the initial and converted TGO layers during the holding phases.Specifically, normal growth strains are applied during the holding phase when the layer is transforming to TGO.Subsequently, parallel growth strains are applied in the remaining holding phases after the layer has completed the conversion.

UMAT Implementation of the Constitutive Model
The developed unified viscoplastic constitutive model, which incorporates TGO growth, was implemented in a user-defined material subroutine (UMAT) for subsequent integration into a viscoplasticity FE modelling framework.The numerical implementation of this constitutive model employed the return mapping algorithm [30].The process consists of two steps: elastic prediction and viscoplastic correction.In the first step, the total strain is assumed to be entirely elastic, and therefore, the viscoplastic strain and other internal variables are not considered.In the second step, known as viscoplastic correction, the total strain is held fixed while the internal variables evolve.This step involves numerical integrations of the rate-form equation of the viscoplastic strain, typically performed using the Newton-Raphson algorithm.
By employing Hooke's law and the unified viscoplasticity framework, the following expressions can be derived: Φ := where ∆ε vp denotes the increment in the viscoplastic strain, G is the shear modulus, ∆t is the time increment, and σ tr , J 2 tr , N tr and Φ are the trial stress, trial second invariant, trial flow direction and viscoplastic residue, respectively.The nullity of the residual function in Equation ( 16) can be used to solve the accumulated plastic strain when f > 0. This is iteratively performed through the Newton-Raphson scheme.After rearrangement, the Newton increment in the accumulated plastic strain rate can be expressed as follows: where ,T .The Newton-Raphson increment in the accumulated plastic strain, δp, can be determined and updated by δp = • p∆t.When the norm of Φ is less than the tolerance, the iterative scheme stops and the time increment in the viscoplastic strain, ∆ε vp , is updated by ∆ε vp = N • p∆t, while the total stress σ is corrected by Equation ( 14).

Analytical Solution for the Non-Uniform Temperature Field
From Equations ( 11)- (13), it is evident that the TGO growth rate and growth strains exhibit temperature dependence.This section derives the analytical solution for the nonuniform temperature field during the high-temperature holding phase, taking into account the thermal gradient.

Problem Description
Figure 1a illustrates the three-dimensional axisymmetric structure of a TBC system with a central air hole.The system consists of four layers: substrate, BC, initial TGO, and TC.The top and bottom surfaces represent the hot and cooling sides, respectively, resulting in a temperature gradient along the thickness direction.The TBC under investigation experiences axisymmetric thermo-mechanical loads.Hence, we establish an axisymmetric unit model using a cylindrical coordinate system, depicted in Figure 1b.Here, the z-axis aligns with the centerline of the air hole, while the r-axis lies in the plane of the bottom.The surface of each layer is represented by z i , where the subscript i ranges from 1 to 4. Failures in TBCs commonly occur at or near the undulations of the TGO layer [31].The model accounts for the irregular TC/TGO and TGO/BC interfaces by employing sinusoidal curves, as shown in Figure 1b.
exhibit temperature dependence.This section derives the analytical solution for the nonuniform temperature field during the high-temperature holding phase, taking into account the thermal gradient.

Problem Description
Figure 1a illustrates the three-dimensional axisymmetric structure of a TBC system with a central air hole.The system consists of four layers: substrate, BC, initial TGO, and TC.The top and bottom surfaces represent the hot and cooling sides, respectively, resulting in a temperature gradient along the thickness direction.The TBC under investigation experiences axisymmetric thermo-mechanical loads.Hence, we establish an axisymmetric unit model using a cylindrical coordinate system, depicted in Figure 1b.Here, the z-axis aligns with the centerline of the air hole, while the r-axis lies in the plane of the bottom.The surface of each layer is represented by  , where the subscript i ranges from 1 to 4. Failures in TBCs commonly occur at or near the undulations of the TGO layer [31].The model accounts for the irregular TC/TGO and TGO/BC interfaces by employing sinusoidal curves, as shown in Figure 1b.

Analytical Solution for the Non-Uniform Temperature Field
The temperature of the TBC remains constant during the high-temperature holding phase.Heat transfer within the TBC system along the thickness direction can be modeled as one-dimensional steady state.Therefore, the temperature in the ith layer of the TBC (see Figure 1b), denoted as  , satisfies where  is the conductivity coefficient for the ith layer.The corresponding heat flux  can be expressed by Next, by integrating Equation ( 18) twice, the general solution for  is obtained as where  and  are integration constants.
To identify  and  in Equation ( 20), the continuity conditions of the temperature and heat flux at the interfaces for the adjacent layers are employed.Since the amplitude of the interface roughness is merely 10 µm, it is reasonable to treat the TC/TGO and

Analytical Solution for the Non-Uniform Temperature Field
The temperature of the TBC remains constant during the high-temperature holding phase.Heat transfer within the TBC system along the thickness direction can be modeled as one-dimensional steady state.Therefore, the temperature in the ith layer of the TBC (see Figure 1b), denoted as T i , satisfies where k i is the conductivity coefficient for the ith layer.The corresponding heat flux q i can be expressed by Next, by integrating Equation ( 18) twice, the general solution for T i is obtained as where A i and B i are integration constants.
To identify A i and B i in Equation ( 20), the continuity conditions of the temperature and heat flux at the interfaces for the adjacent layers are employed.Since the amplitude of the interface roughness is merely 10 µm, it is reasonable to treat the TC/TGO and TGO/BC interfaces as flat surfaces when applying continuity conditions in the analytical derivation, which are and In addition, the temperature values on the top and bottom surfaces are also used as Coatings 2024, 14, 362 Substituting Equations ( 21)-( 24) into Equation ( 20) leads to A i and B i : and The accuracy of the derived analytical solution for the temperature distribution of the TBC system during the high-temperature holding phase will be verified in Section 4.4.1.

Application of the Constitutive Model to Stress Analysis near an Air Hole within TBCs
Complex stresses develop and evolve near the air hole within TBCs owing to the freeedge effect, dynamic TGO growth, and severe cyclic thermo-mechanical loading.Moreover, creep and plastic deformations exist and interact.To analyze stress distribution and hole edge delamination in TBCs near air holes, the developed constitutive model is implemented and integrated into a FE modelling framework.

Cyclic Thermo-Mechanical Loading
The TBC system surrounding an air hole underwent cyclic thermo-mechanical loading for ten cycles.Each cycle lasted 8400 s and comprised a 10-minute heating phase, a 2hour high-temperature holding phase, and a 10-minute cooling phase, as illustrated in Figure 2. In the heating phase, the temperatures of the TBC system gradually increased from an initial value of 25 • C to their maximum values during the high-temperature holding phase.Due to their effective thermal insulation properties, the TBCs exhibited a temperature gradient along their thickness direction.During the high-temperature holding phase, the temperatures on the top and bottom surfaces of the TBC system were 1200 • C and 1000 • C, respectively.In the cooling phase, the temperatures of the system gradually decreased from their maximum values to the final value of 25 • C. Thermal analyses determined the time-varying temperature distribution and thermal gradients during the heating and cooling phases.Apart from cyclic thermal loads, TBCs employed in hot-section components frequently encounter a pressure difference between their outer and inner surfaces.Acquiring a realistic time-dependent pressure load for TBCs in an aero-engine presents a challenge and it has often been either disregarded [32] or assumed to be constant [33].Nevertheless, the pressure difference escalates with the increase in thermal load, peaking during engine operation [34].Hence, in this study, the pressure difference ∆p(t) was presumed to be in sync with and linearly correlated to the temperature, as depicted in Figure 2. Within each cycle, the pressure difference increased from zero to its maximum value during the heating phase, remained constant throughout the holding phase, and subsequently decreased to zero during the cooling phase.A magnitude of 0.2 MPa was employed for the maximum pressure difference [34].
Coatings 2024, 14,362 was presumed to be in sync with and linearly correlated to the temperatur in Figure 2. Within each cycle, the pressure difference increased from zero to value during the heating phase, remained constant throughout the holdi subsequently decreased to zero during the cooling phase.A magnitude o employed for the maximum pressure difference [34].

Finite Element Model
The sequentially coupled thermo-mechanical analysis comprised two s a heat transfer simulation was performed to acquire the thermal input for t structural analysis and validate the derived analytical solution for the tem subsequently, the structural analysis was executed utilizing the develope model from Section 2. Both heat transfer and structural analyses employ with different settings, encompassing the input of material properties, elem ing, and boundary conditions.
An axisymmetric FE model of the TBC system surrounding an air structed based on the geometry depicted in Figure 1b.The materials comp strate, BC, initial TGO, and TC were Inconel 718, NiCoCrAlY, α-Al2O3, and stabilized ZrO2 (YSZ), respectively.The respective thicknesses were 1.2 m µm, and 300 µm, with an additional air hole radius of 1.5 mm.The unit m passed one and a half wavelengths, incorporating features of peaks, valleys as illustrated in Figure 3.The rough interface had a wavelength of 0.06 mm tude of 0.01 mm.The growth thickness of TGO in each cycle could be co Equation (11), utilizing the analytically determined temperature and the e fitted growth coefficient A from [35].Layers undergoing conversion from loading cycles, denoted as growth layers l ranging from 1 to 10, were geome lished in the FE model, as shown in Figure 3.

Finite Element Model
The sequentially coupled thermo-mechanical analysis comprised two steps: initially, a heat transfer simulation was performed to acquire the thermal input for the subsequent structural analysis and validate the derived analytical solution for the temperature field; subsequently, the structural analysis was executed utilizing the developed constitutive model from Section 2. Both heat transfer and structural analyses employed FE models with different settings, encompassing the input of material properties, element type, loading, and boundary conditions.
An axisymmetric FE model of the TBC system surrounding an air hole was constructed based on the geometry depicted in Figure 1b.The materials composing the substrate, BC, initial TGO, and TC were Inconel 718, NiCoCrAlY, α-Al 2 O 3 , and 8 wt.% yttria-stabilized ZrO2 (YSZ), respectively.The respective thicknesses were 1.2 mm, 150 µm, 1 µm, and 300 µm, with an additional air hole radius of 1.5 mm.The unit model encompassed one and a half wavelengths, incorporating features of peaks, valleys, and middles, as illustrated in Figure 3.The rough interface had a wavelength of 0.06 mm and an amplitude of 0.01 mm.The growth thickness of TGO in each cycle could be computed using Equation (11), utilizing the analytically determined temperature and the experimentally fitted growth coefficient A from [35].Layers undergoing conversion from BC to TGO in loading cycles, denoted as growth layers l ranging from 1 to 10, were geometrically established in the FE model, as shown in Figure 3.In the heat transfer simulation, the temperature histories depicted in Figure 2 were applied as boundary conditions to the top and bottom surfaces of the TBC system.In the structural simulation, the outer edge of the model could move radially but remained in a plane using the equation constraint [20,36], as illustrated in Figure 1b.The displacement in the z-direction of the bottom point at the outer edge, as depicted in Figure 1b, was con- In the heat transfer simulation, the temperature histories depicted in Figure 2 were applied as boundary conditions to the top and bottom surfaces of the TBC system.In the structural simulation, the outer edge of the model could move radially but remained in a plane using the equation constraint [20,36], as illustrated in Figure 1b.The displacement in the z-direction of the bottom point at the outer edge, as depicted in Figure 1b, was constrained to position the FE model.Additionally, the top and bottom surfaces, as well as the hole edge, were unconstrained in the structural analysis.
Perfect bonding was assumed between adjacent layers.The geometry of the model was partitioned appropriately prior to meshing.The FE model was meshed using eightnode axisymmetric elements with reduced integration schemes.DCAX8 and CAX8R element types were employed for heat transfer and structural analyses, respectively.Mesh refinements were locally applied near the hole edge, initial TGO, and TGO growth layers, as depicted in Figure 4. Sensitivity analyses of the meshing strategy were conducted to detect potential issues associated with element mixing and size.Thermo-mechanical responses, including temperature and stress, remained stable with a minimum element size of 0.03 µm in local regions.Additionally, there was no evidence of element mixing at the boundaries in the sensitivity analyses.The FE model comprised 20,504 elements and 61,963 nodes.This approach ensures adequate computational accuracy at reasonable costs.In the heat transfer simulation, the temperature histories depicted in Figure applied as boundary conditions to the top and bottom surfaces of the TBC system structural simulation, the outer edge of the model could move radially but remai plane using the equation constraint [20,36], as illustrated in Figure 1b.The displa in the z-direction of the bottom point at the outer edge, as depicted in Figure 1b, w strained to position the FE model.Additionally, the top and bottom surfaces, as the hole edge, were unconstrained in the structural analysis.
Perfect bonding was assumed between adjacent layers.The geometry of th was partitioned appropriately prior to meshing.The FE model was meshed usin node axisymmetric elements with reduced integration schemes.DCAX8 and CAX ment types were employed for heat transfer and structural analyses, respectivel refinements were locally applied near the hole edge, initial TGO, and TGO growth as depicted in Figure 4. Sensitivity analyses of the meshing strategy were condu detect potential issues associated with element mixing and size.Thermo-mecha sponses, including temperature and stress, remained stable with a minimum elem of 0.03 µm in local regions.Additionally, there was no evidence of element mixin boundaries in the sensitivity analyses.The FE model comprised 20,504 eleme 61,963 nodes.This approach ensures adequate computational accuracy at rea costs.

Material Properties
Table 1 lists the thermophysical properties, such as thermal conductivity k and specific heat C, of the four constituent layers [37][38][39][40], while Table 2 presents the thermal and mechanical properties, including the thermal expansion coefficient, Poisson's ratio, elastic modulus, and yield stress [40][41][42][43][44].The TGO, BC, and substrate layers utilized the embedded viscoplastic constitutive model with the parameters detailed in Table 3 [41], while the TC layer was assumed to be elastic [42].The unified viscoplasticity constitutive model for the superalloy Inconel 718, commonly employed in TBC systems, has been validated in prior research [20], demonstrating good agreement with experimental results [45].Elevated temperatures can induce sintering within the TC layer.The mechanical property evolution of the TC layer due to sintering was considered.During the high-temperature holding phase, micro-pores and micro-cracks underwent healing, while sintering contacts formed to bridge the gaps.The porosity of the TC layer decreased, leading to an increase in Young's modulus over time.Experimental observations indicated an increase in the Young's modulus of TC at a rate of 0.2% per hour under comparable conditions [46].The initial Young's modulus of the TC, denoted as E TC,0 , is provided in Table 2, while the evolving Young's modulus E TC was determined by the following equation: where a TC is the evolution rate of Young's modulus of the TC layer.

Temperature Fields: Analytical and Numerical Solutions
Heat transfer FE analysis was performed to validate the analytically predicted temperature field and to obtain the thermal load input for the subsequent thermo-mechanical coupling analysis.Prior to the formal heat transfer FE analysis, a preliminary heat transfer simulation accounting for material conversion occurring over multiple cycles was conducted, aiming to assess its impact on temperature values.This simulation utilized the user subroutine USDFLD, enabling the redefinition of field variables at material points, to facilitate the conversion of thermophysical properties from the BC to the TGO within the transforming layers.The results indicate that the material conversion resulting from TGO growth minimally affected the temperature distribution within the TGO layer.Therefore, neglecting TGO growth in the heat transfer analysis was justified.
As described in Section 4.1, each cycle comprised three phases: the heating phase, high-temperature holding phase, and cooling phase.Figure 5 illustrates the simulated temperature fields of the TBC system in various phases throughout a cycle.As shown in Figure 5a,e, the initial and final temperatures of the system were both 25 • C. Throughout the heating phase, the temperature gradually increased, accompanied by the emergence of a thermal gradient along the thickness direction.In the heating phase, the temperature difference between the top and bottom surfaces was 100 • C at 300 s, subsequently escalating to its peak of 200 • C during the high-temperature holding phase, as illustrated in Figure 5b,c.The thermal gradient exacerbated the thermal expansion mismatch, resulting in a more pronounced stress field near the air hole and interfaces.The proposed modelling framework utilizes the time-varying thermal outputs from the heat transfer FE analysis as inputs for the subsequent structural analysis, facilitating the incorporation of thermal gradient and time-varying temperature fields into the framework.Figure 6 presents the simulated and analytically calculated temperatures along the thickness direction of the TBC during the high-temperature holding phase, demonstrating excellent agreement between the two.The analytical formulae accurately predicted the non-uniform temperature field of the TBC during the high-temperature holding phase.Utilizing the analytically predicted temperature of the TGO layer to calculate the TGO growth rate in each cycle using Equation (11) was reliable.Moreover, the temperature of the TGO layer during the holding phase was 1044 • C. Figure 6 also shows that the temperature increased from 1000 • C at the bottom to 1200 • C at the top surface.The thermal gradient within the TC layer was much higher than in other layers due to its lower thermal conductivity.
a more pronounced stress field near the air hole and interfaces.The proposed modelling framework utilizes the time-varying thermal outputs from the heat transfer FE analysis as inputs for the subsequent structural analysis, facilitating the incorporation of thermal gradient and time-varying temperature fields into the framework.Figure 6 presents the simulated and analytically calculated temperatures along the thickness direction of the TBC during the high-temperature holding phase, demonstrating excellent agreement between the two.The analytical formulae accurately predicted the non-uniform temperature field of the TBC during the high-temperature holding phase.Utilizing the analytically predicted temperature of the TGO layer to calculate the TGO growth rate in each cycle using Equation ( 11) was reliable.Moreover, the temperature of the TGO layer during the holding phase was 1044 °C. Figure 6 also shows that the temperature increased from 1000 °C at the bottom to 1200 °C at the top surface.The thermal gradient within the TC layer was much higher than in other layers due to its lower thermal conductivity.

TGO Growth Thickness
Utilizing the temperature of the TGO layer during the holding phase, along with the experimentally determined growth parameters [35], enabled the calculation of the TGO growth rate and thickness for each cycle using Equations ( 11)- (13).TGO underwent growth during the high-temperature holding phase, with the TGO thickness over the holding time for ten cycles depicted in Figure 7.Over ten cycles, the TGO thickness increased from an initial value of 1 µm to 3.88 µm.The growth rate progressively decreased ulated and analytically calculated temperatures along the thickness direct during the high-temperature holding phase, demonstrating excellent agree the two.The analytical formulae accurately predicted the non-uniform tem of the TBC during the high-temperature holding phase.Utilizing the an dicted temperature of the TGO layer to calculate the TGO growth rate in ea Equation ( 11) was reliable.Moreover, the temperature of the TGO layer du ing phase was 1044 °C. Figure 6 also shows that the temperature increased at the bottom to 1200 °C at the top surface.The thermal gradient within th much higher than in other layers due to its lower thermal conductivity.

TGO Growth Thickness
Utilizing the temperature of the TGO layer during the holding phase, experimentally determined growth parameters [35], enabled the calculati growth rate and thickness for each cycle using Equations ( 11)- (13).TG growth during the high-temperature holding phase, with the TGO thick holding time for ten cycles depicted in Figure 7.Over ten cycles, the TGO creased from an initial value of 1 µm to 3.88 µm.The growth rate progressi

TGO Growth Thickness
Utilizing the temperature of the TGO layer during the holding phase, along with the experimentally determined growth parameters [35], enabled the calculation of the TGO growth rate and thickness for each cycle using Equations ( 11)- (13).TGO underwent growth during the high-temperature holding phase, with the TGO thickness over the holding time for ten cycles depicted in Figure 7.Over ten cycles, the TGO thickness increased from an initial value of 1 µm to 3.88 µm.The growth rate progressively decreased throughout the holding time.The decline in the TGO growth rate over time in TBCs can be ascribed to several factors, including aluminum consumption, the formation of protective oxide layers, and diffusion limitations, among others.throughout the holding time.The decline in the TGO growth rate over time in T be ascribed to several factors, including aluminum consumption, the formation of tive oxide layers, and diffusion limitations, among others.

Interfacial Stresses
Prior to examining the interfacial out-of-plane stresses, in-plane radial stresse the TGO were initially derived from the thermo-mechanical analysis outputs.Th mum compressive radial stress after the final cycle was 2.832 GPa, aligning with from previous studies in the literature [47,48].These compressive in-plane stress erbate the propensity for rumpling in the TGO layer, consequently fostering buck interface delamination.
Within TBCs, stress redistribution near the hole edge induces interfacial stress  perpendicular to the interface and shear stress  parallel to the in Tensile peeling stress and shear stress facilitate normal and shear interfacial sepa respectively.The combined influence of these interfacial stresses and interface im tions, such as undulations and voids, can lead to the initiation of interfacial crac the hole edge.Throughout continuous thermo-mechanical cycles, interfacial crack agate and merge with interfacial cracks or vertical cracks in the TC layer, ultim sulting in delamination and spallation failures [1,49].Figure 8 illustrates the evol  at the interface between growth layers 1 and 2 during the first cycle.During t ing phase,  increased to 625.5 MPa at 536 s due to rising temperature and therm match.Meanwhile, the yield strengths of the BC, TGO, and substrate materials de with rising temperature.After 536 s,  began to decrease due to the decrease strengths, reaching 597.6 MPa at 600 s.During the high-temperature holding pha time-dependent deformation occurred in and near the TGO layer due to the creep ity) effect, leading to a gradual decrease in interfacial stress, as depicted in Figu particular, the interfacial  gradually diminished to 495.5 MPa at the end of th temperature holding phase.The inhibitory effect on interfacial stresses caused b can delay the initiation of cracks.Figure 8 also illustrates a significant increase in cial  during the cooling phase, reaching 1625.6 MPa at the end of the first cy rise is primarily attributed to two factors.Firstly, completion of the conversion of layer 1 from the BC to the TGO by the end of the high-temperature holding phas first cycle positioned the investigated  at the interface between the TGO and BC where thermal mismatch was particularly severe.Secondly, stresses near the in entered a yielding state during the high-temperature holding phase, and as the t ture decreased, the yield strengths of the materials increased.The combined effec vere thermal mismatch and increasing yield strengths contributed to the signific

Interfacial Stresses
Prior to examining the interfacial out-of-plane stresses, in-plane radial stresses within the TGO were initially derived from the thermo-mechanical analysis outputs.maximum compressive radial stress after the final cycle was 2.832 GPa, aligning with findings from previous studies in the literature [47,48].These compressive in-plane stresses exacerbate the propensity for rumpling in the TGO layer, consequently fostering buckling and interface delamination.
Within TBCs, stress redistribution near the hole edge induces interfacial peeling stress σ zz perpendicular to the interface and shear stress τ zr parallel to the interface.Tensile peeling stress and shear stress facilitate normal and shear interfacial separations, respectively.The combined influence of these interfacial stresses and interface imperfections, such as undulations and voids, can lead to the initiation of interfacial cracks near the hole edge.Throughout continuous thermo-mechanical cycles, interfacial cracks propagate and merge with interfacial cracks or vertical cracks in the TC layer, ultimately resulting in delamination and spallation failures [1,49].Figure 8 illustrates the evolution of σ zz at the interface between growth layers 1 and 2 during the first cycle.During the heating phase, σ zz increased to 625.5 MPa at 536 s due to rising temperature and thermal mismatch.Meanwhile, the yield strengths of the BC, TGO, and substrate materials decreased with rising temperature.After 536 s, σ zz began to decrease due to the decrease in yield strengths, reaching 597.6 MPa at 600 s.During the high-temperature holding phase, slow time-dependent deformation occurred in and near the TGO layer due to the creep (viscosity) effect, leading to a gradual decrease in interfacial stress, as depicted in Figure 8.In particular, the interfacial σ zz gradually diminished to 495.5 MPa at the end of the high-temperature holding phase.The inhibitory effect on interfacial stresses caused by creep can delay the initiation of cracks.Figure 8 also illustrates a significant increase in interfacial σ zz during the cooling phase, reaching 1625.6 MPa at the end of the first cycle.This rise is primarily attributed to two factors.Firstly, completion of the conversion of growth layer 1 from the BC to the TGO by the end of the high-temperature holding phase of the first cycle positioned the investigated σ zz at the interface between the TGO and BC layers, where thermal mismatch was particularly severe.Secondly, stresses near the interfaces entered a yielding state during the high-temperature holding phase, and as the temperature decreased, the yield strengths of the materials increased.The combined effects of severe thermal mismatch and increasing yield strengths contributed to the significant rise in interfacial σ zz after the high-temperature holding phase and during the cooling phase.The evolutions of  near the hole edge at four interfaces are depicted in Figure 10. Figure 10b,d show a change in the direction of interfacial  from compression to tension after the layer above underwent oxidation, with the maximum interfacial  occurring simultaneously.This change in stress direction and magnitude is attributed to oxide growth, which introduced severe thermal mismatch at the evolving TGO/BC interface.Consequently, tensile peeling stress increased, enhancing the likelihood of interfacial crack initiation and propagation near the hole edge.Subsequently, the interfacial  The evolutions of σ zz near the hole edge at four interfaces are depicted in Figure 10. Figure 10b,d show a change in the direction of interfacial σ zz from compression to tension after the layer above underwent oxidation, with the maximum interfacial σ zz occurring simultaneously.This change in stress direction and magnitude is attributed to oxide growth, which introduced severe thermal mismatch at the evolving TGO/BC interface.Consequently, tensile peeling stress increased, enhancing the likelihood of interfacial crack initiation and propagation near the hole edge.Subsequently, the interfacial σ zz gradually decreased with cycle progression, as the layer below oxidized to the same material as the TGO above, moving away from the evolving TGO/BC interface and the associated thermal mismatch.This observation aligns with prior research [50].Furthermore, the maximum interfacial shear stress τ zr was observed in the middle region of the wavy interface, near the TC/TGO interface.This suggests that mode-II crack initiation due to interfacial τ zr is more likely to occur at the TC/TGO interface in regions between the peak and valley.gradually decreased with cycle progression, as the layer below oxidized to the same material as the TGO above, moving away from the evolving TGO/BC interface and the associated thermal mismatch.This observation aligns with prior research [50].Furthermore, the maximum interfacial shear stress  was observed in the middle region of the wavy interface, near the TC/TGO interface.This suggests that mode-II crack initiation due to interfacial  is more likely to occur at the TC/TGO interface in regions between the peak and valley.

Strain Energy Density
The strain energy represents the energy stored within the TBC due to deformations during thermo-mechanical loading cycles, while the strain energy density quantifies the amount of strain energy per unit volume.Irreversible plastic strain energies exacerbate the damage in TBCs caused by cyclic thermo-mechanical loading.The findings indicate that the total and plastic strain energy densities, denoted as  and  , respectively, were higher in the growth layers near the hole edge compared to other regions, with their maximum values increasing with each cycle.Particularly, location competition for the maximum  was found between the first growth layer and the evolving TGO/BC interface.The  near the first growth layer increased to 3.177 × 10 −5 J/m 3 after the first cycle, gradually rising to 3.279 × 10 −5 J/m 3 by the final cycle.Meanwhile, the maximum  near the evolving TGO/BC interface rose from 1.53 × 10 −5 J/m 3 after the first cycle to 3.42 × 10 −5 J/m 3 by the final cycle, surpassing the value near the first growth layer.This location competition for the maximum plastic strain energy implies a competition for potential failure initiation between the first oxide growth layer and the evolving TGO/BC interface in TBCs with air holes.Initially, energy-driven cracks were prone to initiate near the first oxide growth layer, while in later cycles, they tended to initiate near the evolving TGO/BC interface.

Strain Energy Density
The strain energy represents the energy stored within the TBC due to deformations during thermo-mechanical loading cycles, while the strain energy density quantifies the amount of strain energy per unit volume.Irreversible plastic strain energies exacerbate the damage in TBCs caused by cyclic thermo-mechanical loading.The findings indicate that the total and plastic strain energy densities, denoted as u t and u p , respectively, were higher in the growth layers near the hole edge compared to other regions, with their maximum values increasing with each cycle.Particularly, location competition for the maximum u p was found between the first growth layer and the evolving TGO/BC interface.The u p near the first growth layer increased to 3.177 × 10 −5 J/m 3 after the first cycle, gradually rising to 3.279 × 10 −5 J/m 3 by the final cycle.Meanwhile, the maximum u p near the evolving TGO/BC interface rose from 1.53 × 10 −5 J/m 3 after the first cycle to 3.42 × 10 −5 J/m 3 by the final cycle, surpassing the value near the first growth layer.This location competition for the maximum plastic strain energy implies a competition for potential failure initiation between the first oxide growth layer and the evolving TGO/BC interface in TBCs with air holes.Initially, energy-driven cracks were prone to initiate near the first oxide growth layer, while in later cycles, they tended to initiate near the evolving TGO/BC interface.
The evolutions of u t and u p in growth layers near the hole edge are plotted in Figure 11.They demonstrate that u p increases in steps before the layer above oxidizes to TGO; then, it experiences a more significant increment during the cooling phase of the cycle when the layer above oxidizes to TGO.Subsequently, it stabilizes.The findings concerning the strain energies demonstrate that, due to the combined influences of the hole edge, TGO growth, creep-plasticity interaction, interface undulation, cyclic thermo-mechanical loading, and temperature gradient, damage accumulates and intensifies near the hole edge and TGO interfaces over the cycles, rendering this region a potential failure zone susceptible to edge-delamination.In particular, prolonged TGO growth can markedly escalate damage accumulation and the likelihood of failure.Moreover, the results underscore the importance, in the design of TBCs for hot-section components, of not only ensuring that air holes meet the requirements for engine combustion and cooling but also evaluating the strength of the TBC system near the hole edge and TGO interfaces.The evolutions of  and  in growth layers near the hole edge are plotted in Figure 11.They demonstrate that  increases in steps before the layer above oxidizes to TGO; then, it experiences a more significant increment during the cooling phase of the cycle when the layer above oxidizes to TGO.Subsequently, it stabilizes.The findings concerning the strain energies demonstrate that, due to the combined influences of the hole edge, TGO growth, creep-plasticity interaction, interface undulation, cyclic thermo-mechanical loading, and temperature gradient, damage accumulates and intensifies near the hole edge and TGO interfaces over the cycles, rendering this region a potential failure zone susceptible to edge-delamination.In particular, prolonged TGO growth can markedly escalate damage accumulation and the likelihood of failure.Moreover, the results underscore the importance, in the design of TBCs for hot-section components, of not only ensuring that air holes meet the requirements for engine combustion and cooling but also evaluating the strength of the TBC system near the hole edge and TGO interfaces.

Summary and Concluding Remarks
A unified viscoplastic constitutive model incorporating TGO growth is developed and integrated into a finite element modelling framework for TBCs near air holes.A reliable analytical solution for the non-uniform temperature field is obtained.The framework comprehensively accounts for various complex factors such as dynamic TGO growth, creep-plasticity interaction, interface undulation, temperature gradient, and sintering, facilitating accurate modelling of the dynamic behaviors leading to hole edge delamination failure.Several key findings and conclusions are summarized as follows: (1) The derived analytical solution for the non-uniform temperature field of the TBC in the high-temperature holding phase agrees well with the simulation.Thus, utilizing

Summary and Concluding Remarks
A unified viscoplastic constitutive model incorporating TGO growth is developed and integrated into a finite element modelling framework for TBCs near air holes.A reliable analytical solution for the non-uniform temperature field is obtained.The framework comprehensively accounts for various complex factors such as dynamic TGO growth, creep-plasticity interaction, interface undulation, temperature gradient, and sintering, facilitating accurate modelling of the dynamic behaviors leading to hole edge delamination failure.Several key findings and conclusions are summarized as follows: (1) The derived analytical solution for the non-uniform temperature field of the TBC in the high-temperature holding phase agrees well with the simulation.Thus, utilizing the analytically predicted temperature of the TGO layer to compute the TGO growth thickness and growth strains is reliable.(2) The calculated in-plane stresses in the TGO layer, determined using the developed viscoplastic constitutive model integrating the dynamic growth of TGO, align well with experimental findings reported in the literature [47,48].(3) Integrating TGO growth into the viscoplastic constitutive model facilitates comprehensive analyses of stress evolution near the evolving TGO/BC interface, enhancing the comprehension of hole edge delamination mechanisms under the combined effects of multiple factors.TGO growth can shift the direction of peeling stress near an air hole from compression to tension, thus facilitating interface debonding.(4) Irreversible plastic strain energy concentrates in the TGO growth layers near the hole edge, intensifying damage induced by cyclic thermo-mechanical loading.Importantly, competition for the potential failure initiation site is revealed between the first TGO growth layer and the evolving TGO/BC interface.
Future experiments are anticipated to acquire and calibrate the material parameters in the developed viscoplastic constitutive model for the coating layers.Nonetheless, this study establishes a foundation for future modelling developments that may explore interfacial delamination, crack propagation, thermo-mechano-chemical coupling, and non-uniform oxidation during TGO growth, alongside the influences of grain size and orientation in TBCs with air holes.

Figure 1 .
Figure 1.Schematic diagrams of the TBC with a central air hole: (a) three-dimensional axisymmetric structure; (b) the axisymmetric unit model.

Figure 1 .
Figure 1.Schematic diagrams of the TBC with a central air hole: (a) three-dimensional axisymmetric structure; (b) the axisymmetric unit model.

Figure 2 .
Figure 2. Cyclic thermo-mechanical loading for the TBC system near an air hole.

Figure 2 .
Figure 2. Cyclic thermo-mechanical loading for the TBC system near an air hole.

Figure 3 .
Figure 3. Local geometry of the axisymmetric FE model of the TBC system around an air hole.

Figure 3 .
Figure 3. Local geometry of the axisymmetric FE model of the TBC system around an air hole.

Figure 3 .
Figure 3. Local geometry of the axisymmetric FE model of the TBC system around an air h

Figure 4 .
Figure 4. Local mesh refinements of the TBC FE model.

Figure 4 .
Figure 4. Local mesh refinements of the TBC FE model.

Figure 5 .Figure 6 .
Figure 5.The FE simulated temperature fields of the TBC in different phases during a cycle: (a) the initial state at t = 0 s; (b) the heating phase at t = 300 s; (c) the high-temperature holding phase at t = 7800 s; (d) the cooling phase at t = 8100 s; (e) the final state at t = 8400 s.

Figure 5 .
Figure 5.The FE simulated temperature fields of the TBC in different phases during a cycle: (a) the initial state at t = 0 s; (b) the heating phase at t = 300 s; (c) the high-temperature holding phase at t = 7800 s; (d) the cooling phase at t = 8100 s; (e) the final state at t = 8400 s.

Figure 5 .Figure 6 .
Figure 5.The FE simulated temperature fields of the TBC in different phases durin initial state at t = 0 s; (b) the heating phase at t = 300 s; (c) the high-temperature hold 7800 s; (d) the cooling phase at t = 8100 s; (e) the final state at t = 8400 s.

Figure 6 .
Figure 6.Simulated and analytically predicted temperatures along the thickness direction of the TBC in the high-temperature holding phase.

Figure 7 .
Figure 7. TGO thickness with respect to the holding time in ten cycles.

Figure 7 .
Figure 7. TGO thickness with respect to the holding time in ten cycles.

Figure 8 .
Figure 8.The evolution of interfacial peeling stress  at the growth layer 1/2 interface du first cycle.

Figure 9 Figure 9 .
Figure9displays the distributions of  near the hole edge and TGO after cycles.The results indicate that tensile  primarily resided in the area neighbo hole edge and the evolving TGO/BC interface.This suggests that mode-I crack an delamination, induced by tensile  , are more likely to occur at the TGO/BC i near the hole edge rather than at the TC/TGO interface.With cycle progression, hig nitude tensile peeling stresses gradually shifted downward along the hole edge the evolving TGO/BC interface.Additionally, the maximum residual  after o initially trended upward in early cycles, followed by a slight decrease in later cy tributed to a viscous effect induced by high-temperature environments.The spatia ment and magnitude escalation of interfacial  under cyclic thermo-mechanical exacerbate the likelihood of crack initiation at the TGO/BC interface near the hole TBC systems.

Figure 8 .
Figure 8.The evolution of interfacial peeling stress σ zz at the growth layer 1/2 interface during the first cycle.

Figure 9 Figure 8 .
Figure9displays the distributions of σ zz near the hole edge and TGO after various cycles.The results indicate that tensile σ zz primarily resided in the area neighboring the hole edge and the evolving TGO/BC interface.This suggests that mode-I crack and edge delamination, induced by tensile σ zz , are more likely to occur at the TGO/BC interface near the hole edge rather than at the TC/TGO interface.With cycle progression, highmagnitude tensile peeling stresses gradually shifted downward along the hole edge toward the evolving TGO/BC interface.Additionally, the maximum residual σ zz after one cycle initially trended upward in early cycles, followed by a slight decrease in later cycles, attributed to a viscous effect induced by high-temperature environments.The spatial movement and magnitude escalation of interfacial σ zz under cyclic thermo-mechanical loading exacerbate the likelihood of crack initiation at the TGO/BC interface near the hole edge in TBC systems.

Figure 9 Figure 9 .
Figure9displays the distributions of  near the hole edge and TGO after various cycles.The results indicate that tensile  primarily resided in the area neighboring the hole edge and the evolving TGO/BC interface.This suggests that mode-I crack and edge delamination, induced by tensile  , are more likely to occur at the TGO/BC interface near the hole edge rather than at the TC/TGO interface.With cycle progression, high-magnitude tensile peeling stresses gradually shifted downward along the hole edge toward the evolving TGO/BC interface.Additionally, the maximum residual  after one cycle initially trended upward in early cycles, followed by a slight decrease in later cycles, attributed to a viscous effect induced by high-temperature environments.The spatial movement and magnitude escalation of interfacial  under cyclic thermo-mechanical loading exacerbate the likelihood of crack initiation at the TGO/BC interface near the hole edge in TBC systems.

Figure 11 .
Figure 11.Evolutions of total and plastic strain energy densities near the hole edge in (a) initial TGO; (b) growth layer 4; (c) growth layer 8; (d) growth layer 10.

Figure 11 .
Figure 11.Evolutions of total and plastic strain energy densities near the hole edge in (a) initial TGO; (b) growth layer 4; (c) growth layer 8; (d) growth layer 10.

Table 1 .
Thermophysical properties of the constituent layers in the TBC system.

Table 2 .
Thermal and mechanical properties of the constituent layers in the TBC system.