Discrete Element Modeling of Fracture Behavior and Stress Analyses in Thermal Barrier Coatings During Wear Tests Open Access

Thermal barrier coating (TBC) systems typically consist of a dual-layer structure, comprising a ceramic top layer (TC) that protects the metal from high temperatures, excessive wear, and corrosion, and a metallic bond coat (BC) that safeguards the metal substrate from oxidation and corrosion while enhancing adhesion to the ceramic layer. TBCs are widely used in the hot components of gas turbine engines (e.g., turbine blades and combustion chambers) due to their low thermal conductivity and high resistance to thermal corrosion [1]. During engine operation, thermal barrier coatings are subject to thinning due to factors such as high-temperature friction and sand particle erosion, resulting in a significant decline in their thermal protection capabilities, which, in severe cases, may lead to catastrophic failures [2]. Therefore, a comprehensive investigation into the wear mechanisms of thermal barrier coatings is of significant practical importance and engineering value for enhancing the service life and reliability of gas turbine engines and preventing catastrophic failures.

Thermal barrier coatings have been successfully applied in wear and friction control for mechanical components and turbine blades, protecting high-temperature alloy substrates from mechanical fatigue stresses and foreign particle erosion [3]. Previous studies have primarily focused on assessing the microstructure and wear resistance of TBCs. Research indicates that the tribological behavior of thermal barrier coatings is influenced by three key factors [4,5]: (1) plasma spraying parameters, such as spraying distance and plasma power; (2) wear testing conditions, including sliding speed, applied load, lubricant medium, and temperature; (3) coating microstructure, such as multilayer and nanostructured coatings. Yunus and Alsoufi employed a Taguchi L16 factorial mixed-level experimental design to assess the wear-rate and coefficient of friction (COF) of ceramic coatings applied via atmospheric plasma spraying. They also utilized Taguchi-based gray relational analysis to optimize the process parameters for various tribological characteristics of these coatings [6]. Borgaonkar and Syed investigated the impact of various input parameters (coating deposition process, coating material properties, and thickness) on the rolling contact fatigue life and tribological performance of coated rolling and sliding contact elements, and summarized potential combinations of input parameters that can enhance the performance of coated contact elements [7]. Chen et al. conducted friction and wear experiments on coatings using a wear testing machine, revealing that nanostructured thermal barrier coatings exhibit enhanced resistance to plastic deformation and superior wear resistance compared to traditional coatings, attributed to the optimization of their microstructure and improved mechanical properties [8]. Rakesh et al. utilized an air-jet erosion wear testing apparatus to investigate the wear performance of coatings sprayed at different substrate preheating temperatures. The study indicated that coatings exhibit higher scratch and wear resistance at lower preheating temperatures, and the pseudoelastic behavior of the austenitic phase (NiTi-B₂) plays a significant role in wear resistance [9]. Sun et al. used a ball-on-flat rotating friction and wear testing machine to compare different compositions of Al₂O₃ coatings, finding that under thermomechanical coupling, composite coatings with NiAl and Bi₂O₃ components developed a friction oxide layer on the wear surface, exhibiting superior lubrication and wear resistance [10]. Moridi et al. employed the finite element method (FEM) to study the effect of coating roughness on the stress distribution of the specimen [11]. Wu et al. deposited a series of Ni60 alloy coatings on copper surfaces at different preheating temperatures and found that as the preheating temperature increased, the dilution rate of the coatings significantly increased, resulting in compositional inhomogeneity and a decline in friction performance [12]. However, research on the friction damage mechanisms of coatings and their relationship with external factors, particularly in-depth analysis at the microscopic level, remains limited. Most existing studies primarily focus on the evaluation of macroscopic performance, while research on the microscopic damage evolution processes, stress distribution, and crack propagation of coatings under various external conditions (such as load and sliding speed) is still insufficient. Therefore, further investigation of the microscopic damage behavior of coatings under complex working conditions is crucial for optimizing coating design and enhancing their wear resistance.

Typically, wear behaviors occurring at the microscale, such as local structural detachment and debris dislodgement, are difficult to observe and analyze experimentally. Additionally, the FEM in continuum mechanics is not well suited for modeling noncontinuous phenomena such as the evolution of heterogeneous structures and repeated dynamic fracture. The discrete element method (DEM), a numerical simulation method that has been extensively utilized in analyzing materials with discontinuous media, offers distinct advantages in the investigation of structural damage in materials, particularly in the context of impact and friction. Li et al. used a discrete element model to simulate the effects of interface bond coat roughness and preexisting cracks on crack propagation induced by particle erosion in thermal barrier coatings. The results indicate that as the bond coat roughness increases, the propagation length of delamination cracks decreases, and the delamination cracks are suppressed. Initial vertical TC defects can effectively inhibit the nucleation and propagation of new cracks [13]. Najafabadi et al. employed DEM and analytical modeling to conduct a numerical simulation analysis of abrasive wear on iron ore pellets. Comparing the results with experimental data, the numerical analysis predicted errors in abrasive wear of iron ore pellets of +21.7% for optimized quality and −34.4% for low-quality core blocks [14]. Gao et al. utilized DEM to simulate particle interactions and, in conjunction with FEM equations, successfully investigated the effects of different normal loads on seismic moment, macroscopic friction coefficient, and sliding behavior in a shear particle fault system containing fault gouge and plates, demonstrating the capability of DEM in studying the stick-slip dynamics of particle fault systems [15].

This study focuses on the DEM simulation of the wear process of thermal barrier coatings. First, the DEM parameters for the abrasives and coatings were calibrated and validated using the experimentally measured macroscopic properties. Second, the wear behavior of the TBCs was investigated in depth through the simulated wear test and was compared with the ball-on-flat wear experiment when possible. The damage processes, crack formation, and energy aspects of the wear system were evaluated. The Hertzian contact theory was employed to explore the wear mechanisms. The coefficients of friction of the coatings under various normal loads, sliding speeds, and number of cycles were computed. Finally, the effects of normal load and the number of cycles on the wear behavior of the coatings were evaluated.

This study employed the particle flow code (PFC) program for modeling. The PFC employs a bidirectional coupling iterative algorithm based on Newtonian dynamics and contact constitutive relations to numerically simulate the mechanical behavior of the discrete medium. At each computational time-step, the system first updates the kinematic parameters (such as acceleration, velocity, and spatial position) of the particles based on Newton's second law, using the resultant force and torque acting on the particles. Subsequently, based on the updated particle displacements, the normal and tangential stiffness of the contact surface, as well as the relative displacement, are computed through the force–displacement relation, allowing the determination of the contact force and the identification of contact failure states. The aforementioned mechanical parameters are driven by the time-steps, forming a dynamic iterative process. The computation loop terminates when the contact force residual satisfies the preset convergence criteria. This algorithm achieves a self-consistent solution of the macroscopic and microscopic mechanical responses of the particle system through the iterative update of the microscopic parameters. In the two-dimensional (2D) model, particle shapes are represented as 2D discs, and the material constitutive relationships are established through microscale contacts and bonding between particles [16]. The forces and moments acting at the contact points are related to the maximum normal and shear stresses in the surrounding bonded material; if any of these maximum stresses exceed the corresponding limit bond strength, the bonded material, along with its associated forces, moments, and stiffness, will be removed from the model.

Different constitutive models result in various mechanical properties. In this work, the simulated ceramic layer material lanthanum zirconate (LZO), being a brittle material, was modeled using a parallel bond model, while the bond layer material NiCrAlY and the high-temperature alloy substrate GH4169 (Inconel 718), being ductile materials, were modeled using a contact bond model. Additionally, the contacts between particles in the model were randomly categorized into strong and weak contacts to better simulate the behavior of real composite materials. The models were implemented in the PFC program.

In this study, three different initial values for each microparameter were established as control groups, arranged in an orthogonal layout for numerical experiments. The orthogonally arranged numerical experimental groups of the TC layer are detailed in Table 1. Each group of microvalues was then substituted into four calibration tests, and the corresponding macrovalues for each group were calculated based on the principles outlined in Sec. 2.1. Subsequently, a fitting equation based on the correlation between microparameters and macroparameters was used to determine the final values of the microparameters.

It is considered that a Pearson correlation coefficient $r≥0.3$ indicates a significant correlation between the two parameters, and the Pearson correlation coefficients for the TC layer are presented in Table 2.

Finally, the obtained microparameter values were substituted back into the calibration tests to calculate the macroparameters. The values of the macroparameters were then observed, and the microparameter values were fine-tuned accordingly to obtain the final values of the microparameters. The final microstructural parameters of the TC layer and abrader are shown in Table 3, while those of the BC layer and the substrate are presented in Table 4.

Following the aforementioned steps, the microparameters and macroparameters obtained from the numerical simulation tests and experiments for the coating system are presented in Table 5.

The calibration of the interlayer strength of the thermal barrier coating followed the ASTM C633-2001 international standard. The specific two-dimensional discrete element numerical experiment is shown in Fig. 2. The simplification principle lay in replacing the adhesive bonding force with the contact force between the designated wall and particles, while the loading speed of the loading block was substituted by the assigned wall velocity. Additionally, the bonding strength within the discrete element coating model was determined by the microscopic strength values (⁠$Tσ$⁠, $Sσ$⁠) assigned between the coating-bonding layer and the bonding layer-substrate (i.e., different types of particles). In this study, the microscopic strength values $Tσ$ and $Sσ$ were set to 200 MPa.

To reduce errors, three repeated tests were conducted on the specimens in the numerical experiment, with the bonding strength taken as the average value from the three tests. Table 6 presents the tensile test results for the NiCrAlY bonding layer system of the thermal barrier coating. The average bonding strength measured in the numerical experiment was 28.73 MPa, while the actual value measured in the experiment was 29.69 MPa, resulting in an error of only 3.233%.

This simulation established a 2D DEM model for the wear process of thermal barrier coatings based on the data in Table 5. The model consists of two parts: the abrader and the thermal barrier coating system. The abrader used is tungsten carbide (WC), which has the strength and stiffness significantly greater than that of the ceramic layer: E = 550 GPa, v = 0.3, $σb=510MPa$⁠, $σbc=5GPa$⁠, $KIC=6.5MPa⋅m1/2$⁠, and the ball radius $Rc=0.15mm$⁠.

The overall dimensions of the thermal barrier coating system are 3.2 mm in length and 1.2 mm in height, with the TC layer measuring 0.25 mm in height and the BC layer measuring 0.15 mm. The entire model consists of 29,255 particles, with a total of 76,795 contacts between the particles. The particle size in the TC layer is 8–10 µm, in the BC layer is 5–6 µm, and in the substrate is 5–7.5 µm [21–23]. The principle of the numerical simulation friction test was the same as that of the laboratory friction and wear experiments, as shown in Fig. 3. First, the longitudinal and lateral displacements of the substrate are constrained by three walls; second, the abrader is subjected to a normal load while being made to slide laterally in a reciprocating manner to achieve cyclic friction. As the number of cycles increases, the friction performance and wear condition of the thermal barrier coating are observed.

In this simulation, numerical friction tests were conducted at various sliding speeds with a normal load P ranging from 5 to 25 N. The sliding velocity V includes three levels of 0.10, 0.13, and 0.16 m/s.

In the DEM wear model, the wear debris appears in the form of individual particles and particle clusters. The vertical displacement of the particles indicates the wear depth of the surface. To investigate the abrasive wear behavior of the coating, displacements of the abrader under a load of 25 N at different numbers of cycles were selected, as shown in Fig. 4. In the figure, the dashed lines represent cracks, while the yellow particles denote wear debris. With the reciprocating motion of WC, concentrated stress and adhesive wear behavior occur at the surface contact due to tensile, compressive, and shear stresses, causing surrounding particles to undergo plastic deformation, resulting in spallation of the coating particles. The wear surface exhibits pitting, leading to the formation of numerous pits and consequently resulting in plenty of grooves, scratches, and ploughings, which render the fracture surface of the coating uneven.

As shown in Fig. 4(a), the initial wear of the coating is relatively mild, characterized mainly by microcracks formed from plastic deformation, with no spallation occurring. Figure 4(b) illustrates that the damage at the contact surface intensifies, with existing cracks and damaged areas leading to the onset of particle spallation and wear debris formation, resulting in the appearance of small pits on the coating surface. In the mid-friction stage, Fig. 4(c) shows that a large number of coating particles have spalled off and transformed into wear debris, with numerous pits present on the wear track surface; however, no deep mechanical ploughings were observed. In the late friction stage, Fig. 4(d) illustrates that the extrusion phenomenon caused by the motion of the abrader exacerbates the wear of the coating surface, with some wear debris filling the pits on the rough surface, while another portion leaves the damage zone with the movement of the abrader. While moving within the damage zone, asperities are formed by the accumulation of wear debris. After the wear debris and asperities are destroyed, the surface forms plenty of grooves, scratches, and ploughings, leading to further expansion of the damaged area.

The simulated damage morphology qualitatively agrees with experimental observations [24], as shown in Fig. 5. In the figure, the solid line represents the outer contour of the coating, while the dashed line indicates cracks. Since the strength and stiffness of the abrader are much greater than those of the ceramic layer, wear primarily occurs on the surface of the coating. In the initial stage of friction, the wear depth is shallow, with numerous microcracks accumulating beneath the contact interface between the abrader and the coating. As the number of cycles increases, the wear depth gradually increases, with some microcracks extending into the coating, forming macroscopic cracks. The spalled particles either contribute to the formation of wear debris within the wear layer or move away from the contact surface.

From Eq. (12), it can be seen that in contact pairs with a higher equivalent elastic modulus, the contact radius becomes smaller, resulting in higher tensile, compressive, and shear stresses. Thus, it can be inferred that the maximum stress leading to crack formation is greater at the contact edge. It can be observed from Eqs. (14)–(16) that the maximum stresses affecting crack initiation and propagation are inversely proportional to the square of the contact radius. From Eq. (18), the quantitative relationship between the contact load, material fracture toughness, and the dimensionless constant can be observed. Specifically, the crack length is positively correlated with the load but exhibits a nonlinear relationship; moreover, higher fracture toughness helps reduce crack propagation.

The simulation results show that the depth of the maximum shear stress in the coating is much lower than the thickness of the coating, and the shear stress distribution is entirely within the coating and does not significantly affect the substrate. However, as the load increases, the maximum tensile stress and maximum shear stress also increase, exacerbating wear.

For the maximum shear stress depth Z, when the abrader moves from left to right, the DEM simulation results are shown in Fig. 6. Since the software package could not directly identify the distribution of shear stress between particles, this study suggests that the maximum shear stress depth Z of the coating by measuring the maximum lateral interaction force between two particles using the software package.

As shown in Fig. 6, when the load is small, the maximum lateral force between particles at each load (indicated by the arrows) occurs at a distance of 15 µm from the coating surface. The numerical simulation results and theoretical calculation results for the maximum shear stress depth Z are shown in Fig. 7, and both results are consistent. As the load is increased, the lateral force at greater depths gradually increases (boxed area), with the maximum lateral force extending downward. When the load is increased to 25 N, the maximum lateral force extends downward to 30 µm, as shown in Fig. 6(e). The reason is that the simulation program can only identify the interaction force between two particles, and the minimum height between the two particles is 15 µm; therefore, the maximum shear stress depth Z varies in increments of 15 µm.

Figure 8 presents the localized cone cracks generated by the abrader sliding from left to right in the DEM wear simulation, along with experimental evidence of cone cracks produced in Ref. [26]. Among them, Fig. 8(a) depicts the crack damage map of the discrete element simulation model at a sliding speed of 0.10 m/s under a load of 20 N, while Fig. 8(b) shows a microscopic photo of the lateral sliding failure of soda-lime glass caused by a steel ball with a radius of 3.17 mm under a normal load of P = 20 N and a friction coefficient of 0.1. The selected scenario in Fig. 8 represents the condition where the coating experienced multiple cyclic movements of the abrader before a significant number of cracks appeared. Initially, these cycles did not result in noticeable crack formation on the coating surface. However, when the abrader moved from left to right in a subsequent cycle, a greater number of left-oriented cracks emerged. As shown in the schematic in Fig. 8, when the abrader slides from left to right, localized cracks oriented to the left (marked in black) are generated. Conversely, as the abrader slides back and forth, when it moves from right to left, cracks oriented to the right (marked in yellow) are generated. These cracks result from traces left by the abrader's right-to-left motion in previous cycles. This further confirms that crack formation is closely related to the reciprocating movement of the abrader. Ultimately, due to the reciprocating sliding motion of the abrader, cracks are formed in both directions and connect to create an inverted conical shape; these continuous localized cracks, oriented in opposite directions, intersect and generate an inverted conical damage region with each contact. The particles in the damaged area undergo plastic deformation during the sliding wear process, resulting in pitting and subsequent spallation.

In DEM wear simulations, according to wear criteria, when the parallel bonding exceeds the strength limit, the bond fails, and each contact failure generates a new crack that passes through the contact point between two particles. The direction of each crack is defined by the unit vector perpendicular to the line connecting the two particles, while the crack angle is calculated using the unit direction vector of the crack. As the number of cycles increases, more contact areas reach the critical strain limit, leading to crack initiation at these weak points and extending toward even weaker planes, ultimately resulting in the formation of primary cracks. The high stresses generated by continuous sliding contacts promote the interconnection of more cracks, leading to outward propagation and the formation of abraders that separate from the coating. Cracks arise from the failure of the bond; therefore, the energy of bond failure should be included in the discussion of the wear mechanism of the coating, as illustrated in Fig. 9.

The bond failure released bond failure energy, leading to the formation of cracks. Therefore, each crack generated represents the release of an equivalent amount of bond failure energy. The cumulative crack propagation rate followed the same growth trend as the accumulation of bond failure energy, meaning the cumulative crack propagation rate could be observed through the cumulative growth rate of bond failure energy. When the P value was set to different values, the growth trends of the cumulative crack propagation rate and bond failure energy were similar, resembling an “S” shape. However, after stabilization, the maximum values differed for each case. Figure 9 illustrates the number of newly formed cracks per cycle and the generation of adhesive failure energy in the coating as the number of cycles increases, using P = 5 N as an example at a sliding speed of 0.10 m/s with the abrader. The numbers (I)–(III) correspond to the damage stages of the coating's cracks. As shown in the figure, the wear mechanism during the microcrack initiation stage (stage (I)) is characterized by slight abrasive wear, with fine cracks beginning to appear within or on the surface of the coating material. At this stage, crack propagation occurs at a relatively slow rate, corresponding to Fig. 4(a). In the accelerated crack growth stage (stage (II)), a large number of coating particles undergo spallation and turn into debris, with numerous pits observed on the wear track surface. The crack propagation rate increases significantly, often accompanied by localized spallation or fracture of the coating, corresponding to Figs. 4(b) and 4(c). In the stabilized damage stage (stage (III)), the coating surface develops plenty of grooves, scratches, and ploughings. The crack growth rate slows down once again. During this stage, cracks in the remaining coating propagate slowly under relatively low stress levels, typically after localized coating failure has already occurred. Damage accumulates progressively, but the propagation rate is lower than in the previous stage, corresponding to Fig. 4(d).

The COF is a dimensionless parameter that reflects the friction characteristics of one material sliding against another. Figure 10 illustrates the variation of the coating's coefficient of friction with increasing friction cycles under different normal loads. It was observed that during the initial friction stage, the COF for all materials tended to increase, likely due to enhanced interactions between the asperities of the two materials. Subsequently, the growth trend of the friction coefficients slowed down. If the simulation continued, it was expected that the COF would peak when the maximum adhesion force is reached and then fluctuate around a certain value (if the asperities began to recover after reaching maximum deformation, the COF would decrease), until a stable state was reached. This is a typical behavior associated with the stable state of tribological conditions [30].

Figure 11 illustrates the evolution of the average friction coefficient (COF) of the tribo-pair with respect to load and sliding speed and compares it with the ball-on-flat wear test results from Ref. [31] (test conditions: disc rotation speed of 250 r/min, normal loads of 8, 12, 16, 20, and 24 N) for validation. The results indicate that the simulation results are closest to the experimental results when the sliding speed is 0.10 m/s. At a sliding speed of 0.10 m/s, the COF at a 5 N load is the highest, approximately 0.58; however, as the applied load is increased, the COF gradually decreases to 0.51. Analyzing the reason, the compacted friction layer acts as an effective solid lubricant, resulting in a lower coefficient of friction. Furthermore, the increase in load leads to an increase in shear stress in the contact area, making it easier for the friction layer to undergo plastic deformation under high loads, thus reducing surface roughness. The findings of this study corroborate the aforementioned explanation.

As illustrated in Fig. 11, the COF for the friction pair is less than 0.58, with an overall variation range between 0.47 and 0.58. Under varying working conditions, the COF exhibits slight variations, with an average decrease of 11.92% and 5.93% observed as load and abrader sliding speed are increased, respectively. When the applied load is held constant, the COF gradually decreases with the increasing sliding speed, which is consistent with the trends reported in Refs. [32–34]. It is believed that this decrease may result from the higher frictional velocity, which potentially induces microscopic plastic deformation or wear on the material surface. This deformation or wear reduces the true contact area of the surfaces in contact, consequently decreasing both the frictional force and the coefficient of friction. Alternatively, the increase in frictional velocity may facilitate the formation of a third-body layer of wear debris between the contact surfaces. The sliding behavior of this third-body layer can also lead to a reduction in the coefficient of friction.

As observed in Fig. 11, the variation trends of the COF influenced by changes in sliding speed and load exhibit significant similarities between simulations and experiments. Moreover, it confirms that the influence of sliding speed and normal load on the COF is consistent across both macroscale and microscale. Therefore, the DEM wear model can qualitatively investigate the effects of sliding velocity and normal load on the coating's COF.

In this work, a discrete element model has been developed to simulate the ball-on-flat wear test of the thermal barrier coatings. The main conclusions are summarized as follows:

1. During the wear test process, concentrated stress occurs at the contact surface of the coating due to the combined stresses. This leads to coating fragmentation and spallation of the surrounding material, resulting in pitting on the worn surface and the formation of multiple pits, which subsequently create grooves, scratches, and ploughings.

2. The maximum tensile stress occurs at the contact edge, while the shear stress distribution is entirely located within the coating. The DEM results are consistent with the Hertzian contact theory, in which the maximum tensile stress and maximum shear stress are inversely proportional to the square of the contact radius and directly proportional to the one-third power of the normal load.

3. As the friction test continues, cracks form and propagate, extending deeper along the thickness of the material. The energy generated by crack formation and bond failure is increased slowly during the initial stage of friction; however, due to the accumulation of damage and the intensification of bond fracture, the number of cracks and the energy of bond failure increased rapidly until it reaches a stable saturation value.

4. The study shows that the coefficient of friction generally decreases with increasing load and friction speed. However, under low-speed and low-load conditions, the instability in lubrication caused by the accumulation of wear debris may lead to fluctuations in the COF.

These preliminary results provide important support for the optimization of the wear resistance of thermal barrier coatings from a microscopic mechanism perspective, particularly in understanding wear behavior under high stress. A detailed analysis of stress distribution and crack propagation during friction offers a scientific basis for coating design, aiming to improve its performance stability under extreme temperatures and frictional environments. Meanwhile, the study identifies the variation of the friction coefficient with load and speed, laying the foundation for optimizing coating performance under practical application conditions. This helps in formulating effective operational strategies to reduce wear-rates and failure risks, thereby enhancing the overall reliability and service life of the coating.

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

• The Natural Science Foundation of Tianjin (21JCYBJC01400).

• 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.