Shakedown Analysis and Experimental Study of Thermal Barrier Coatings

Abstract

Interfacial stress–strain fields become complicated in thermal barrier coatings (TBCs) under cyclic thermal loading, which affects the stability and spalling failure of TBCs directly. The convex and concave interfacial structures of TBCs were approximated as a multilayer cylinder model, and an analytical method of TBCs for shakedown analysis was established. A series of 8-YSZ TBC specimens were prepared by the plasma spraying process, followed by isothermal and thermal shock tests. The results showed that the stability limit is significantly greater than the elastic limit, the limit for the convex model was higher than that in the concave model, the first failure occurs in the concave area, and the main failure mode of a thermal barrier coating is the appearance of cracks at the interface layer during a thermal shock test. For the coating samples prepared in this study, the stability limits were between 950 °C and 1050 °C, and the validity of the stability limit analysis model of a multilayer structure was verified.

1. Introduction

Thermal barrier coatings (TBCs) are widely used in turbines for propulsion and power generation. The ultimate objective of TBC applications is an improvement in fuel efficiency by increasing the operating temperature. This protective coating is also used to prolong the life of turbine components and to improve the reliability of engine components [1,2,3]. There are four primary constituents in a thermal barrier coating system [4]. They comprise (i) the thermal barrier coating (Top coat, TC) itself; (ii) the super-alloy component, treated here as the substrate (Sub); (iii) an aluminum-containing bond coat (BC) located between the Sub and the TC; and (iv) a thermally grown oxide (TGO), predominantly α-alumina, that forms between the TC and the BC. Because of the multilayer features and undulated interfacial microstructures, a complex stress and strain state appears around the interfacial area under thermal loadings, which induces instability in the system and eventually affects the system’s durability [5,6,7,8].

When a structure is subject to cyclic loads, three structural responses may occur: (i) purely elastic deformation, when the magnitude of the load applied is lower than the elastic limit of the structure; (ii) shakedown status, when the load is larger than the elastic limit of the structure but less than the critical limit-shakedown limit; and (iii) non-shakedown status, when the load applied is higher than the shakedown limit and non-restricted plastic flow occurs. If the non-shakedown condition (iii) occurs, the structure may undergo failure as either an incremental collapse (ratcheting) or alternating plasticity (low cycle fatigue). Generally, there are two ways to calculate the shakedown limit: incremental elastic–plastic analysis and shakedown analysis. The incremental elastic–plastic method is based on a calculation of the full history of the stress–strain curve, which involves a large computation cost. The method with the shakedown analysis is a direct approach to calculating the limit load and the stability of the structure. Classical shakedown theorems were initially proposed just for elastic and perfectly plastic materials. Due to the diversity and complex plastic behavior of engineering materials, these shakedown theorems should be extended before application. The stability analysis for a structure under thermal loading is affected by the following factors, a change in the stress field caused by thermal strain, a change in the yield strength and elastic parameters with temperature, and the creep and relaxation under elevated temperatures [9,10]. Therefore, for a TBC stability analysis under cyclic thermal load, considering the limitations of classical stability theory on materials and the complexity of the stress and strain field under cyclic thermal load, the classical stability theorems need to be developed. Additionally, the influence of temperature on the material parameters and yield strength should be considered.

For the relatively simple structures, the calculation methods and extended formulations of stability analysis are given in the literature [11]. The linear programming method was adopted by most researchers to solve the stability limit based on the essence of the stability theorem. Considering the temperature-dependent yield strength, the stability analysis for symmetrical cylindrical and spherical shell structures under an external force and temperature was calculated using the linear programming method in the literature [12]. Chen [13] analyzed the upper and lower limits of a few examples with a temperature-dependent yield strength using the linear matching method. For cylinder structures, Liu [14] performed a stability analysis under the combination of internal pressure and thermal loading, considering the temperature-dependent yield strength. Du [15] proposed a kinematic shakedown theorem derived from the concept of an assumed fictitious yield surface and with the temperature implicitly contained in the theorem. Then, Du [16] presented a new kinematic shakedown theorem without the concept of an assumed yield surface and with the temperature explicitly contained in the new formulation; the relationship between yield strength and temperature was also simplified as being linear. The influence of temperature on yield strength was obviously at the elevated temperature stage but had little effect at low temperatures, which indicates that approximating this relationship as being linear in the whole temperature range may lead to a significant error. For the stability of TBCs, some researchers [6,17,18,19,20,21,22,23] have proposed the concept of “displacement instability” and implemented experimental and simulation analyses on TBCs under cyclic thermal loading; the results showed that oxide growth is the main factor that induces displacement instability. Moreover, the main unstable mechanisms are the ratcheting effect and a volume decrease in the BC layer. Other researchers have also made explored TBC displacement instability. For example, Xu [24] utilized the elastic–plastic finite element method to analyze the stability of TBCs with different irregular interfacial morphologies, and the influence of the yield stress of the BC layer on coating stability was analyzed. Huang [25] performed a finite element analysis of interface displacement instability for the Fecralloy substrate with human-made grooves embedded into its surface, based upon FEA software, and the simulation result showed great agreement with the experimental result. TBC stability is not only related to the interface stability but also closely connected to the whole structure’s performance.

In this study, an analytical method for performing a shakedown analysis of TBCs is proposed. To verify its effectiveness and the rationality of the results, the convex and concave microstructures near the TBC interface are approximately treated as multi-cylinder models, and 8YSZ thermal barrier coatings are prepared by the plasma spraying process. The isothermal thermal shock test is carried out. The comparative test results show that the shakedown limit can easily be solved using the shakedown limit analysis method based on the multi-layer cylinder model and can easily be used for engineering applications. The stability limit of the thermal barrier coating is higher than the elastic design value, and the stability limit of the convex region of the interface is higher than that of the concave region. The structure first fails in the concave region. The larger the base curvature and coating thickness of the cylinder model, the higher the stability limit of the structure. The analysis results are consistent with the existing test results.

2. Shakedown Analysis of Cylinder Model

2.1. Cylinder Model Derived from the Interfacial Region of TBCs

The morphology of complex concave and convex shapes around the TBC interface can be approximately treated as a semicircle, as shown in Figure 1a. Extracting those local characteristics, the cylinder geometry model can be set up. The concave and convex models were built for areas A1 and A2, respectively. The material of each layer is illustrated as shown in Figure 1b. Generally, the thicknesses of TC, BC, and TGO range from 120 to 200 $\mathsf{\mu }\mathrm{m}$, from 30 to 200 $\mathsf{\mu }\mathrm{m}$, and from 2 to 10 $\mathsf{\mu }\mathrm{m}$, respectively [6].

2.2. Shakedown Analysis of Multilayer Cylinder Model

For an elasto-perfectly plastic solid body, the volume $V$, traction ${\widehat{p}}_i$, and displacement ${\widehat{u}}_i$ are given on boundaries $S_P$ and $S_u$. The body force is represented by $Q_i$. The basic governing equations and boundary conditions are as follows:

$$\begin{array}{c}{\sigma }_{ij,j}+Q_i=0 (\mathrm{in} V)\\ {\epsilon }_{ij}=\left(u_{i,j}+u_{j,i}\right)/2 (\mathrm{in} V)\\ {\sigma }_{ij}{n}_j={\widehat{p}}_i (\mathrm{on} S_P)\\ u_i={\widehat{u}}_i (\mathrm{on} S_u)\end{array}$$

(1)

Considering the geometry and working conditions for TBCs, such as no traction and displacement constraints on the outside surface, the body force, displacement, and strain continuity conditions on the interface can be neglected. Assume that the volume of each layer is $V_k$, that the surface of each layer is ${\left(S_u\right)}_k$, and that the displacement and strain are ${\left({\widehat{u}}_i\right)}_k$ and ${\left({\widehat{\epsilon }}_i\right)}_k$. Additionally, k = 1, 2, 3, 4, which stands for different layers in TBCs, as shown in Figure 1b. Then, the basic governing equations and boundary conditions for TBCs are as follows:

$$\begin{array}{c}{\left({\sigma }_{ij,j}\right)}_k=0 (\mathrm{in} V_k)\\ {\left({\epsilon }_{ij}\right)}_k=\left({\left(u_{i,j}\right)}_k+{\left(u_{j,i}\right)}_k\right)/2 (\mathrm{in} V_k)\\ {\left({\epsilon }_{ij}\right)}_k={\left({\widehat{\epsilon }}_{ij}\right)}_k (\mathrm{on} {\left(S_u\right)}_k)\\ {\left(u_i\right)}_k={\left({\widehat{u}}_i\right)}_k (\mathrm{on} {\left(S_u\right)}_k)\end{array}$$

(2)

Taking one layer in the multilayer cylinder model as an analysis object, elastic strain ${\epsilon }_{ij}^E$ and thermal strain ${\epsilon }_{ij}^T$ are given by

$${\epsilon }_{ij}^E=A_{ijkl}{\sigma }_{kl},{\epsilon }_{ij}^T=a_{ij}T$$

(3)

where $A_{ijkl}$ is a symmetric and positive definite elastic tensor, $a_{ij}$ is a tensor that stands for thermal expansion coefficient, and T stands for temperature.

Ignoring the viscosity factor during deformation, the total strain ${\epsilon }_{ij}$ is the sum of three components in the following form:

$${\epsilon }_{ij}={\epsilon }_{ij}^E+{\epsilon }_{ij}^P+{\epsilon }_{ij}^T$$

(4)

where ${\epsilon }_{ij}^P$ is the plastic strain. The yield surface equation is given by

$$F=f\left({\sigma }_{ij}\right)-K\left(T\right)=0$$

(5)

where $f\left({\sigma }_{ij}\right)$ is the stress field and $K\left(T\right)$ is the yield stress.

The flow rules are given by [15]

$$\{\begin{array}{l}{\dot{\epsilon }}_{ij}^P=\dot{\lambda }\partial F/\partial {\sigma }_{ij}\hfill \\ \dot{\lambda }\ge 0,if F=0 and \frac{\partial f}{\partial {\sigma }_{ij}}{\dot{\sigma }}_{ij}-\frac{dK}{dT}\dot{T}=0\hfill \\ \dot{\lambda }=0,if F<0 or F=0 and \frac{\partial f}{\partial {\sigma }_{ij}}{\dot{\sigma }}_{ij}-\frac{dK}{dT}\dot{T}<0\hfill \end{array}$$

(6)

where $\dot{\lambda }$ is the plastic multiplier ratio. Then,

$$\dot{\lambda }\ge 0,\dot{\lambda }F=0$$

(7)

The plastic power is defined as

$$D={\sigma }_{ij}{\dot{\epsilon }}_{ij}^P=D\left({\dot{\epsilon }}_{ij}^P,T\right)$$

(8)

Define

$$\{\begin{array}{l}K\left(T\right)=K_{0}h\left(T\right)\hfill \\ K_0=K\left(0\right)\hfill \end{array}$$

(9)

where $h\left(T\right)$ is a function related to temperature and $h\left(0\right)$ = 1.

From Equations (6)–(9), we have

$$D=D_0\left({\dot{\epsilon }}_{ij}^P\right)h\left(T\right)$$

(10)

where $D_0$ is the specific dissipation rate when T = 0.

The total stress ${\sigma }_{ij}^{}$ can be divided into two parts:

$${\sigma }_{ij}={\sigma }_{ij}^E+{\rho }_{ij}$$

(11)

where ${\sigma }_{ij}^E$ is the thermal elastic stress when the structure is treated as an ideal elastic solid under the same external loading and temperature. ${\rho }_{ij}$ is the residual stress, which is self-equalization stress and only determined by the instantaneous plastic strain field in the structure.

$${\rho }_{ij}=L\left({\epsilon }_{ij}^P\right)$$

(12)

where $L$ is a given linear integral operator.

For simplicity, only consider the quasi-static process here. Assuming that loads and temperature fields are linearized forms of load-temperature coefficients $a_s$ (s = 1, 2, …, n_s), then

$${\sigma }_{ij}^E\left(r,t\right)={\displaystyle \sum _{s=1}^{n_s}{a}_s\left(t\right){\sigma }_{ij}^{ES}\left(r\right)}$$

(13)

where ${\sigma }_{ij}^{ES}\left(x\right)$ is the thermal elastic stress field when $a_s=1$, $a_j=0\left(j\ne s\right)$, $r$ is spatial coordinate, and t is the time variable. Thus, total strain ${\epsilon }_{ij}$ can also be written as

$${\epsilon }_{ij}=A_{ijkl}{\sigma }_{kl}^E+{\epsilon }_{ij}^T+{\epsilon }_{ij}^P+A_{ijkl}{\rho }_{kl}$$

(14)

where

$$A_{ijkl}{\sigma }_{kl}^E+{\epsilon }_{ij}^T=\left(u_{i,j}^E+u_{j,i}^E\right)/2$$

(15)

$${\epsilon }_{ij}^P+A_{ijkl}{\rho }_{kl}=\left(u_{i,j}^R+u_{ij,i}^R\right)/2$$

(16)

$$u_i=u_i^E+u_i^R$$

(17)

where $u_i^E$ is the thermal elastic displacement field and $u_i^R$ is the residual displacement.

$u_i^E$ can be expressed as

$$u_i^E\left(r,t\right)={\displaystyle \sum _{s=1}^n{a}_s\left(t\right)u_i^{ES}\left(r\right)}$$

(18)

where $u_i^{ES}\left(r\right)$ is the thermal elastic displacement field when $a_s=1$, $a_j=0\left(j\ne s\right)$.

Employing this description, the change region of the coefficient $a_s$ decides the change limit of external loads and temperature.

$$p_i\left(r,t\right)={\displaystyle \sum _{s=1}^n{a}_s\left(t\right)p_i^s\left(r\right)},r\in S_p$$

(19)

$$Q_i\left(r,t\right)={\displaystyle \sum _{s=1}^n{a}_s\left(t\right)Q_i^s\left(r\right)},r\in V$$

(20)

$$T\left(r,t\right)={\displaystyle \sum _{s=1}^n{a}_s\left(t\right)T^s\left(r\right)},r\in V$$

(21)

Assuming that ${\dot{\epsilon }}_{ij}^{*P}$ is a kinematically admissible plastic strain rate in a cycle. Then, the plastic strain increment is given by

$$\Delta {\dot{\epsilon }}_{ij}^{*P}={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\dot{\epsilon }}_{ij}^{*P}dt=\frac{1}{2}\left(u_{i,j}^{*}+u_{j,i}^{*}\right)}$$

(22)

where $u_i^{*}=0$(on $S_u$).

The kinematic shakedown theorem can be described as follows: if there exists a load-temperature route $a_s\left(t\right)\in G$ and a plastic strain rate in the period cycle ${\dot{\epsilon }}_{ij}^{*P}$, such that Equation (22) is satisfied and [16]

$$\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_0+{\mathrm{\tau }}_s}\left\{{\displaystyle \sum _{s=1}^n{a}_s\left(t\right)\left[{\displaystyle {\int }_V{Q}_i^s{\dot{u}}_i^{*}dV+{\displaystyle {\int }_{S_P}{\widehat{p}}_i^s{\dot{u}}_i^{*}dS}}\right]+{\displaystyle {\int }_V{a}_{ij}T{\dot{\rho }}_{ij}^{*}dV}}\right\}dt}\\ >{\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_{V}D\left({\dot{\epsilon }}_{ij}^{*P},T\right)dVdt}}\end{array}$$

(23)

Then, shakedown will not occur. Additionally, G is the variation range of $a_s$ in n-dimensional space. ${\rho }_{ij}^{*}$ denotes the residual stress field associated with the plastic strain field ${\dot{\epsilon }}_{ij}^{*P}$.

Equation (23) is further simplified. Additionally, Equation (13) is rewritten:

$${\sigma }_{ij}^E={\sigma }_{ij}^{EE}+{\rho }_{ij}^T={\displaystyle \sum _{s=1}^{n_s}{a}_s\left(t\right)\left[{\sigma }_{ij}^{EEs}\left(x\right)+{\rho }_{ij}^{Ts}\left(x\right)\right]}$$

(24)

where ${\sigma }_{ij}^{EEs}$ is the stress, and it is only determined by external loads ${\widehat{p}}_i$ and $Q_i$. The elastic stress field ${\sigma }_{ij}^{EE}$ is determined if $a_s=1$, $a_j=0\left(j\ne s\right)$. ${\rho }_{ij}^{Ts}$ equals the residual stress field caused by temperature $T$ when $a_s=1$, $a_j=0\left(j\ne s\right)$.

With the principle of virtual power, the third term in the left-hand side of Equation (23) can be rewritten as

$${\displaystyle {\int }_V{a}_{ij}T{\dot{\rho }}_{ij}^{*}dV={\displaystyle {\int }_V{\epsilon }_{ij}^T{\dot{\rho }}_{ij}^{*}dV}}=-{\displaystyle {\int }_V{\rho }_{ij}^T{A}_{ijkl}{\dot{\rho }}_{kl}^{*}dV}={\displaystyle {\int }_V{\rho }_{ij}^T{\dot{\epsilon }}_{ij}^{*P}dV}$$

(25)

In Equation (23), the preceding two terms in the left-hand side can be transformed. Taking into account Equations (8), (24) and (25),

$${\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_V{\sigma }_{ij}^E\left(r,t\right){\dot{\epsilon }}_{ij}^{*P}\left(r,t\right)dVdt>}}{\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_{V}D\left({\dot{\epsilon }}_{ij}^{*P},t\right)}}dVdt$$

(26)

According to Equation (26), the shakedown criteria for four-layer TBCs can be obtained, namely, four inequations similar to Equation (26), corresponding to TC, TGO, BC, and Sub.

Considering that $h\left(T\right)$ is the linear function of T for simplicity, we let

$$h\left(T\right)=A-BT$$

(27)

where A and B are material constants; thus, we have

$$D\left({\dot{\epsilon }}_{ij}^{*P},T\right)=D_0\left({\dot{\epsilon }}_{ij}^{*P}\right)h\left(T\right)=A{D}_0\left({\dot{\epsilon }}_{ij}^{*P}\right)-BT{D}_0\left({\dot{\epsilon }}_{ij}^{*P}\right)$$

(28)

From the shakedown theorem above, the necessary condition for the shakedown state is immediately obtained. If shakedown occurs, then

$$\mathrm{\eta }{\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_V\left[{\sigma }_{ij}^E\left(r,t\right){\dot{\epsilon }}_{ij}^{*P}\left(r,t\right)+BT{D}_0\left({\dot{\epsilon }}_{ij}^{*P}\right)\right]dVdt}}\le {\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_{V}A{D}_0\left({\dot{\epsilon }}_{ij}^{*P}\right)dVdt}}$$

(29)

where $\mathrm{\eta }$ is the safety factor and $\mathrm{\eta }\ge 1$.

Next, the time integration is solved by the incremental failure criterion.

According to Equation (22),

$$\Delta {\dot{\epsilon }}_{ij}^{*P}\left(r\right)={\epsilon }_{ij}^{*P}\left(r,t_0+\mathrm{\tau }\right)-{\epsilon }_{ij}^{*P}\left(r\right)$$

(30)

According to the process and concept of incremental failure, there exists a scalar function $\lambda \left(r,t\right)$ such that

$${\epsilon }_{ij}^{*P}\left(r,t\right)=\lambda \left(r,t\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)$$

(31)

where

$$\dot{\lambda }\left(r,t\right)\ge 0, \lambda \left(r,t_0\right)=0,\lambda \left(r,t_0+\mathrm{\tau }\right)=1$$

(32)

Noticing that the specific dissipation rate $D_0\left({\dot{\epsilon }}_{ij}^{*P}\right)$ is the homogenous function of plastic strain rate ${\dot{\epsilon }}_{ij}^{*P}$, from Equations (31) and (32), we obtain

$${\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{D}_0\left({\dot{\epsilon }}_{ij}^{*P}\right)dt={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{D}_0\left(\dot{\lambda }\Delta {\epsilon }_{ij}^{*P}\right)}}dt=D_0\left(\Delta {\epsilon }_{ij}^{*P}\right)$$

(33)

From Equations (13), (21) and (23), taking $\mathrm{\eta }=1$, Equation (29) can be written as

$${\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_V{\displaystyle \sum _{s=1}^n{a}_s\left(t\right)}\left[{\sigma }_{ij}^{Es}\left(r\right){\dot{\epsilon }}_{ij}^{*P}\left(r,t\right)+B{T}^s\left(r\right)D_0\left({\dot{\epsilon }}_{ij}^{*P}\right)\right]dVdt}}\le {\displaystyle {\int }_{V}A{D}_0\left(\Delta {\epsilon }_{ij}^{*P}\right)dV}$$

(34)

Taking into account Equations (31) and (32), and $D_0\left(\dot{\lambda }\Delta {\epsilon }_{ij}^{*P}\right)=\dot{\lambda }{D}_0\left(\Delta {\epsilon }_{ij}^{*P}\right)$, the left-hand side of Equation (34) can be written as

$${\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\displaystyle {\int }_V{\displaystyle \sum _{s=1}^n{a}_s\left(t\right)}\left[{\sigma }_{ij}^{Es}\left(r\right)\dot{\lambda }\left(r,t\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)+B{T}^s\left(r\right)\dot{\lambda }\left(r,t\right)D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)\right]dVdt}}$$

(35)

For the given incremental failure mechanism $\Delta {\epsilon }_{ij}^{*P}\left(x\right)$, the function $\lambda \left(r,t\right)$ is chosen by the following method. In region G of $a_s\left(t\right)$, only if

$${\displaystyle \sum _{s=1}^n{a}_s\left(t\right)}\left[{\sigma }_{ij}^{Es}\left(r\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)+B{T}^s\left(r\right)\dot{\lambda }\left(r,t\right)D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)\right]$$

(36)

achieves admissible maximum can we have $\dot{\lambda }\left(x,t\right)\ne 0$. Now, the left-hand side of Equation (34) can achieve its maximum. Letting

$$f\left(r\right)=\underset{a_s\in G}{\max }\left\{{\displaystyle \sum _{s=1}^n{a}_s\left(t\right)}\left[{\sigma }_{ij}^{Es}\left(r\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)+B{T}^s\left(r\right)D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)\right]\right\}$$

(37)

be equal to $D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)$, the incremental failure criterion is obtained. Lastly, the incremental failure criterion can be expressed in the following form:

$${\displaystyle {\int }_{V}f\left(r\right)dV={\displaystyle {\int }_{V}A{D}_0\left(\Delta {\epsilon }_{ij}^{*P}\right)}}dV$$

(38)

The variation range $G$ of $a_s$ is defined as a group of inequalities:

$$a_s^{-}\le a_s^{}\le a_s^{+}$$

(39)

Then, Equation (38) can be written as

$${\displaystyle {\int }_V{\displaystyle \sum _{s=1}^n{\mathrm{\beta }}_s\left(r\right)y_s\left(r\right)}dV=}{\displaystyle {\int }_{V}A{D}_0\left(\Delta {\mathrm{\epsilon }}_{ij}^{*P}\right)}dV$$

(40)

where ${\mathrm{\beta }}_s\left(r\right)$ and $y_s\left(r\right)$ are given by

$$\{\begin{array}{l}{\mathrm{\beta }}_s\left(r\right)=\{\begin{array}{l}{a}_s^{+},y_s\left(r\right)>0\\ a_s^{-},y_s\left(r\right)<0\end{array}\hfill \\ y_s\left(r\right)={\mathrm{\sigma }}_{ij}^{Es}\left(r\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)+B{T}^s\left(r\right)D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)\hfill \end{array}$$

(41)

Then, the shakedown limit can be obtained from Equations (40) and (41).

3. Shakedown Analysis for Typical TBCs

3.1. Material Model and Loading

The bilinear plastic properties for BC are listed in

Table 1. The bilinear plastic date for BC [26,27].

T/°C	400	400	600	600	800	800	900	900	1000	1000
σ/MPa	1100	2500	1100	2200	300	380	45	60	10	15
εp	0	0.24	0	0.30	0	0.02	0	0.02	0	0.01

, which appears to be temperature-dependent [26,27]. The yield strength of TGO remains at 10 GPa between 0 and 800 °C but decreases linearly to 2 GPa at 1000 °C [28]. Temperature shows less effects on yield strength for BC and TGO under low temperatures, but the yield strength is significantly reduced above 800 °C. Therefore, it is not accurate to simplify the relationship between yield strength and temperature as being linear in the whole load range. The TC and Sub layers are treated as elastic, and the specified parameters can be found in the literature [29]. Considering the temperature-dependent property of TBCs, the various physical parameters are converted to the function of external load T by curve fitting, so that the analytic mathematical model derived above can be used for the stability limit solution. The radiuses of the interface from inside to outside are set to 5 mm, 5.13 mm, 5.135 mm, and 5.165 mm.

The thermal load is applied to the surface of TC, and the cool gas flows inside the Sub to reduce the system temperature. The thermal stress calculating condition is set at 1000 °C on TC and at 25 °C on the Sub inner side.

Therefore, for a load-temperature coefficient $a_s$, the thermal load ranges from 25 °C to 1000 °C and $\Delta T$ lies within 0–975 °C, which corresponds to $a_s$ being 0–1.

Assuming that the material of TGO and BC conforms to the Tresca yield criterion, the yield stress decreases as the temperature increases.

$$\max \left\{\left|{\mathrm{\sigma }}_r\right|,\left|{\mathrm{\sigma }}_{\mathrm{\phi }}\right|,|{\mathrm{\sigma }}_{\mathrm{\phi }}-{\mathrm{\sigma }}_r|\right\}=2K\left(T\right)$$

(42)

where

$$K\left(T\right)=K_0\left(A-BT\right), \left(2K_0={\mathrm{\sigma }}_s\right)$$

(43)

3.2. Bilinear Model for Yield Strength and Temperature

The radial velocity for small deformation can be described as $\dot{u}\left(r\right)=\frac{\dot{c}}{r}$.

When the incremental collapse is studied, we have

$$\{\begin{array}{l}\Delta {\mathrm{\epsilon }}_r={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\dot{\epsilon }}_{r}dt}={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}\frac{\partial \dot{u}}{\partial r}dt}=-\frac{\Delta c}{r^2}\hfill \\ \Delta {\mathrm{\epsilon }}_{\phi }={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}{\dot{\epsilon }}_{\mathrm{\phi }}dt}={\displaystyle {\int }_{t_0}^{t_0+\mathrm{\tau }}\frac{\dot{u}}{r}dt}=\frac{\Delta c}{r^2}\hfill \end{array}$$

(44)

where $\Delta c>0$. $D_0\left(\Delta {\dot{\epsilon }}_{ij}^{*P}\right)$ refers to the above plastic deformation state and can be expressed as

$$D_0={\mathrm{\sigma }}_r\left(-\frac{\Delta c}{r^2}\right)+{\mathrm{\sigma }}_{\mathrm{\phi }}\left(\frac{\Delta c}{r^2}\right)=\frac{\Delta c}{r^2}\left({\mathrm{\sigma }}_{\mathrm{\phi }}-{\mathrm{\sigma }}_r\right)=2K_0\frac{\Delta c}{r^2}$$

(45)

It is necessary to evaluate ${\mathrm{\sigma }}_{ij}^{Es}\left(x\right)\Delta {\epsilon }_{ij}^{*P}\left(x\right)$ in $y_s\left(r\right)$ to calculate the left-hand side in Equation (40).

$${\mathrm{\sigma }}_{ij}^{Es}\left(r\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)={\mathrm{\sigma }}_r^E\left(r\right)\Delta {\mathrm{\epsilon }}_r\left(r\right)+{\mathrm{\sigma }}_{\mathrm{\phi }}^E\left(r\right)\Delta {\mathrm{\epsilon }}_{\mathrm{\phi }}\left(r\right)$$

(46)

where ${\mathrm{\sigma }}_r^E\left(r\right)$ is the radial stress and ${\mathrm{\sigma }}_{\mathrm{\phi }}^E\left(r\right)$ is the circumferential stress, which can be obtained from the literature [30,31].

Since $y_s\left(r\right)={\mathrm{\sigma }}_{ij}^{Es}\left(r\right)\Delta {\epsilon }_{ij}^{*P}\left(r\right)+B{T}^s\left(r\right)D_0\left(\Delta {\epsilon }_{ij}^{*P}\left(r\right)\right)>0$ is always satisfied, ${\mathrm{\beta }}_s\left(r\right)=1$.

Transforming Equation (40), we obtain

$$A_D={\displaystyle {\int }_V{D}_0\left(\Delta {\dot{\epsilon }}_{ij}^{*P}\right)}dV-{\displaystyle {\int }_V{\displaystyle \sum _{s=1}^n{\mathrm{\beta }}_s\left(r\right)y_s\left(r\right)}dV}={\displaystyle {\int }_{a_i}^{a_{i+1}}{D}_0\left(\Delta {\dot{\epsilon }}_{ij}^{*P}\right)rdr}-{\displaystyle {\int }_{a_i}^{a_{i+1}}{\displaystyle \sum _{s=1}^n{\mathrm{\beta }}_s\left(r\right)y_s\left(r\right)}rdr}$$

(47)

The expression for each variable in Equation (47) can be obtained, and substituting thermal load $T_w$ into Equation (47), a function with only one variable ($T_w$) can be derived. Due to A_D being characterized as a decreasing function, when the temperature is set initially at 25 °C, the function’s value is greater than zero. Therefore, $A_D$ decreases with increasing $T_w$, and when $A_D$ is equal to zero, the corresponding temperature is the stability limit.

3.3. Stability Analysis Result

Combined with the traditional elastic design method, the stability limits of the concave and convex models for TBCs were analyzed with different radius curvatures and layer thicknesses based on the cylinder model.

3.3.1. Effect of Curvature on Stability Limit

The effects of the Sub radius curvature on the elastic limit and the stability limit for the concave and convex model are shown in

Table 2. Limit load for TBCs with different curvatures of Sub.

Radius of Curvature/mm	Concave Elastic Limit Load/°C	Concave Stability Limit Load/°C	Convex Elastic Limit Load/°C	Convex Stability Limit Load/°C
15.165	605	797.8	722	949.3
17.165	609	799.8	724	950.1
19.165	615	802.8	726	950.9
21.165	623	806.8	728	951.6
23.165	635	812.8	730	952.3

. The stability limits are greater than the elastic limit for both the concave and convex models. The elastic limit and stability limit increase with increasing radius curvature. The elastic and stability limits for the convex model are higher than those in the concave model under the same radius curvature. With the development of research and manufacturing levels for TBCs, the actual working temperature has reached more than 1000 °C; however, the stability limit obtained in this study is less than the actual temperature, as shown in Table 2 (just between 800 and 900 °C). It is mainly due to stress release and relaxation under high temperatures as an influence of material creep, which releases system stress effectively. However, creep is not considered in this study, which results in the calculated thermal stresses being higher than those in actual conditions and having smaller stability limits.

3.3.2. Effect of Layer Thickness on Stability Limit

The effects of layer thickness on elastic limit and stability limit for the concave and convex model are shown in

Table 3. Relationship between limit load and thickness of TBCs for the concave model.

Thickness Incremental Ratio	TC Thickening Elastic Limit Load/°C	BC Thickening Elastic Limit Load/°C	TGO Thickening Elastic Limit load/°C	TC Thickening Stability Limit Load/°C	BC Thickening Stability Limit Load/°C	TGO Thickening Stability Limit Load/°C
0	605	605	605	797.8	797.8	797.8
0.1	622	605.1	605.4	805.3	798	798.14
0.2	636	605.2	606	813.6	798.2	798.4
0.3	648	605.3	606.5	822.7	798.35	798.7

and

Table 4. Relationship between limit load and thickness of TBCs for the convex model.

Thickness Incremental Ratio	TC Thickening Elastic Limit Load/°C	BC Thickening Elastic Limit Load/°C	TGO Thickening Elastic Limit load/°C	TC Thickening Stability Limit Load/°C	BC Thickening Stability Limit Load/°C	TGO Thickening Stability Limit Load/°C
0	730	730	730	949.3	949.3	949.3
0.1	743	731.2	731.9	952.8	949.5	950.1
0.2	754	732.4	733.5	956.4	949.7	950.9
0.3	763	733.6	735.5	960	949.9	951.4

, respectively. First, the sensitive relationship between the layer thickness and elastic limit, and stability limit can be obtained, and the effect of thickness of the TC layer on the elastic and stability limits is greater than that of the BC layer thickness variation. Second, the elastic limit and stability limit increase with a thickening layer, and the slope of the stability limit curve decreases gradually, which indicates that, with an increase in TC thickness, the incremental limit load of TBCs decrease. Even with the TC thickness increasing to the threshold, the incremental value is very small. This indicates that trying the limit load cannot be improved by increasing the thickness of TC.

Hence, the elastic and stability limits for the convex model are higher than those in the concave model, as shown in Table 2, Table 3 and Table 4. Namely, compared with the outer convex model, the concave model has less ultimate bearing capacity, and the TBCs are usually damaged at this position, which is in agreement with previous experimental results, as shown in Figure 1. The crack initiation around the concave interface after a long working time leads to spalling failure.

4. Constant Temperature Oxidation Test and Thermal Shock Test

To verify the rationality of the stability limit analysis model and the results of multilayer structures, according to the testing methods for thermal cycle and thermal shock resistance of thermal barrier coatings (standard NO: JIS H8451-2008), a series of 8-YSZ TBC specimens were prepared by the plasma spraying process and combined with the constant temperature oxidation test and thermal shock test. The TBC specimens were stable when the thermal shock temperatures were 850 °C and 950 °C, and the structure can be considered to be in a stable state. When the thermal shock temperature was 1050 °C, the stability of the TBC samples was significantly reduced and the stability limits were between 950 °C and 1050 °C.

4.1. Specimen Preparation

Classical 8-YSZ TBC specimens were prepared, and the dimensions and materials for the TBCs are shown in

Table 5. The dimensions and materials for the specimens.

	SUB	BC	TC
Material	DZ 125	Ni-Cr-Al-Yttrium	8 wt%Y2O3-ZrO2
Radius/mm	10	10	10
Thickness/mm	2	0.1	0.25

.

The specimen preparation process is shown in Figure 2, and the specimen after deposition is shown in Figure 3.

4.2. Testing Plan and Processing

For the constant temperature oxidation test, 24 TBC samples were heated to 1000 °C in an SRJX-8-13 box-type resistance furnace, and then, at 1 h, 3 h, 6 h, 12 h, 24 h, 48 h, 72 h, and 96 h, three samples were taken out each time and air-cooled. To study the stability of TBCs under a thermal cycle, thermal shock tests were carried out at different temperatures based on the failure evaluation of coat peeling. The thermal shock temperatures were 850 °C (three samples, numbered as T1-1, T1-2, and T1-3), 950 °C (three samples, numbered as T2-1, T2-2, and T2-3), and 1050 °C (five samples, numbered as T3-1, T3-2, T3-3, T3-4, and T3-5), for a total of 11 samples. The heat preservation time of each cycle was 1 h, followed by air cooling. Based on the number of thermal shocks and the thermal shock temperature used [32,33], the upper limit of the thermal shock number was set as 300 times; after that, if still no cracking or spalling failure occurs, we argue that the coating structure is stable. When the thermal shock temperature is 1050 °C, a sample is taken out at the 5th and 30th hours to investigate the difference in interface deformation and to verify the stability of interface displacement.

4.3. Test Results and Verification

4.3.1. Verification of Interface Stability

The TGO layer formed at the TC/BC interface and was distributed in a strip at the beginning of the thermal cycle at 1050 °C. After 30 thermal cycles, despite the TC layer not spalling, micro-cracks appeared in the coating internal. Once the micro-cracks were connected and extended, an interface crack appeared. As the heat cycle continued, the sample started to peel [32], failures were located at the TC/TGO and TGO/BC interfaces, and the cracking direction was parallel to the interface direction. The failure mechanism was determined to be in accordance with the conclusions of similar literature [33], which verified previous conclusions about the cracking location of the coating.

The TBCs prepared by plasma spraying were not flat at the initial SUB/BC/TC interfaces. Under thermal shock conditions, the TGO layer formed and continuously thickened at the BC/TC interface due to BC oxidation. At the same time, the unfitness of the interface was further aggravated (the amplitude of the curve interface changes). The interface of the TBCs showed significant instability (displacement instability) until internal cracks appeared and further expanded the coating, eventually leading to peeling failure of the coating. Since the thermal cycle test is a process with alternating cold and heat, the alternating stress generated in the thermal cycle process makes it easier for the TC layer and TGO layer to generate defects such as holes and micro-cracks. After several thermal cycles, defects such as holes in the coating gradually grow and connect to form small cracks. The micro-cracks further expand and connect under the actions of stress and penetrate through the coating, resulting in peeling of the coating. During thermal shock, some cracks appear at the TC/TGO/BC interface, and the failure can be observed after the peeling layer is removed. This is mainly due to the presence of fluctuation in the interface topography of TC/TGO/BC. Under thermal cycling conditions, BC layer oxidation, thermal mismatch, phase transition, and other factors cause a very complex stress state near the interface. Especially for the convex part of the interface morphology, its thermal stress state is more complex. Under the condition of alternating heat and cold during the thermal cycling process, the stress concentration in the above parts is more likely to be generated, and the residual stresses are usually large, which further promote the generation and extension of the crack defect in the coating until the final peeling failure. Therefore, the interfacial oxidation growth brought by the thermal cycle will cause a deformation in the coating interface to reduce the stability and then will lead to coating peeling failure, which is consistent with the above conclusions about the overall stability of the coating and the interfacial stability.

4.3.2. Verification of the Analysis Method for Structural Stability Limit

The stability limit range of TBCs can be verified by a thermal shock test. As shown in

Table 6. Failure thermal cycles of coatings at different thermal shock temperatures.

Sample	850 °C			950 °C			1050 °C
	T1-1	T1-2	T1-3	T2-1	T2-2	T2-3	T3-1	T3-2	T3-3
Thermal cycles	300 *)	300 *)	300 *)	300 *)	154	300 *)	300 *)	77	139

Note: *⁾ No peeling or cracking was observed after 300 thermal cycles.

, due to the relatively low-temperature oxidation at 850 °C, the three samples numbered T1-1, T1-2, and T1-3 show no significant peeling or cracks on the coating surface after 300 thermal shocks. At the thermal shock temperature of 950 °C, no obvious coating spalling or significant cracks appear in samples T2-1 and T2-3 after 300 thermal shocks. However, sample T2-2 endured thermal shocks only 154 times before cracks appear on the coating surface. At the thermal shock temperature of 1050 °C, sample numbered T3-1 showed no obvious coating peeling or significant cracks after 300 thermal shocks, but significant cracks appeared on the coating surface after 77 thermal shocks for sample T3-2 and after 139 thermal shocks for sample T3-3.

From the test results, the effects of thermal shock temperature on the failure of TBCs are more obvious at the high-temperature stage (such as 1050 °C). For the thermal shock test conditions, the TBCs have a thermal stability limit. Additionally, for the thermal shock temperatures 850 °C, 950 °C and 1050 °C, regardless of the coating preparation, test conditions, and precision of the instrument, the temperatures 850 °C and 950 °C can be assumed to be stable temperatures, while the thermal barrier coating exhibited poor stability under thermal cycling at 1050 °C. Therefore, it can be considered that the stability limit of the thermal barrier coating prepared in this study is between 950 °C and 1050 °C, slightly higher than the results obtained from the structural stability limit analysis model in the previous section (800–950 °C). This is mainly due to the simplification of some factors in the stability limit analytical model, which causes specific errors. Although the thermal shock tests are only carried out at three temperature conditions, it verified the validity of the stability limit analysis model presented in this study.

5. Conclusions

Considering the characteristics of the microstructure interface morphology of TBCs, an analytical method for multiplayer structure shakedown analysis was established, and the isothermal and thermal shock tests for 8-YSZ TBCs specimen were carried out to verify the validity of the proposed method. The main conclusions included the following:

For classical TBCs, the stability limit is significantly greater than the elastic limit, and the limit for the convex model was higher than that in the concave model. Additionally, the larger the radius of curvature of the Sub and the thickness of the TBCs, the higher the stability limit. This indicates that the first failure occurs in the concave area, and a smooth interface for the coating is beneficial to the stability of the multiplayer structure. According to the experimental results, (i) the interface oxidation at high temperature would cause interface instability, which will further affect the overall stability and lead to coat peeling failure; (ii) he mathematical approach developed for the stability limit evaluation is reliable and reasonable; (iii) the stability of the TBCs was significantly reduced with the increase in temperature; and (iv) for the coating sample prepared this time, when the thermal shock temperatures were 850 °C and 950 °C, the TBC sample had good stability and the structure could be considered to be in a stable state, and when the thermal shock temperature was 1050 °C, the stability of the TBC samples was significantly reduced, and the shakedown limit was between 950 °C and 1050 °C, which verifies the effectiveness of the shakedown limit analysis model of the multi-layer structures proposed in this paper.