Discrete Element Modeling of Fracture Behavior of Thermal Barrier Coatings Under Bending Conditions Open Access

Coating the surface of high-temperature alloy turbine blades with thermal barrier coatings (TBCs) can reduce the metal substrate temperature by 100–300 °C, significantly extending the blade's service life while enhancing engine efficiency and fuel economy [1]. In operational environments, TBCs are subjected to mechanical loads such as tension and bending, which can induce surface microcracks or delamination. These defects accelerate coating detachment and degradation. The fracture and delamination of the coating compromise its protective function, exposing the substrate material to high temperatures or corrosive media, ultimately increasing the risk of substrate failure [2]. Studying the damage behavior and fracture mechanisms of TBCs is essential for enhancing the service life and durability of components, as well as for gaining a deeper understanding of the damage patterns in coating systems.

In recent years, numerous researchers have focused on predicting the lifetime and understanding the failure mechanisms of TBCs. Studies have revealed that defects often make the interface a weak point in the mechanical performance of TBCs. The primary fracture modes observed include surface cracking, internal interface cracking within the coating, and cracking at the coating–substrate interface, particularly under high temperatures, erosion, or other extreme conditions. These fracture phenomena are largely attributed to residual stresses induced by oxidation or thermal cycling, as well as stress concentrations resulting from external loads. Previous studies [3–6] have discussed the effects of geometric and physical factors on coating failure caused by residual stresses. For example, Yu et al. investigated residual stresses in TBCs during the cooling process using numerical simulations and analyzed how the material properties of each layer influence these stresses [4]. Bhattacharyya et al. analyzed residual stress distribution in TBCs via plate bending analysis, comparing two material systems [5]. In addition, Refs. [7–10] provide studies on fracture modes induced by external loads. For instance, Li et al. analyzed the crack propagation behavior of TBCs by simulating particle erosion [7]. Zhang et al. conducted tensile tests on TBCs and studied the microscopic morphology of their cross-sectional fractures [8]. Šulák et al. conducted low cycle fatigue (LCF) tests in a symmetrical push-pull cycle under strain control at 900 °C. Cyclic hardening/softening curves, cyclic stress–strain curves, and fatigue life curves of both the TBCs-coated and uncoated materials were assessed [10].

Three-point bending (3PB) and four-point bending (4PB) tests are commonly used methods for studying the fracture behavior of brittle coatings, employed to analyze and evaluate the mechanical properties of the coating/substrate [11–14], damage characteristics [15–18], and influencing factors [19–21]. For example, Zhao et al. improved the 4PB specimen and combined it with finite element analysis to accurately evaluate the fracture toughness of the TBCs from a fracture mechanics perspective. Through experimental analysis and theoretical calculations, they derived the bond strength required for the initiation of interface edge cracks and the fracture toughness of the system [11]. Chen et al. employed a sandwiched four-point bend specimen to evaluate the crack growth resistance in plasma-sprayed TBCs, discovering a rising crack growth resistance curve and obtaining the steady-state strain energy release rate [14]. Wang et al. employed the acoustic emission (AE) technique to characterize the failure behavior of TBCs under three-point bending tests. The results indicated that the main failure mechanisms of TBCs under three-point bending conditions include the debonding of the metal coating from the substrate and the propagation of the horizontal crack along the bond coat (BC)/substrate interface under the action of the flexural moment [17]. Eriksson et al. conducted tensile bending fatigue tests on TBCs and observed that under low-load conditions, the TBCs exhibited extensive delamination damage, while under high-load conditions, the TBCs demonstrated relatively better integrity. At the same time, they used finite element (FE) analysis to investigate potential transverse crack structures [18]. Frommherz et al. measured the delamination toughness and its sintering behavior using an improved 4PB apparatus. The research results indicate that the increase in coating stiffness is related to changes in its microstructure. The specific spray structure of TBCs has a great influence on the sinter stability which results in a characteristic temperature dependency [19]. In addition, Martins et al. conducted tests on coatings using an improved 4PB test to study the effect of the morphology curvature of the bond coat on the fracture mechanical properties of TBCs [20]. Despite the extensive research on TBCs in the aforementioned studies, the relationship between the performance of these complex coatings and their microstructure remains unclear. The discrete element method (DEM), a widely used numerical simulation technique for discontinuous media materials, can effectively represent the discrete characteristics and interactions of bulk materials. This method, by simulating particle collisions and frictional interactions, can capture the behavior of materials under complex loading conditions, making it particularly suitable for performance prediction and microstructural analysis of coating materials.

To investigate the effects of localized damage and cracks on the overall mechanical behavior of coatings, this study employs DEM and conducts a four-point bending (4PB) numerical simulation by constructing a rectangular DEM model of the coating system. The failure modes of TBCs under bending loads were analyzed at the microscale, and the simulation results were validated against existing experimental data. Additionally, the study examines two coatings with different thicknesses to evaluate the influence of coating geometry on structural mechanical properties and crack propagation behavior. The TBCs analyzed consist of a ceramic top layer (TC) made of lanthanum zirconate (LZO) and a bond coat (BC) composed of NiCrAlY, both prepared using atmospheric plasma spraying (APS), with a high-temperature alloy GH4169 (Inconel 718, Zhejiang Flysun Special Steel Co., Ltd, China) serving as the substrate.

The DEM model consists of a collection of discrete circular particles, where the contact forces between particles follow the force-displacement law, and the particles move based on Newton's second law. As a DEM, the bonded particle model includes a contact bond model (CBM) and a parallel bond model (PBM) [22]. Different constitutive models give rise to distinct mechanical properties in the medium. In this study, the ceramic layer material, LZO, is brittle, so the parallel bond model is used, whereas the bond coat material, NiCrAlY, and the high-temperature alloy substrate, GH4169 (Inconel 718), are plastic, and thus, the contact bond model is applied.

Since the DEM model reflects the macroscopic mechanical properties of the coating based on the mechanical properties of particles and their interaction behaviors, a systematic calibration of the model's microscopic parameters is required prior to numerical simulations to ensure accurate prediction of the macroscopic mechanical performance of TBCs. The symbols and definitions of the parameters involved in the calibration are detailed in Table 1. During the calibration process, initial values are assigned to various microscopic parameters. Numerical simulations of typical mechanical behaviors, including uniaxial compression, three-point bending, fracture toughness, and uniaxial tension, are conducted for each functional layer of the coating, as shown in Fig. 1. Subsequently, based on the obtained stress–strain curves and load–displacement curves, the equivalent macroscopic mechanical parameters of the coating are extracted, following the calculation methods provided in Eqs. (1)–(5).

On this basis, the correlation between microscopic and macroscopic parameters is analyzed using Eq. (6), taking the relationship between the microscopic and macroscopic parameters of the TC layer in Table 2 as an example. Key microscopic parameters that are highly correlated with the target macroscopic properties are identified and fitted using multivariate linear regression to establish multivariate linear regression equations between microscopic parameters and macroscopic mechanical properties (e.g., Eqs. (7)–(11) represent the regression models for the TC layer). The final microscopic parameters obtained from the regression analysis are subsequently incorporated into the DEM numerical model for validation simulations. The resulting macroscopic parameters from the simulations are compared with experimental data for each material to evaluate the rationality of the calibrated parameters and the accuracy of the model predictions.

To systematically investigate the effect of coating thickness on the mechanical properties and damage behavior of the TBCs, two different coating thicknesses were designed in this study. The first group (G1) has a TC layer thickness of 300 µm, a BC layer thickness of 150 µm, and a substrate thickness of 750 µm; the second group (G2) has a TC layer thickness of 600 µm, a BC layer thickness of 150 µm, and a substrate thickness of 750 µm. The length of both sets of samples is 12 mm. The Particle Flow Code (PFC) software was used to build DEM models for simulating the 4PB test. The model consists of a collection of particles bonded together, with the particle positions randomly distributed in an irregular pattern. The G1 model contains 83,587 particles, while the G2 model contains 93,890 particles. In the simulation, the upper end of the model is fixed by a support disk, while the lower end is moved upward at a speed of 16 µm/s by a loading disk to replicate the 4PB test process. The schematic of the loading and support arrangement of the model is shown in Fig. 2.

Figure 3 shows the crack propagation and fracture process of G1 during the 4PB test. Multiple cracks formed during the bending process, with the main crack located near the loading disk. The enlarged view in Fig. 3 further reveals the details of a typical crack. Initially, when the stress on the TC layer surface exceeds its fracture limit, transverse cracks perpendicular to the interface appear at the interface. Due to the sufficiently high fracture toughness of the BC layer, these cracks are blocked at the TC/BC interface (Fig. 3(a)), corresponding to point $A(A′)$ on the curve in Fig. 9.

Then, as the bending continues, new transverse cracks initiate on the surface of the TC layer. At the same time, the deflection tendency of the transverse cracks gradually increases, leading to the formation of interface microcracks at the TC/BC interface (Fig. 3(b)), corresponding to point $B(B′)$ on the curve in Fig. 9. However, such transverse cracks do not propagate before reaching steady-state expansion conditions. In such a way, multiple transverse cracks rapidly propagate (some transverse cracks initiate from the TC/BC interface and extend upward to the surface of the TC layer), causing further damage (Fig. 3(c)), corresponding to point $C(C′)$ on the curve in Fig. 9.

When the transverse cracks in the TC layer approach saturation, interface cracks at the TC/BC interface begin to expand rapidly. Additionally, transverse cracks initiate from this interface and propagate through the BC layer, eventually being blocked at the BC/substrate interface (Fig. 3(d)), corresponding to point $D(D′)$ on the curve in Fig. 9. Then, steady-state expansion of the interface cracks at the TC/BC interface leads to the convergence of these cracks, resulting in delamination at the TC/BC interface (Fig. 3(e)), corresponding to point $E(E′)$ on the curve in Fig. 9. Finally, when the transverse cracks in the BC layer reach the BC/substrate interface, the cracks deflect into this interface and begin propagating along it, eventually causing delamination between the BC layer and the substrate (Fig. 3(f)), corresponding to point $F(F′)$ on the curve in Fig. 9.

For G2, the cracking process is similar to that of G1. However, compared to G1, G2 cracks at lower deflection values and strain rates (Fig. 11). By comparing the fracture conditions of two coating thicknesses under the same strain conditions (Fig. 4), it can be observed that the thick coating exhibits a wider TC/BC interface delamination zone and a narrower BC/substrate interface delamination zone compared to the thin coating. Furthermore, the results show that as the TC layer thickness increases, the delamination location shifts toward the TC/BC interface, which is similar to the experimental results of Jiang et al. [27].

As previously mentioned, when there are no pore defects in the coating, cracks propagate along a direction perpendicular to the induced maximum tensile stress (Fig. 5(a)). However, due to the presence of microcracks and pores ahead of the main crack, in the TC layer (Fig. 5(b)), cracks nucleate at the valley regions of the top coat surface and around the preexisting pores. When the crack extends to an existing pore, phenomena such as crack penetration through the pore, crack path deflection, and crack bifurcation (indicated by the arrows) occur. The crack propagation behavior is consistent with the experimental observations reported in Ref. [31] (Fig. 5(c)). The crack branching mechanism is a complex phenomenon, influenced by various factors, including the crack tip velocity, stress intensity factor, and their rates [32].

During the 4PB test, the overall particle contact state of the TBCs transitions from tensile stress in the center region of the coating to compressive stress in the shoulder region (near the end constraints). The force chain distribution characteristics of G1 are shown in Fig. 6. Initially, tensile contact primarily occurs on the upper surface of the TC layer and the lower surface of the substrate, while the shoulder region and BC layer are mainly in a compressed state (Fig. 6(a)).

When the transverse cracks in the TC layer nucleate, the contact state at the crack changes to compression, and tensile contact moves outward from the crack, concentrating near the BC/substrate interface (Fig. 6(b)), corresponding to Fig. 3(a). As the transverse cracks in the TC layer approach saturation, the contact state in the BC layer changes to tensile, with tensile contact primarily concentrated in the BC layer (Fig. 6(c)), corresponding to Figs. 3(b) and 3(c). When transverse cracks nucleate in the BC layer, tensile contact in both the TC and BC layers decreases, with tensile contact concentrated near the BC/substrate interface (Fig. 6(d)), corresponding to Fig. 3(d). When delamination occurs at the TC/BC interface, the contact state of the TC layer changes to compression, and tensile contact in the BC layer decreases (Fig. 6(e)), corresponding to Fig. 3(e). When delamination occurs at the BC/substrate interface, a significant amount of tensile contact appears in the fractured delaminated part of the TC layer, while tensile contact is primarily concentrated near the substrate at the BC/substrate interface (Fig. 6(f)), corresponding to Fig. 3(f).

Additionally, it is worth noting that the study found that before crack initiation, the particles at the crack initiation location were all in tensile contact, and crack formation was primarily caused by tensile instability fracture. Tensile failure dominates the damage mechanism. At the same time, the shoulder region (end constraints) remains in a compressed contact state, with no significant damage observed. Therefore, it can be inferred that the coating exhibits strong damage resistance under compressive conditions. For G2, similar evolutionary behavior and the same damage pattern are observed under the 4PB load. This is because, in both coatings, the damage mechanism is dominated by tensile failure.

Figure 7 presents the stress–strain curve for the elastic deformation phase in this simulation. For G1, the yield strain during the elastic deformation phase is 0.05%, the yield strength is 187 MPa, and the bending modulus is 0.37 MPa. In contrast, G2 exhibits a relatively smaller yield strain of 0.04% during the elastic deformation phase, with a corresponding yield strength of 185 MPa and a bending modulus of 0.46 MPa. The results indicate that the thicker coating has a higher bending modulus, demonstrating greater resistance to deformation, whereas the thinner coating may experience a substrate effect, where part of the stress is borne by the substrate, leading to relatively lower stress under the same strain conditions.

From Eq. (24), it can be seen that the displacement change is decomposed into elastic response $(ΔF/ke)$ and damage displacement $(Δud)$⁠, thereby quantitatively distinguishing between elastic deformation and nonlinear damage effects. By calculating the variation of damage displacement with load, the accumulation pattern of material damage, the evolution trend of nonlinear behavior, and the characteristics of stiffness degradation and energy dissipation can be revealed.

Figure 9 shows the normalized load–displacement curves obtained from the G1 and G2 models. To facilitate comparison with the data from Ref. [37] and make the curves more intuitive, the load–displacement curves were normalized according to Eqs. (25) and (26). All the curves in the figure show an initial linear range, which then changed to a nonlinear portion, indicative of an elastic-plastic response. In the $OA(OA′)$ stage, as the load applied by the loading disk increases, the TBCs exhibits linear deformation, and the accumulated strain energy provides the energy required for subsequent crack initiation. When the load reaches point $A(A′)$⁠, it reaches the critical value, and the initial nonlinear “toe” region can be observed. In $AB(A′B′)$ stage, the load decreases with increasing displacement. This may be due to cracks in the TC layer being saturated and then blocked at the TC/BC interface, resulting in the instantaneous release of accumulated strain energy [38]. Subsequently, in $BE(B′E′)$ stage, as the loading displacement further increases, the load continues to increase. Finally, in $EF(E′F′)$ stage, delamination occurs at the TC/BC interface and the BC/substrate interface, leading to the release of strain energy again. As displacement increases, the load growth slows down.

Figure 10 presents the fracture modes of the TBCs under 4PB loading, and the fracture characteristics are consistent with those in Ref. [39] (Fig. 10(a)). Under a strain rate of 0.6%, the crack density of the thin coating G1 (Fig. 10(b)) is significantly higher than that of the thick coating G2 (Fig. 10(c)). Additionally, it can be observed that the spacing between transverse cracks is approximately one–two times the coating thickness. These values align well with the shear lag model used to determine crack spacing in thin films on ductile substrates [40]. According to this model, the typical crack spacing $lc$ should satisfy $3<(σYSlc)/(σbChc)<23$⁠, where $hc$ is the coating thickness, $σYS$ is the yield strength of the substrate, and $σbC$ is the tensile strength of the coating. Using $σYS=718.5MPa$ from Ref. [41] and $σbC=413.6MPa$ as calibrated in Table 4 (with the tensile strength of the TC layer being much smaller than that of the BC layer, it is neglected), the relationship $1<lc/hc<2$ is obtained, which is in good agreement with the observed values from this simulation.

Figure 11 shows the relationship between the crack density and strain in the TBCs, where crack density is defined as the ratio of the number of transverse cracks to their maximum distribution distance. From the figure, it can be seen that the crack density in G1 is generally higher than that in G2, but during the initial stage of crack nucleation, the crack density in G2 is slightly higher than in G1. The growth trend of crack density is in good agreement with the 4PB results for YSZ TBCs presented in Ref. [42], and the values are similar. Moreover, for cracks in the BC layer, G2 produces cracks earlier than G1 at lower strain rates. Combining the comparative analysis in Fig. 3, it is found that once the strain reaches 0.5%, further increasing the strain does not significantly increase the number of cracks in the crack density but instead more easily leads to delamination of the coating. This result is consistent with the experimental observations in Ref. [35], indicating that delamination of the coating is the primary failure mode at high strain stages.

The AE characteristics of TBCs are related closely to the formation of microcracks. In PFC, the formation of a microcrack is accompanied by a release of strain energy, which causes an AE event. The AE characteristics can be simulated from the statistics of microcrack numbers over a certain interval during loading [43]. Figure 12 shows the relationship between strain energy and AE counts as a function of strain in the TBCs. The results indicate that as the bending strain of the specimen increases, the strain energy of the material continuously accumulates, with part of the strain energy being converted into the energy required for crack propagation. Crack propagation is accompanied by local damage, leading to the formation of crack surfaces and changes in the internal stress state of the material, thereby affecting the release and transfer of strain energy. During this process, the rate of crack propagation is closely related to the applied load and the accumulation rate of the material's strain energy.

Furthermore, the experiment further revealed the significant impact of coating thickness on crack formation and strain energy variation. The thick coating G2, due to its larger volume and complex structural characteristics, exhibits significantly higher strain energy and AE counts than the thin coating G1 under the same strain conditions. This is primarily because the thick coating can withstand larger stress concentration areas during loading, thereby accumulating more elastic energy and generating more crack sources.

Through simulation experiments and theoretical analysis, the damage evolution behavior and mechanical properties of TBCs with two different thicknesses under 4PB tests were clarified. The study confirmed that the DEM model can effectively simulate the mechanical response of the TBCs, providing a reliable tool for coating performance prediction and design optimization. The main conclusions are summarized as follows:

1. Under continuous bending, the damage process exhibits six stages: The transverse cracks first initiate and propagate in the TC layer. The cracks then multiply and become saturated; local interface cracks propagate between the TC and BC layers. Transverse cracks initiate and propagate in the BC layer. Local interface cracks propagate between the BC layer and the substrate; Finally, when some adjacent interface cracks connect with each other, interface delamination occurs.

2. During the bending process, the overall particle contact state of the TBCs transitions from the tensile region at the center of the coating to the compressive region at the shoulder area (near the end constraints). Both coating thicknesses experience tensile instability fracture during the bending process, with tensile failure being the dominant mechanism of coating damage.

3. This study systematically analyzes the relationship between load and displacement variation based on experimental load–displacement curves, exploring the underlying mechanisms and the relationship between different loading stages and the damage evolution of TBCs. The study also indicates that coating thickness has a significant impact on mechanical properties: thicker coatings have higher load-bearing capacity, but they also carry a greater risk of crack accumulation.

4. By observing the fracture modes of the TBCs under 4PB, it was found that the crack density in the thin coating is significantly higher than that in the thick coating, with the spacing between transverse cracks being approximately one–two times the coating thickness. Furthermore, the thick coating generates cracks earlier than the thin coating at lower strain rates. The average crack spacing in the thick coating is larger than that in the thin coating, and its damage evolution rate is faster, exhibiting more pronounced interface delamination. Based on the comparative analysis of damage evolution, it can be observed that when the strain reaches 0.5%, further increasing the strain does not significantly increase the crack density but rather leads to easier delamination of the coating. These results indicate that coating delamination is the primary failure mode at high strain stages.

5. By analyzing the variations of strain energy and AE counts with strain in the TBCs, it was found that as the bending strain increases, the strain energy of the material continuously accumulates and is partially converted into crack propagation energy. Crack propagation is accompanied by local damage, which alters the internal stress state and affects the release of strain energy. Compared to the thin coating G1, the thick coating G2 exhibits higher strain energy and AE counts under the same strain conditions. The strain energy release rate is closely related to the crack propagation rate, and when the release rate exceeds the critical value, cracks propagate rapidly. The thick coating G2 has a higher critical release rate, allowing it to release more energy during propagation while maintaining higher strain energy storage, demonstrating stronger energy-bearing capacity.

This study employs the DEM to analyze, at the microscale, the damage evolution mechanism of TBCs under 4PB conditions, providing new theoretical insights into the failure modes of TBCs. The study reveals the significant impact of coating thickness on damage behavior and mechanical properties, particularly in terms of differences in crack propagation and interface delamination damage evolution. This research provides a theoretical basis for optimizing coating design, enhancing coating load-bearing capacity, and predicting coating failure, offering important academic and engineering significance for the application and performance enhancement of TBCs.

Yafeng Li greatly appreciates the support from the Natural Science Foundation of Tianjin (21JCYBJC01400) and the China Scholarship Council (CSC) Scholarship.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.