Reliability Evaluation of EB-PVD Thermal Barrier Coatings in High-Speed Rotation and Gas Thermal Shock

Abstract

The uncertain service life of thermal barrier coatings (TBCs) imposes constraints on their secure application. In addressing this uncertainty, this study employs the Monte Carlo simulation method for reliability evaluation, quantifying the risk of TBC peeling. For reliability evaluation, the failure mode needs to be studied to determine failure criteria. The failure mode of high-speed rotating TBCs under gas thermal shock was studied by combining fluid dynamics simulations and experiments. Based on the main failure mode, the corresponding failure criterion was established using the energy release rate, and its limit state equation was derived. After considering the dispersion of parameters, the reliability of TBCs was quantitatively evaluated using failure probability and sensitivity analysis methods. The results show that the main mode is the fracture of the ceramic layer itself, exhibiting a distinctive top-down “step-like” thinning and peeling morphology. The centrifugal force emerges as the main driving force for this failure mode. The failure probability value on the top side of the blade is higher, signifying that coating failure is more likely at this location, aligning with the experimental findings. The key parameters influencing the reliability of TBCs are rotation speed, temperature, and the thermal expansion coefficient. This study offers a valuable strategy for the secure and reliable application of TBCs on aeroengine turbine blades.

1. Introduction

The aeroengine serves as the core of an aircraft, with the thrust-to-weight ratio being a pivotal parameter for the engine. In recent years, the engine’s gas inlet temperature has risen due to an increase in the thrust-to-weight ratio, reaching 1900 K in advanced turbofan engines [1,2,3]. Single-crystal materials and high-efficiency cooling air film technology are no longer sufficient for the demands of advanced hot-end engine components. Currently, the most effective method for enhancing the service temperature of the engine is the heat insulation protection technology provided by thermal barrier coatings (TBCs). TBCs consist of a heat-insulating ceramic layer, an oxidation-resistant and bond-enhancing intermediate layer, and a substrate made of nickel-based superalloy [4].

The TBCs are susceptible to peeling off in unforeseen circumstances, leading to turbine blade failure, owing to performance parameter mismatches across layers, intricate geometries, and challenging service conditions [5]. Thus, the study of the service life of TBCs is of paramount importance. Developing an accurate life prediction model is beneficial to the safe and reliable application of TBCs. This study revealed that the fracture toughness of TBCs follows the Weibull distribution [6]. Both the thermal expansion coefficient and Young’s modulus of TBCs follow a normal distribution [7]. Additionally, uncertainties exist in the service conditions and thickness of TBCs. The stochastic nature of these parameters will inevitably result in the dispersion of TBCs’ lifespan. Nevertheless, conventional life models for TBCs fail to account for uncertainties in material parameters, geometries, and loading conditions [8,9]. Predicted results can only offer the average lifespan and lack information on the degree of lifespan dispersion [10]. For instance, a singular life prediction model, such as maximum stress or certain specific damage parameters, is insufficient to assess the safety status of TBCs [11]. Consequently, a substantial gap persists between the actual and predicted lifespan.

Compared with life prediction, the service reliability evaluation considers the uncertainty of these factors from a probabilistic and statistical perspective [12]. It is first necessary to study the main failure modes that cause TBCs to peel off. The corresponding failure criteria are established based on the main failure modes. Subsequently, the statistical characteristics of each random variable influencing these failure modes are analyzed. Next, a suitable reliability calculation method is chosen to assess the reliability of TBCs during varied operational periods and under diverse working conditions. Meanwhile, the sensitivity of each parameter to the service life of the TBCs can also be analyzed to determine the key performance parameters of the service life of the TBCs. The reliability evaluation method considers the dispersion of material parameters, the randomness of geometric structures, and the uncertainty of loading conditions. This is used to quantify the peeling risk of the TBCs under variable service conditions.

At present, the main causes of TBC spalling failure include CMAS corrosion, high-speed rotation, erosion, gas thermal shock, and other service environments [13]. In these service environments, the interaction between high-speed rotation and gas thermal shock is the key factor causing TBC peeling [14]. It is also the most common service environment for TBCs on aeroengine turbine blades. However, existing research mainly focuses on the static service performance of the TBCs, without considering the effect of high-speed rotation. Therefore, it is crucial to study the service behavior and life prediction of the TBCs on turbine blades under high-temperature and high-speed rotation.

In this work, the reliability evaluation method of the Monte Carlo simulation is used to quantify the risk of spalling of the TBCs caused by high-speed rotation and gas thermal shock. For reliability evaluation, the failure mode first needs to be studied to determine failure criteria. The failure mode of high-speed rotating TBCs under gas thermal shock was studied by combining fluid dynamics simulations and experiments. Secondly, based on the main failure mode, the corresponding failure criterion was established using the energy release rate. And the limit state equation of TBC failure was obtained to determine the boundary between failure and safe regions. Then, by analyzing the random characteristics of the corresponding parameters in the failure criteria, the reliability of TBCs was quantitatively evaluated using failure probability and sensitivity analysis methods and compared with the experimental results. Finally, the key parameters affecting the reliability of high-speed rotating TBCs under gas thermal shock were obtained. The results are used to guide the application and optimization of the TBCs on aeroengine turbine blades.

2. Reliability Evaluation Method

2.1. Basic Theory

For reliability evaluation, it is first necessary to establish corresponding failure criteria based on the main failure modes of TBC peeling. The method draws on the reliability evaluation theory of engineering structures. This can be generalized into a functional formula, called the functional function [15]:

$$Z=g\left(X\right)=R\left(X\right)-S\left(X\right)$$

(1)

where X = X₁, X₂, ‧‧‧, X_n are random variables; R represents the ability of the TBCs to withstand load; S refers to the generalized load that causes the internal force and deformation in the TBCs. If Z > 0, the TBC is in a safe state; if Z < 0, the TBC is in a failure state; and when Z = 0, the TBC is in a critical failure state. Therefore, the equation Z = 0 is also called the limit state equation.

The probability of failure can be described by the following integral:

$$P_f={\displaystyle {\int }_{-\infty }^0{f}_Z\left(\mathcal{z}\right)d\mathcal{z}}={\displaystyle \int \cdots {\displaystyle {\int }_{\mathsf{\Omega }}{f}_X(x_1,x_2,\dots ,x_n)}}\mathrm{d}{x}_1\cdots \mathrm{d}{x}_n$$

(2)

where f_X(x₁, x₂,…, x_n) is the joint probability density function of X; $f_Z\left(\mathcal{z}\right)$ is the probability density function of Z; Ω represents the failure domain Z ≤ 0. If the random variables X_i are independent of each other, Formula (2) can be expressed as

$$P_f={\displaystyle {\int }_{\mathsf{\Omega }}{f}_X\left(x\right)}\mathrm{d}x={\displaystyle \int \cdots {\displaystyle {\int }_{\mathsf{\Omega }}{f}_{X_1}\left(x_1\right)}}{f}_{X_2}(x_2)\dots f_{X_n}(x_n)\mathrm{d}{x}_1\mathrm{d}{x}_2\dots \mathrm{d}{x}_n$$

(3)

where $f_{X_i}\left(x_i\right)$ is the probability density function of the random variable X_i. However, the solutions of these individual density functions and multiple integrals are often difficult to determine. Therefore, in recent years, many scholars have solved this problem through indirect methods [16,17,18]. The Monte Carlo method is used for calculations in this study [19].

2.2. Monte Carlo Simulation Method

In the Monte Carlo method, the failure probability is calculated by dividing the failure frequency by the total number of samples. Assuming that the total number of sampling times is N, X is first sampled according to $f_X\left(x\right)$, and then the sample value x is substituted into the functional function g(X) for calculation. If g(X) < 0, the TBC has failed, and the number of failures (N_f) will automatically increase by 1. According to Bernoulli’s theorem, the frequency (N_f/N) of random events g(X) < 0 converges to the probability of occurrence P_f. Hence, the predicted value is

$${\widehat{P}}_f=\frac{N_f}{N}$$

(4)

Based on Formula (2), the failure probability of the TBCs is

$$P_f={\displaystyle {\int }_{\mathsf{\Omega }}{f}_X\left(x\right)}\mathrm{d}x={\displaystyle {\int }_{-\infty }^{+\infty }I\left[g_X\right(x\left)\right]}{f}_X(x)\mathrm{d}x=E\left\{I\left[g_X\right(x\left)\right]\right\}$$

(5)

where $I\left[g_X\right(x\left)\right]$ is the indicator function of $g_X\left(x\right)$:

$$I[g_X(x)]=\left\{\begin{array}{cc}1& g_X(x)<0\\ 0& g_X\left(x\right)\ge 0\end{array}\right.$$

(6)

Therefore, according to Formula (5), the predicted value can be described as

$${\widehat{P}}_f=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left\{I\left[g_X\left({}^{i}x\right)\right]\right\}}$$

(7)

where $I[g_X{(}^{i}x\left)\right]$ is the subsample drawn from $I\left[g_X\right(x\left)\right]$. According to Formula (7), the mean of these samples is ${\widehat{P}}_f$. For large samples, regardless of the distribution of $I\left[g_X\right(x\left)\right]$, ${\widehat{P}}_f$ follows a normal distribution.

According to Monte Carlo principle, the specific steps for calculating the model are as follows:

(1) Assuming that n random variables X = [X₁, X₂,…, X_n]^T are independent of each other, the N-group data are generated according to the distribution law of the variable X;

(2) Each set of data represents a sample. If the samples falling into the failure domain Ω are N_f, then the failure probability is P_f = N_f/N.

2.3. Sensitivity Factor

To obtain the degree of influence on reliability when each parameter changes, we calculate the sensitivity factor of each parameter by fitting the failure probability and the quadratic function of each parameter. In this way, the key factors affecting the peeling failure of TBCs are extracted.

Assume the fitting function between the failure probability P_f of the obtained TBCs and any random parameter x is as follows:

$$P_f=g\left(x\right)$$

(8)

Then, the expression of the sensitivity factor can be as follows:

$$w_x=\frac{\partial g\left(x\right)}{\partial x}\left|{}_{x=x_r}\right.\times \frac{x_r}{P_{fr}}$$

(9)

where x represents a certain random parameter, x_r represents the value of the parameter under the reference condition, and P_fr is the failure probability under the reference condition.

3. Failure Modes of EB-PVD TBCs with High-Speed Rotation and Gas Thermal Shock

3.1. Fluid–Structure Interaction

3.1.1. Basic Governing Equations

The service environment of high-speed rotating EB-PVD (electron beam-physical vapor deposition) TBCs under gas thermal shock is extremely harsh. It is difficult to detect corresponding data, such as temperature fields and flow fields, in high-speed rotating real service environments. Computational fluid dynamics (CFD) can simulate the state of high-speed rotation under the thermal shock of gas in the engine. This study uses the commercial software ANSYS 2022R1 for simulation. The Navier–Stokes equation for compressible gases averaged by Reynolds is used to describe the fluid flow and heat transfer processes within an engine [20,21]. Gas air and cold air satisfy mass conservation, momentum conservation, and energy conservation in the flow process. The basic equation is as follows:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)=0$$

(10)

where ρ is the fluid density, t is the time, and u is the fluid velocity vector.

$$\frac{\partial \left(\rho u\right)}{\partial t}+\nabla \cdot (\rho uu)=-\nabla p+\nabla \cdot \tau +f$$

(11)

where p is the fluid pressure, f is the body force vector, τ is the viscous stress tensor, and its stress component τ_ij can be expressed as

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}{\delta }_{ij}$$

(12)

where μ is the hydrodynamic viscosity and δ_ij is the Kronecker delta function.

$$\frac{\partial \left(\rho h\right)}{\partial t}-\frac{\partial p}{\partial t}+\nabla \cdot (\rho uh)=\nabla \cdot (\lambda \nabla T)+\nabla \cdot (u\cdot \tau )+u\cdot \rho f+S_E$$

(13)

where h is the specific total enthalpy, λ is the thermal conductivity, T refers to the temperature, and S_E refers to the energy source. In addition, the gas in this study is assumed to be an ideal gas, and the expression of the ideal state equation is as follows:

$$p=\rho RT$$

(14)

where R is the ideal gas constant, R = 8.314 J/(mol∙K).

For the solid part, both the substrate and coating materials are assumed to be homogeneous elastic and isotropic materials. The governing equation is as follows:

$$\frac{\partial }{\partial x_j}({\sigma }_{ij})+F_i=0$$

(15)

$${\epsilon }_{ij}=\frac{1}{2}(\frac{\partial w_i}{\partial x_j}+\frac{\partial w_j}{\partial x_i})$$

(16)

$${\epsilon }_{ij}=\frac{1}{2G}({\sigma }_{ij}-\frac{\lambda }{2G+3\lambda }{\sigma }_{kk}{\delta }_{ij})+\alpha \Delta T{\delta }_{ij}$$

(17)

where σ_ij and ε_ij are the stress and strain components, respectively; F_ij is the body force component; w_i is the displacement component; and α is the thermal expansion coefficient. The shear modulus G and Lame constant λ can be expressed by Young’s modulus E and Poisson’s ratio ν, respectively:

$$G=\frac{E}{2(1+v)}$$

(18)

$$\lambda =\frac{vE}{(1+v\left)\right(1-2v)}$$

(19)

The fluid–structure interaction, similarly, adheres to fundamental conservation principles. At the fluid–solid coupling interface, the equality or conservation of variables such as fluid and solid stress τ, displacement d, heat flow q, temperature T, and others must be maintained. This requirement is expressed through the following set of four equations:

$$\left\{\begin{array}{l}{\tau }_f\cdot n_f={\tau }_s\cdot n_s\\ d_f=d_s\\ q_f=q_s\\ T_f=T_s\end{array}\right.$$

(20)

where the subscript f represents fluid, and the subscript s represents solid.

This paper adopts the two-equation turbulence model proposed by Menter [22], namely the SST k-ω model. This model combines the characteristics of the k-ω model used in the near-wall region and the k-ε model used in the far-wall region. The model uses a hybrid model to make the transition between the two models smooth. Therefore, the SST model has been widely used in engineering. The transport equation of the SST k-ω turbulence model is

$$\frac{\partial }{\partial t}(\rho k)+\frac{\partial }{\partial x_i}(\rho k{u}_i)=\frac{\partial }{\partial x_j}\left[(\mu +{\sigma }_k{\mu }_t)\frac{\partial k}{\partial x_j}\right]+P_k-{\beta }^{\ast }\rho \omega k+S_{k3}$$

(21)

$$\begin{array}{c}\frac{\partial }{\partial t}(\rho \omega )+\frac{\partial }{\partial x_i}(\rho \omega u_i)=\frac{\partial }{\partial x_j}\left[(\mu +{\sigma }_{\omega }{\mu }_t)\frac{\partial \omega }{\partial x_j}\right]+\frac{\gamma }{{\nu }_t}{P}_k-\beta \rho {\omega }^2+S_{\omega 3}\\ +2(1-F_1)\frac{\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\end{array}$$

(22)

where k is the turbulent kinetic energy, μ_t is the turbulent viscosity coefficient, P_k is the generation item of turbulent kinetic energy, ω is the specific dissipation rate, and β*, σ_k, γ, β, σ_ω, and σ_ω2 are the model constants. On the right side of the transport equation are the diffusion term, the generation term of turbulent kinetic energy, the dissipation term, and the source term. The fifth term on the right side of Formula (22) is the staggered diffusion term.

3.1.2. Geometric Modeling and Parameter Setting

To realize the service environment of high-speed rotating TBCs under gas thermal shock, a dynamic simulation test platform for TBCs was independently developed, as shown in Figure 1 [23]. The platform was mainly composed of a high-speed rotation module (Figure 1a), a gas thermal shock module (Figure 1b), and other auxiliary modules. The sample coated with the TBCs was fixed on the turbine disk at a circumference of 120 degrees, and the turbine disk was connected to the high-speed rotating shaft to realize the state of high-speed rotation of the TBCs, as shown in Figure 1a. Six supersonic flame spray guns were symmetrically arranged in the annular direction in this study, as shown in Figure 1b. Atomized aviation kerosene was used as the combustion medium, and oxygen and compressed air were used as combustion aids to fully burn to generate high-temperature airflow, and finally realize the simulation of high-temperature and high-speed gas.

To better simulate and verify the macro-failure location and micro-failure mode of the TBCs, the geometric model designed in this paper highly restores the main structure in Figure 1, which is closer to the test results. The geometric model of the turbine with a high-speed rotating TBCs is shown in Figure 2. The geometric modeling consists of the fluid domain and the solid domain. The fluid domain includes the inlet domain, the rotation domain, and the outlet domain. Six high-temperature gas inlets and one low-temperature cold air inlet are designed in the inlet domain. An outlet is set up in the outlet domain. The solid domain consists of an integrated simplified turbine model and simulated blades with TBCs arranged at 120° intervals around them. It should be noted that to study the influence of turbine blade curvature radius, this work uses cylindrical specimens instead of turbine blades.

Six gas inlets are designed on the circumference of the inlet domain. The inlet of the cold air is arranged at the center of the inlet domain. The outlet is set directly above the outlet area to form a pressure row. The inlet domain adopts the boundary conditions of total temperature and total pressure, and the outlet domain adopts the static pressure form. The boundary conditions in the fluid domain are shown in

Table 1. The specific boundary parameters for the fluid domain.

Pg,inlet/MPa	Tg,inlet/K	Pc,inlet/MPa	Tc,inlet/K	Poutlet/MPa	Rp/rpm
1.2	1273	0.8	373	0.5	0/3000/4000/5000

. P_g,inlet and T_g,inlet are total inlet pressure and total inlet temperature of high-temperature gas, respectively. P_c,inlet and T_c,inlet are the total inlet pressure and the total inlet temperature of the low-temperature cold air, respectively. P_outlet is the outlet static pressure. The rotation domain rotates around the z-axis at a velocity of R_p. In addition, the displacement in the x and y directions is free, but there is no displacement in the z direction.

The composite TBCs are simplified as a single-layer ceramic layer with a thickness of 200 μm and a radius of curvature of 5 mm. The internal defects and microstructure of the EB-PVD TBCs are ignored. The roughness of the TBC surface is a fixed value. Meanwhile, both the substrate and TBCs are assumed to be homogeneous elastic and isotropic materials. These material parameters are not related to temperature. The specific material parameters of the substrate and coating are shown in

Table 2. The specific material parameters of the substrate and TBCs.

Material Parameters	Substrate	YSZ-TBCs
Density (kg/m3)	8700	5500
Young’s modulus (GPa)	200	100
Poisson’s ratio	0.35	0.2
Thermal Conductivity (W/(m∙K))	88	1.65
Thermal expansion coefficient (1/°C)	18 × 10−6	12 × 10−6
Specific heat capacity (J/(kg∙K))	440	418

[23]. The boundary conditions for the surface of the rotating blade and the wall of the fluid domain are set as no-slip coupled heat transfer.

3.2. Result Analysis and Discussion

3.2.1. Analysis of Flow Field Characteristics

To better study the effect of high-speed rotation on the flow field characteristics around the TBCs, the flow field characteristics of the TBCs in the static state were also analyzed under the same conditions. In Figure 3a, the high-temperature and high-speed airflow flows in from the inlet, and after being hindered by the front of the simulated blade, it flows out through the two sides of the blade without changing the direction of the airflow significantly. From Figure 3b–d, it can be intuitively observed that the interaction between high-speed rotation and high-temperature gas leads to a significant change in the angle of airflow reaching the coating surface. At the same time, the deflection angle (the angle between the coating side and the airflow) is directly related to the speed of the rotation speed. As the rotational speed increases, the deflection angle increases to a certain extent. Also, the direction of deflection is always opposite to the direction of rotation (counterclockwise). Compared with the static state, the high-speed rotating state of the columnar coating directly changes the overall direction of the airflow in the flow field, and the deflection of the airflow further accelerates the fluid. Its formation principle is similar to the wake effect, which strongly indicates that the high-speed rotation state has a more significant impact on the flow field characteristics of the TBCs.

Figure 4 shows the Mach number nephogram of the TBC section at different rotational speeds. The front in the picture shows the location of the thermal shock. Figure 4a corresponds to the static condition, and the Mach number reaches its maximum at the left and right sides of the TBCs of the simulated blade. And the airflow velocity near the TBCs is as high as 1.03 Mach. Since the simulated blade is columnar rather than a real curved surface structure, it will affect the drainage effect, resulting in a low-velocity zone at the position directly behind the simulated blade. However, the low-velocity region of the real blade usually appears at the leading edge of the blade [24]. Figure 4b–d show the simulation results in the rotating state. It can be found that after the rotation, the high-speed airflow is deflected to the right, resulting in a shift in the positions of the high-speed zone and the low-speed zone. And the higher the speed is, the higher the Mach number on the right side of the simulated blade. It should be noted that, to ensure the comparability of the corresponding results, all the cross-sections in this paper are selected at the same height position, that is, the position above the middle of the blade.

The pressure nephogram of the TBC section under static and rotating conditions is shown in Figure 5. The position with the highest pressure under the static state in Figure 5a appears directly in front of the blade section, while the lower position appears on both sides of the blade. The high-pressure zone forms primarily due to the obstruction of high-speed airflow at the front end of the blade, creating a stagnation point where the pressure is generally the highest and the airflow velocity is the smallest. Conversely, the low-pressure area exhibits the opposite characteristics. The positions of both the high-pressure and low-pressure areas shift under the rotational state, as depicted in Figure 5b–d. With an increase in rotational speed, the high-pressure zone gradually shifts to the left, reaching a pressure as high as 1.2 MPa. However, as the speed continues to rise, the low-pressure area on the left diminishes, while the low-pressure area on the right expands. This phenomenon is associated with the flow velocity after the deflection of the flow field, aligning with the findings in Figure 4. Therefore, the change in rotational speed influences the size of the deflection angle of the flow field, affecting the position distribution of the velocity field and the pressure field.

3.2.2. Temperature Field of EB-PVD TBCs

The flow characteristics of the fluid are closely related to the distribution of the surface temperature of the TBCs. Therefore, based on the above velocity and pressure analysis, it is necessary to deeply analyze the temperature field of the TBCs. Figure 6 shows the temperature field cloud map of the surface of the TBCs in both stationary and rotating states. By selecting planes in four directions, the temperature distribution on the surface of the TBCs is intuitively and comprehensively displayed, wherein the sequence numbers marked in the figure are the front, left, back, and right in descending order. In Figure 6a, the red high-temperature areas on the surface of the TBCs in the static state are mainly distributed in positions ① directly in front of and ③ directly behind the thermal shock of the gas. The temperature at positions ② and ④ decreases sequentially from top to bottom, and the cooling effect is obvious. At the same time, the same result can be obtained from the flow field temperature distribution of the TBC section in Figure 7a. Subsequent experiments can prove that the position directly in front of the gas thermal shock is also the spalling area of the TBCs in the static state.

The temperature distribution on the TBC surface under a high-speed centrifugal load is quite different from that in Figure 6a, as shown in Figure 6b–d. The high-temperature region’s distribution area is significantly reduced, signifying that high-speed rotation enhances the overall cooling effectiveness of the blade coating surface. And the high rotation speed is more conducive to the overall cooling of the TBCs than the low rotation speed. Secondly, the high-temperature region appears at position ①, and its position shifts to the left, which is related to the trend of the flow field temperature distribution, as shown in Figure 7b–d.

3.2.3. Stress Field of EB-PVD TBCs

The high-speed rotating TBCs are subjected to centrifugal load, as shown in Figure 8. The top position of the TBCs receives the largest centrifugal load, which is directly related to the radius of rotation. In addition, as the rotational speed increases, the centrifugal load on the TBCs increases accordingly. When the rotational speed is 5000 rpm, the centrifugal load is as high as 3.32 × 10⁸ N/m³. According to the test results in Figure 9, the failure occurred at the top of the TBCs, indicating that the failure is closely related to the centrifugal load. Therefore, the experimental results are in good agreement with the simulation results.

In Figure 9, the peeling position of the high-speed rotating TBCs is located on the left side of the top, that is, the side where the counterclockwise rotation is first tangent to the high-temperature airflow inlet. Therefore, in addition to the action of the centrifugal load, it is also subjected to high-speed erosion of the outer wall surface of the coating by high-temperature gas. The load is expressed as the shear stress when the fluid is in contact with the solid wall, as shown in Figure 10. Marks ① and ② refer to the front and left sides of the TBCs, respectively. Figure 10b–d correspond to different rotational states, and it is evident that their high-stress regions are distributed on the top left side of the TBCs. Consequently, under the influence of high-speed centrifugal load and airflow shear, spalling failure of the TBCs occurs at the positions indicated in Figure 9b–d. In Figure 9a, the failure position of the TBCs in a static state is directly in front. An analysis of Figure 10a reveals that the region with the highest wall shear stress is apparent on both sides, suggesting that the peeling of the TBCs in the static state is not necessarily linked to shear stress. Thus, the rotational speed significantly influences TBC failure. Due to the modest rotational speed gradient, the failure of TBCs at different rotational speeds is relatively similar. Additionally, the maximum value of wall shear stress is approximately 2500 Pa, which does not reach the MPa level. This observation indicates that the exfoliation of high-speed rotating TBCs is primarily driven by centrifugal load, with the wall shear stress formed by high-speed airflow serving as the secondary driving force.

3.3. Failure Modes of EB-PVD TBCs

To delve into the specific failure of high-speed rotating EB-PVD TBCs under gas thermal shock, the SEM technique was employed to characterize the microscopic failure morphology of the TBCs, as depicted in Figure 11. The analysis reveals that the bonding interface between the ceramic layer and the bonding layer of the TBCs is relatively intact, with no observable internal oxidation phenomenon. The TGO exhibits slow and dense growth, reaching a maximum thickness of only 1.35 μm. Consequently, the TGO is not the primary cause of TBC spalling, as illustrated in Figure 11a.

Obviously, the columnar crystals are severely damaged or even fractured in the ceramic layer. The fracture and exfoliation of individual columnar crystals were captured, as shown in Figure 11c. The ceramic layer presents a unique top-to-bottom “stepped” thinning and exfoliation morphology, and the fracture positions at different heights are extremely flat. Further magnification of the surface of the exfoliation area reveals that the exfoliation occurs in the columnar crystal region. So, it can be confirmed that the fracture failure occurs inside the ceramic layer, as shown in Figure 11b. Therefore, compared with the static TBCs, the microscopic failure location of the high-speed rotating TBCs occurs inside the ceramic layer.

Building upon the preceding analysis and simulation results, an in-depth examination of the spalling failure of TBCs was conducted, laying the groundwork for the establishment of subsequent failure criteria. Figure 12 illustrates the schematic diagram of the spalling failure of EB-PVD TBCs. The spalling failure of TBCs under high-speed rotation is primarily associated with the centrifugal load. Additionally, high-speed airflow plays a significant role in the erosion of TBCs.

As shown in Figure 12a, the centrifugal load F_c during the high-speed rotation produces a shear stress inside the columnar crystals of the ceramic layer, which causes shear failure of the columnar crystals of the ceramic layer. Meanwhile, the scouring of TBCs by the high-speed airflow also manifests as shear failure. Under the action of these two kinds of shear stress, many shear cracks and transverse cracks are initiated in the columnar crystal region. When the shear stress reaches the fracture toughness of the columnar crystals, the columnar crystals of the ceramic layer fracture. It should be noted that, when establishing the failure criterion of the high-speed rotating TBCs, the centrifugal force is used as the main factor to derive the shear stress borne by the columnar crystals of the ceramic layer. The shear stress caused by the high-speed airflow is ignored.

As shown in Figure 12b, it can be found that the failure mode of the TBCs is mainly the fracture of the ceramic layer itself. The columnar crystals exhibit a “stepped” thinning and exfoliation mode, until the interface of the equiaxed crystal is completely peeled off. Under the joint action of a centrifugal load and a high-speed airflow scouring external load, a small number of microcracks are initiated at the heads of columnar crystals. As the number of cycles increases, the microcracks will further expand and merge. This results in more macroscopic shear cracks and macroscopic transverse cracks in the middle and lower parts of the columnar crystal region. Finally, the columnar grains are distorted and deformed by thinning and exfoliation until they are thinned to the interface between the columnar grain area and the equiaxed grain area. The failure mode presents a similar “step” shape on the interface topography, and an approximate “peak and valley” on the surface topography.

4. Reliability Evaluation of EB-PVD TBCs

4.1. Failure Criterion

Through the preceding analysis, it is determined that the spalling failure of TBCs is primarily induced by centrifugal force, resulting in the fracture of the ceramic layer during service. At this juncture, the cracks formed in the ceramic layer predominantly manifest as transverse cracks within a single columnar crystal unit (with the crack direction perpendicular to the axis of the columnar unit), as illustrated in Figure 11. Consequently, when the energy release rate (G) equals the fracture toughness (Γ_TBC) of the ceramic layer, expressed as G = Γ_TBC, a fracture occurs in a single columnar crystal unit [25]. Subsequently, the limit state equation can be expressed as

$$Z={\Gamma }_{\mathrm{TBC}}-G$$

(23)

If Z > 0, the TBC is in the safe and reliable state; if Z < 0, the TBC is in the failure domain of the failure state; when Z = 0, the TBC is in a critical failure state.

According to Hutchinson and Suo [26], for simple strain conditions, the energy release rate G when the coating cracks from the free boundary satisfies

$$G=H\frac{(1-{\upsilon }_{\mathrm{TBC}}^2){\sigma }_{\mathrm{TBC}}^{2}h}{E_{\mathrm{TBC}}}$$

(24)

where h is the thickness of the TBCs, and E_TBC and υ_TBC are the Young’s modulus and Poisson’s ratio of the TBCs, respectively. H is a dimensionless driving force. The value for H is taken to be 0.343 from Hutchinson and Suo [26], which means that cracks from the edge lead to peeling. The value is a rough approximation based on experimental results. In this study, σ_TBC is the centrifugal stress experienced by the high-speed rotating TBCs under gas thermal shock.

Therefore, the limit state equation can be written as

$$Z={\Gamma }_{\mathrm{TBC}}-G={\Gamma }_{\mathrm{TBC}}-H\frac{(1-{\upsilon }_{\mathrm{TBC}}^2){\sigma }_{\mathrm{TBC}}^{2}h}{E_{\mathrm{TBC}}}$$

(25)

Under gas thermal shock, the centrifugal force F and angular velocity ω of the high-speed rotating TBCs can be expressed as

$$F=m{\omega }^{2}R$$

(26)

$$\omega =2\pi n$$

(27)

where m is the mass of the turbine blade, n is the rotational speed, and R is the radius of rotation.

According to Ghosh’s research [27,28], the critical angular velocity of material failure of high-speed rotation under gas thermal shock is related to the temperature distribution, the elastic constant of materials, the thermal expansion coefficient, and the geometric size, etc., which can be expressed as

$$\omega ={\left[\frac{2T{E}_{\mathrm{TBC}}{\alpha }_{\mathrm{TBC}}}{(3+{\upsilon }_{\mathrm{TBC}}){\rho }_{\mathrm{TBC}}{R}^2}\right]}^{1/2}$$

(28)

where ω is the angular velocity, T is the temperature field distribution on the surface of the TBCs, α_TBC is the thermal expansion coefficient of the TBCs, ρ_TBC is the density of the TBCs, and υ_TBC is the Poisson’s ratio of the TBCs.

From Formulas (26)–(28), it can be concluded that

$$F={\left[\frac{8{\pi }^2{m}^2{n}^{2}T{E}_{\mathrm{TBC}}{\alpha }_{\mathrm{TBC}}}{(3+{\upsilon }_{\mathrm{TBC}}){\rho }_{\mathrm{TBC}}}\right]}^{1/2}$$

(29)

Then, the tensile stress can be estimated according to the centrifugal force, namely:

$${\sigma }_{\mathrm{TBC}}={\left[\frac{8{\pi }^2{m}^2{n}^{2}T{E}_{\mathrm{TBC}}{\alpha }_{\mathrm{TBC}}}{(3+{\upsilon }_{\mathrm{TBC}}){\rho }_{\mathrm{TBC}}}\right]}^{1/2}/S$$

(30)

where S is the cross-sectional area of the columnar sample.

Since it is a cylindrical sample, the area S of the cross-section is

$$S=\pi r^2$$

(31)

where r is the radius of curvature of the cylindrical sample.

From (30) and (31), we know the following:

$${\sigma }_{\mathrm{TBC}}={\left[\frac{8m^2{n}^{2}T{E}_{\mathrm{TBC}}{\alpha }_{\mathrm{TBC}}}{(3+{\upsilon }_{\mathrm{TBC}}){\rho }_{\mathrm{TBC}}{r}^4}\right]}^{1/2}$$

(32)

Substituting Equation (32) into Equation (25), the failure criterion of the high-speed rotating TBCs under gas thermal shock can be obtained as follows:

$$Z={\Gamma }_{\mathrm{TBC}}-G={\Gamma }_{\mathrm{TBC}}-H\frac{8m^2{n}^{2}hT{\alpha }_{\mathrm{TBC}}(1-{\upsilon }_{\mathrm{TBC}}^2)}{(3+{\upsilon }_{\mathrm{TBC}}){\rho }_{\mathrm{TBC}}{r}^4}$$

(33)

4.2. Stochastic Characteristics of Parameters

As mentioned earlier, the fracture of the ceramic layer is mainly due to the stress generated by the centrifugal load and the scour of high-temperature gas. These stresses are influenced by material properties, geometry factors, and loading conditions. Therefore, it can be seen from Equation (33) that there are three variables in the probability analysis of high-speed rotating TBC failure under thermal shock: material parameters (such as fracture toughness, coefficient of thermal expansion, Poisson’s ratio, and density), external loads (such as temperature, rotational speed), and geometric parameters (such as coating thickness, radius of curvature). To determine the boundary between TBC failure and reliability, multiple unrelated variables of each type are included in a limit state function.

The measured material properties, encompassing density, fracture toughness, thermal expansion coefficient, and Poisson’s ratio, exhibit considerable dispersion. It was found that the fracture toughness can be described by a Weibull distribution, while the thermal expansion coefficient and Poisson’s ratio are normally distributed [29,30]. There are also changes and uncertainties in the geometric properties such as the thickness of the TBCs and the radius of curvature of the sample [31]. In addition, the temperature distribution and rotational speed on the TBC surface also show randomness [32]. Therefore, these dispersion characteristics were considered in probability analysis to quantify the peeling risk of TBCs, which is referred to as reliability. The probability characteristics of the random parameters in the model are shown in

Table 3. The probability characteristics of the random parameters in the model.

Stochastic Variables	Mean Value	Standard Deviation	Distribution
ΓTBC (J/m2) [33,34,35,36]	45	9	Weibull
ρTBC (kg/m3)	5500	110	Normal
αTBC (1/°C) [37,38]	7.3 × 10−6	0.37 × 10−6	Normal
T (°C) [19,32]	Tc	Tc/10	Normal
n (rpm)	5000	108.6	Normal
h (μm) [39]	268	11.46	Normal
r (mm)	5	0.2	Normal

. It should be noted that the Poisson’s ratio of the TBCs is taken as a constant value of 0.23. The mass of the turbine blade is taken as a constant value of 134 g.

4.3. Analysis of Reliability and Sensitivity

4.3.1. Reliability Analysis

From the failure criterion in Section 4.1, it can be found that function (33) has strong nonlinearity, so the Monte Carlo method in Section 2.2 is used for the calculation. In this paper, the calculation steps of the Monte Carlo method are realized through MATLAB programming. Therefore, the reliability of the TBCs is analyzed, and the failure probability under different parameter distribution conditions and the relationship between the failure probability and the variation in each parameter are predicted. Specifically, the failure probability of fracture toughness Γ_TBC, coating density ρ_TBC, coating thermal expansion coefficient α_TBC, coating surface temperature T, rotational speed n, coating thickness h, and radius of curvature r are predicted, as shown in Figure 13.

Figure 13a,b,g illustrate that the failure probability of TBCs under high-speed rotation decreases with increasing fracture toughness, coating density, and curvature radius. In the service environment of gas thermal shock and high-speed rotation, a higher fracture toughness, coating density, and curvature radius lead to an increased reliability and reduced susceptibility to spalling failure. This aligns with findings in the literature [8,40]. Conversely, Figure 13c–f reveal an increase in the failure probability of TBCs under high-speed rotation with the rise in thermal expansion coefficient, surface temperature, rotation speed, and coating thickness. In other words, higher values of these parameters make TBCs more prone to peeling failure, posing an increased risk during service. This observation is consistent with experimental results in the literature [41,42,43]. Additionally, the failure probability under specific parameter values can be calculated to assess the reliability and safety of the TBCs. At a speed of 5000 rpm and a temperature of 1000 °C, combined with other random variable values from Table 1, the predicted failure probability of the TBCs is 82.54%. (The MATLAB program for calculating the failure probability using the Monte Carlo method is available in Appendix A.) This suggests a potential risk in this service condition.

To further determine the failure location of the TBCs of the turbine blade, and at the same time consider the non-uniform temperature distribution of the TBC surface, it is necessary to obtain the temperature distribution of the TBC surface. However, the actual temperature distribution on the surface of the TBCs is difficult to obtain through experiments. In this study, computational fluid dynamics is employed to acquire the temperature distribution on the blade surface, as illustrated in Figure 6. The temperature data (T_c) for each node on the TBCs of the blade surface are then obtained. Subsequently, the data from Table 3 and the temperature data are inputted into the Monte Carlo model to calculate the failure probability for each node. The above steps are automatically repeated to calculate the failure probabilities of all nodes. Consequently, the distribution of failure probabilities for TBCs on the blade surface is presented in Figure 14. Marks ① and ② in the figure indicate the front and left sides of the TBCs, respectively.

From the highlighted red area in Figure 14a, it is evident that the uppermost left end of the columnar simulation specimen for TBCs has the highest failure probability value, reaching 79.4%. This signifies a pronounced susceptibility to coating peeling failure in these specific areas. In essence, in the service conditions of gas thermal shock and high-speed rotation, these regions are deemed unreliable, posing a heightened risk during service. Conversely, the failure probability values in other areas are smaller, indicating an elevated reliability and safety of TBCs in those specific regions under this service condition. To further validate the accuracy of the simulation results, experiments were conducted using the self-developed simulation device designed for the TBC service environment of high-speed rotating turbine blades, as depicted in Figure 1. The experimental results in Figure 14b reveal that under the same service environment, the peeling failure of TBCs indeed occurs at the upper left end of the columnar simulation specimen, aligning closely with the simulation results. This demonstrates the appropriateness of the established failure criterion for the reliability calculation under this specific failure mode. Consequently, it presents a valuable strategy for ensuring the safe and reliable service of rotating blades.

4.3.2. Sensitivity Analysis

The quadratic function of the failure probability and the mean value of each random parameter is fitted by the weighted least-squares method, as shown in Figure 13. Then, the average sensitivity factor of each random parameter is calculated according to Formula (9). Hence, the sensitivity analysis for the TBCs is shown in Figure 15. The sensitivity of a random variable is described by the mean. If the absolute value of the sensitivity is larger, it means that the random variable is more sensitive to the failure of TBCs. That is to say, the average value of the random variables Γ_TBC, ρ_TBC, and r increases, and the reliability increases. The opposite takes place for the random variables α_TBC, T, n, and h. This shows that the reliability of the TBCs is proportional to the radius of curvature and inversely proportional to the rotational speed and temperature, which is in good agreement with the experimental results [23]. In addition, n, α_TBC, and T have the highest sensitivity, followed by ρ_TBC, r, and h. The sensitivity of Γ_TBC is relatively small. Therefore, the reliability of TBCs is most sensitive to n, α_TBC, and T.

The insulation technology of the TBCs is internationally recognized as the most practical method for significantly increasing the service temperature of aeroengine turbine blades. However, the premature failure of TBCs represents a bottleneck that restricts its application and development, posing a key issue that urgently requires resolution. In instances where the peeling of TBCs is notably prominent, it becomes imperative to accurately predict the timing, location, and form of the peeling that will occur in the service environment. Subsequently, identifying the key influencing factors and establishing safe service conditions from the mechanism becomes crucial for guiding the secure application of TBCs. Therefore, predicting the service life of TBCs emerges as an urgent necessity in this research field.

This study employs probabilistic and statistical methodologies to assess the reliability of TBCs on rotating blades. The method considers the dispersive nature of TBC spalling, enabling a more accurate reliability assessment. However, the current service life model lacks consideration for the dispersion in the microstructure, performance, and service environment of TBCs, leading to predictions that deviate significantly from actual engineering conditions. Existing models, such as the damage variable life model by Busso et al. [44], the life model based on crack acoustic emission signals proposed by Renusch et al. [45], and the life model based on interface crack growth proposed by He et al. [10], lack practical engineering applications. Furthermore, the interaction between high-speed rotation and gas thermal shock is a critical factor contributing to the peeling of TBCs. Most existing research primarily focuses on the static service reliability of TBCs, neglecting the consideration for high-speed rotation. In contrast, this study incorporates high-speed rotation, providing a comprehensive evaluation of reliability for moving blades.

5. Conclusions

In this paper, the reliability evaluation method is used to quantitatively evaluate the reliability of high-speed rotating TBCs. Meanwhile, the results are compared with the experimental results. The main conclusions are as follows:

(1)

The main driving force causing the failure of the TBCs is the centrifugal force, followed by the wall shear stress caused by the interaction of high-speed rotation and high-speed airflow. Under the influence of this driving force, the failure mode of EB-PVD TBCs is mainly the self-fracture of the ceramic layer. It shows a unique top-down “step-like” thinning and peeling morphology.

(2)

There are three random variables affecting the reliability of high-speed rotating TBCs, namely material parameters (Γ_TBC, α_TBC, and ρ_TBC), external loads (T, n), and geometric parameters (h, r). The reliability of TBCs increases with the increase in fracture toughness Γ_TBC, coating density ρ_TBC, and radius of curvature r. It decreases with the increase in thermal expansion coefficient α_TBC, temperature T, rotation speed n, and coating thickness h.

(3)

Through reliability calculation, the failure probability value at the uppermost left end of the TBCs of the columnar sample is larger, indicating that the coating fails at this location. This is consistent with the experimental results. According to the sensitivity analysis, the most important parameters affecting reliability are n, α_TBC, and T, followed by ρ_TBC, r, and h.