Open Access

ARTICLE

# The Coupled Thermo-Chemo-Mechanical Peridynamics for ZrB_{2} Ceramics Ablation Behavior

1 Department of Engineering Structure and Mechanics, Wuhan University of Technology, Wuhan, 430070, China

2 Hubei Key Laboratory of Theory and Application of Advanced Materials Mechanics, Wuhan University of Technology, Wuhan, 430070, China

* Corresponding Author: Qiwen Liu. Email:

(This article belongs to this Special Issue: Peridynamics and its Current Progress)

*Computer Modeling in Engineering & Sciences* **2023**, *135*(1), 417-439. https://doi.org/10.32604/cmes.2022.021258

**Received** 05 January 2022; **Accepted** 19 May 2022; **Issue published** 29 September 2022

## Abstract

The ablation of ultra-high-temperature ceramics (UTHCs) is a complex physicochemical process including mechanical behavior, temperature effect, and chemical reactions. In order to realize the structural optimization and functional design of ultra-high temperature ceramics, a coupled thermo-chemo-mechanical bond-based peridynamics (PD) model is proposed based on the ZrB_{2}ceramics oxidation kinetics model and coupled thermo-mechanical bond-based peridynamics. Compared with the traditional coupled thermo-mechanical model, the proposed model considers the influence of chemical reaction process on the ablation resistance of ceramic materials. In order to verify the reliability of the proposed model, the thermo-mechanical coupling model, damage model and oxidation kinetic model are established respectively to investigate the applicability of the proposed model proposed in dealing with thermo-mechanical coupling, crack propagation, and chemical reaction, and the results show that the model is reliable. Finally, the coupled thermo-mechanical model and coupled thermo-chemo-mechanical model are used to simulate the crack propagation process of the plate under the thermal shock load, and the results show that the oxide layer plays a good role in preventing heat transfer and protecting the internal materials. Based on the PD fully coupled thermo-mechanical model, this paper innovatively introduces the oxidation kinetic model to analyze the influence of parameter changes caused by oxide layer growth and chemical growth strain on the thermal protection ability of ceramics. The proposed model provides an effective simulation technology for the structural design of UTHCs.

## Keywords

Ultra-high-temperature ceramics (UHTCs) refer to a class of ceramic materials that can maintain physical and chemical stability under high-temperature and oxygen atmospheres. They are mainly multi-component composite ceramics composed of borides, carbides, nitrides and other transition metal compounds based on hafnium, zirconium, and molybdenum, such as ZrB2, HfB2, TaC, HfC, ZrC, and HfN [1]. UHTCs have the characteristics of ultra-high temperature resistance, high thermal conductivity, and high strength. They can be used as heat-bearing structural components such as the nose cone and wing leading edge of reusable spacecraft. In order to achieve the goal of lightweight and high performance of UHTCs, the optimal design of UHTCs ceramic materials must be realized based on the understanding of UHTCs ablation mechanism and ablation resistance.

The research on UHTCs began as early as 1960s [2], but it was difficult to obtain excellent performance due to the technical constraints at that time. In the late 1990s, the urgent demand for ultra-high-temperature non-ablative thermal protective materials has brought them under the limelight [3]. Among currently studied UHTCs, ZrB2-and HfB2-based ceramics render good anti-oxidation and ablation properties, and achieve long-term stability under ultra-high temperature, showing promise as non-ablative ultra-high-temperature heat protection materials. In 2007, Han et al. [4] have investigated the oxidation resistance of ZrB2-SiC composite under the oxyacetylene flame. At greater than 2200°C, the mass loss rate and linear oxidation rate of the composite after 10 min were found to be 0.23 mg/s and 0.66 μm/s, respectively. Also, a dense ZrO2 layer was formed on the oxidized surface, and no obvious cracks and denudation were observed. Then, in 2011, Wang et al. [5] have studied the ablation behavior of C/C-ZrC composites under the oxyacetylene flame and demonstrated that the oxidation reaction generated a dense ZrO2 layer when ZrC was exposed to a high-temperature oxygen environment. When the temperature is increased to 2800°C, the solid ZrO2 transforms into a molten state and covers the material surface. The formation of a dense ZrO2 layer prevents further ablation of C/C-ZrC composites. In 2020, Qing et al. [6] have studied the ablation mechanism of C/C-ZrC-SiC composites containing ZrB2-SiC multiphase ceramic rods. The results demonstrated that the ablation mechanism of Zr-Si-B-C multiphase ceramics in the temperature range of 25°C to 3000°C generates B2O3 film in the low-temperature zone, SiO2 film in the medium-temperature zone and ZrO2 film in the high-temperature zone. This Zr-Si-B-C multiphase ceramic ensures that the material will not be seriously ablated in the low-and high-temperature regimes. Based on these studies, it can be found that a liquid or solid oxide layer is formed on the surface of ablative UHTCs under a high-temperature aerobic environment, which can effectively inhibit or prevent the further reaction between external oxygen and ceramic matrix to reduce the amount of ablation.

At present, the ablation resistance mechanism of ZrB2 ceramic materials has been clearly understood, but the numerical theoretical research on ablation has not been fully studied. Considering the complexity of the actual service environment, a numerical model needs to be established to describe various complex physical and chemical reactions during the ablation process of UHTCs, and this numerical model can be used to study the factors affecting the ablation process, as well as the influence behavior and mathematical quantitative function of each factor against ablation, so as to predict the ablation rate, temperature distribution, structural deformation and denudation damage. These works are helpful to the structural optimization and functional design of UHTCs. In 2005, Fahrenholtz [7] has proposed a theoretical model that can calculate the oxygen partial pressure in the coexistence of ZrB2, ZrO2 and B2O3. The proposed model can also predict the mass fraction of liquid B2O3 remaining in porous solid ZrO2. The model was experimentally verified by Talmy et al. [8]. Then, Fahrenholtz [9] has analyzed the oxidation behavior of ZrB2-SiC UTHCs and proposed the reaction sequence, oxide composition and formation rate during the oxidation process. Based on these studies, Parthasarathy et al. [10,11] have proposed an oxidation kinetic model that can describe oxidation process of diboride UHTCs, such as ZrB2, HfB2 and TiB2. Obviously, the ultra-high-temperature ablation of ceramics needs to comprehensively consider the interaction of temperature field, structural deformation and chemical reactions. Zhou et al. [12] have extended the oxidation kinetics model of ZrB2 and proposed a thermal force coupling model, which could describe the multi-field coupling behavior of ZrB2 during high-temperature oxidation and predict the stress state of ceramic matrix and oxide layer during ablation process. Similarly, Wang et al. [13] have proposed a chemical force coupling model based on Parthasarathy’s oxidation kinetics model, comparing and analyzing the effect of different volume fractions of SiC on ablation resistance of ZrB2-SiC ceramics.

Though these studies can describe the changes in chemical reactions during the ablation process of UHTCs, the ultra-high-temperature ablation of ceramics is also accompanied by some discontinuous problems, such as thermal response, material damage, boundary movement and structural failure. Obviously, the previous studies did not consider these discontinuities. However, the peridynamics (PD) render unique advantages in dealing with discontinuities. The peridynamic theory was first proposed by Silling as a nonlocal theory [14], considering the advantages of meshless method and molecular dynamics, and avoiding the limitations of molecular dynamics on a computational scale. At the same time, material damage is a part of the PD constitutive relationship. So, the damage can develop spontaneously in the PD model [15]. With the continuous development of PD theory research, the coupled thermo-mechanical theory based on PD also appears. In 2014, Oterkus et al. [16] have deduced the fully-coupled peridynamic thermomechanics and given the bond-based PD equation and its dimensionless form. Subsequently, Liao et al. [17] have applied PD to the numerical simulations of transient heat transfer and thermo-mechanical coupling of functionally-graded materials. Wang et al. [18,19] proposed a weak coupled thermo-mechanical ordinary state-based peridynamic model to investigate the thermal cracking behavior of rocks under borehole heating and the crack propagation of ceramic plates in water quenching experiments respectively. The PD’s analysis results are consistent with the experimental observations. Chen et al. [20] have proposed a refined thermo-mechanical fully coupled peridynamics. This method did not require micro-conductivity and could directly represent the temperature deformation term by the integral of displacement difference. Based on the proposed model, a homogeneous concrete cylinder under thermal load was simulated from two-and three-dimensional angles, where the crack propagation results exhibited a good agreement with the experimental data. Crack propagation plays an important role in the failure analysis of brittle materials [21], and the crack initiation of ultra-high temperature ceramic materials is inevitable under extreme thermal load. Therefore, it is very important to characterize the crack initiation and development in numerical simulation. The bond-based PD method has been proved to be reliable in dealing with crack initiation and development [22]. Wang et al. [23,24] proposed an improved bond-based peridynamic coupled thermo-mechanical model based on multi-rate time integration. The model includes crack initiation and development, which can accurately simulate the crack initiation and development of brittle solids under thermal shock load, and accurately predict the radial and circumferential thermal crack growth of fuel rods. In recent years, there are many studies on the fracture mode of PD. Yu et al. [25] re-examined the relation between the critical bond stretch and the energy release rate for different crack types in Peridynamics by using Griffith’s “surface energy approach” and revealed that the choice of the critical bond stretch will be crucial for fracture model. Based on coupled thermo-mechanical peridynamic theory, Song et al. [26] established a coupled thermal mechanical periodic model with modified failure criterion suitable for de-icing technology. Li et al. [27] proposed an improved peridynamic stress expression, which can well predict the stress state at the crack tip.

The thermo-chemo-mechanical coupling needs to comprehensively consider the complex interactions between temperature field, chemical reaction and structural displacement. In recent years, based on the general thermodynamic principle and combined with the experimental data, scholars have studied the thermal force coupling of continuous solid materials [28,29]. Although these coupled thermo-chemo-mechanical models have rendered some good results in their respective research fields, these numerical theories are not suitable for dealing with discontinuities, such as damage and cracking. The ablation of UTHCs is a typical nonlinear and discontinuous problem. The PD theory can better deal with the ablation problem. PD solve various mechanical or thermal problems by using the integral equations based on non-local interaction to replace the differential equations of traditional continuum mechanics.

Based on single-phase ZrB2 ceramic, a novel thermo-chemo-mechanical coupling equation is established by combining the bond-based peridynamics (BB-PD) coupled thermo-mechanical with the high-temperature oxidation kinetics theory of single-phase ZrB2 ceramics in this paper. In Section 2, we will introduce the ablation mechanism, high-temperature oxidation kinetics theory of single-phase ZrB2 ceramics, and the bond-based peridynamic thermo-chemo-mechanical coupling model. Then, we present the initial boundary conditions of the model and damage criteria. Based on the coupled thermo-chemo-mechanical model (Section 2), the numerical simulation method is presented in Section 3. Furthermore, Section 4 shows that the proposed BB-PD coupled thermo-chemo-mechanical model is verified by numerical simulations and examples. Based on ZrB2 ablation, the calculation results of the fully coupled thermo-mechanical model and the coupled thermo-chemo-mechanical model are compared and analyzed. Finally, Section 5 presents the key conclusions of the current work.

2 The Coupled Thermo-Chemo-Mechanical Peridynamic Model

2.1 ZrB2 Ceramic Ablation Mechanism and Oxidation Kinetics Model

2.1.1 Ablation Mechanism of Single-Phase ZrB2 Ceramics

The relevant experimental results [10] reveal that, with the change of surface ablation temperature of ZrB2 ceramics, the following chemical reactions occur during the ablation process:

The abovementioned reaction process is described in Fig. 1. When the temperature exceeds 800 K, the oxidation of ZrB2 in the oxygen environment begins to become obvious and ZrB2 oxidizes to produce solid particles of ZrO2 and liquid B2O3. The former serves as a skeleton and the latter serves as a filling and cover. One should note that the oxygen needs to diffuse through the oxide film composed of ZrO2 and B2O3 in order to react with the underlying matrix material. The oxygen transport rate in the liquid film is very low, which drastically decreases the ablation rate. This is a typical inert oxidation phenomenon. The purpose of protecting the material is achieved by the slight oxidation of the material. When the temperature is higher than 1273 K, the evaporation of B2O3 on the surface becomes significant and, with the increase of temperature, B2O3 gradually shrinks into the pores of ZrO2. When the temperature is higher than 2073 K, the volatilization rate of B2O3 becomes greater than the generation rate and almost only porous ZrO2 remains in the oxide layer.

2.1.2 ZrB2 Ceramic Oxidation Kinetics Model

In this section, based on the oxidation process (Fig. 1), the oxidation kinetics model of single-phase ZrB2 ceramics is established (Fig. 2). This model can represent the ablation mechanism of ZrB2 in the intermediate temperature range (1273–2073 K) in a high-temperature aerobic environment.

As shown in Fig. 2, reaction (1) occurs at the matrix-oxidation interface (interface-1) and reaction (2) occurs at interface-2. The evolution of oxide layer thickness during oxidation can be calculated by the following equation [11]:

where

Similarly,

where R is the gas constant (8.314 J/(mol·K)). For q in Eq. (4), it is a constant that varies with temperature, and related to the liquid B2O3 layer thickness (

For the calculation of

where

2.2 Bond-Based Peridynamic Equation of Motion

The peridynamic theory was first proposed by Silling et al. [14] in 2000, which is commonly referred as the bond-based peridynamics (BB-PD) theory. In BB-PD, the motion equation of mechanics can be expressed as:

where

In the bond-based peridynamic theory, the force response function is the force density vector. Considering the effect of temperature change on deformation, the force density vector can be expressed as:

where

Combined with Eq. (10),

Herein,

Eq. (13) can be re-written as

Combined with Eqs. (9) and (13), the BB-PD motion equation considering temperature effect and chemical growth strain can be written as:

where c refers to the peridynamic parameter and the micro-elastic modulus. In the two-dimensional plane stress problem, the micro elastic modulus of isotropic material can be given as:

2.3 Bond-Based Peridynamic Equation of Heat Conduction

Based on Fourier’s law of heat transfer, the bond-based PD thermal conduction equation can be expressed as:

where cv refers to the specific heat capacity of the substance,

Introducing the temperature change term

where

its time derivative can be written as:

where

In summary, Eqs. (4), (15) and (19) form the bond-based peridynamic thermo-chemo-mechanical coupling equation for ultra-high temperature ablation of ZrB2 ceramics. The proposed model is a simple thermo-chemo-mechanical coupling. The chemical reaction only affects the deformation of the oxide layer, while the effect of the chemical reaction is not considered for the unreacted matrix. However, it should be noted that the model can perform fully coupled thermoelastic analysis. In this model, the structural deformation and temperature field affect each other. The detailed numerical realization method is described in Section 3.

2.4 Initial and Boundary Conditions

2.4.1 Thermal Boundary Conditions

In bond-based peridynamic thermal diffusion theory, an initial temperature is set for all material points. If the boundary temperature distribution

For the heat flow boundary, it is necessary to convert the heat flow into the heat source density and apply the heat source density to the material points of the boundary. Assuming that the heat flux at the boundary

If there is a temperature difference between the material surface (

where h refers to the convection heat transfer coefficient between boundary (

2.4.2 Mechanical Boundary Conditions

The boundary conditions in dynamics’ problem usually include displacement boundary, velocity boundary, acceleration boundary and force boundary. In peridynamic theory, displacement constraints, velocity constraints and acceleration constraints are represented by vector

When distributed pressure

2.4.3 Chemical Initial Conditions

It is assumed that the high-temperature oxidation of ZrB2 only takes place at the interface under the thermal load and only the initial conditions of chemical reaction need to be considered. In the oxidation kinetics model of Section 2.1, only the initial concentration of oxygen

In the solution of peridynamics, the displacement of each material point and elongation between each pair of material points are calculated. A time scalar function

where

where the material’s critical energy release rate

where

When

Combined with the Eqs. (15) and (19), the discrete form of the two-dimensional bond-based peridynamic coupled thermo-chemo-mechanical equation considering damage can be given as:

where the neighborhood volume of each material point

Eq. (32) shows that

The incremental form of the oxide layer thickness calculation using Eq. (4) can be given as:

It is important to note that the solution of the oxidation kinetics model belongs to a semi-independent process, which requires the distribution of temperature field. The bond-based PD coupled thermo-chemo-mechanical equation in this paper is based on the equations of motion and thermal conduction, and the kinetic equation of oxidation is used as an auxiliary. The calculation results of Eq. (33) need to be substituted into Eq. (32) to achieve the coupled solution of thermal conduction (Eq. (31)) and motion (Eq. (32)).

In addition, the material parameters of the newly generated oxide layer during the ablation process are quite different from the original substrate material. As shown in Fig. 4, we ignore the material loss and the oxide layer replaces the base material in the original position from the outside to the inside, and updates the material parameters of the oxide layer at each time step.

For the numerical solution of a classical fully-coupled thermo-mechanical equation, there are two generic methods, i.e., a single-step algorithm and an interleaved algorithm. In the single-step algorithm, also called the ensemble algorithm, time steps are simultaneously applied to the whole system of equations and multiple unknown variables are solved simultaneously. If the time integral of a single-step algorithm is implicit, absolute stability can usually be achieved, but it leads to a massive solution system. The interleaved algorithm partitions the coupled equations, solving the temperature and displacement fields with different time display algorithms. The interleaved algorithm can achieve a stable solution only under certain conditions.

Herein, according to the characteristics of the discrete form of the bond-based peridynamic equation and considering the particularity of the coupled thermo-chemo-mechanical model, an interleaved algorithm is used for the solution of coupled thermo-mechanical equations. The system is automatically divided according to the displacement field and temperature field, where the equation of motion is used for solving the displacement field and the equation of heat conduction is used for solving the temperature field. The solution of both equations is calculated by explicit time integration. Also, the computational equation of oxidation kinetics is relatively simple and can be solved directly.

The numerical calculations of the coupled thermo-chemo-mechanical model during the ablation of ultrahigh-temperature ceramics based on BB-PD mainly includes the following steps:

1. Defining the array, initializing variables and generating the physical numerical model;

2. Searching each material point within the neighborhood particles;

3. Initializing time-dependent function and surface correction factor;

4. Starting the first-time step cycle and calculating the temperature field;

5. Starting oxidation kinetics calculation and calculating the oxide layer film thickness produced due to the chemical reaction;

6. Starting the calculations of displacement field and reassigning material parameters to the material points of the oxide layer;

7. Judging whether the calculations are terminated and, if no, proceeding to Step 4 to start the next time step cycle;

8. Ending the cycle and obtaining the output results at each time step.

Fig. 5 shows the numerical calculations flow chart for solving the coupled thermo-chemo-mechanical model.

3.3 Numerical Stability Condition

The time integral of PD heat conduction equation is calculated by forward difference method, which is conditionally stable. In order to prevent the divergence of numerical solutions, it is necessary to give a stability condition to limit the time step of thermal diffusion. Referring to the existing literature, similar to the method used by Silling et al. [15], in order to meet the stability conditions of thermal diffusion, the critical time step

where,

4 Numerical Examples and Discussion

In this section, first, the correctness and accuracy of the bond-based peridynamic coupled thermo-chemo-mechanical numerical model are verified through three benchmark numerical examples: the fully coupled thermo-mechanical analysis of two-dimensional flat plate under temperature load; crack propagation simulations of the ceramic plate under cold shock, and evolution of oxide layer thickness of ZrB2 ceramic under high-temperature environments. In the fully coupled thermo-mechanical analysis of a two-dimensional flat plate under temperature load, the influence of neighborhood size and particle size on PD calculation results is investigated by comparing the calculation results under different neighborhood sizes and particle sizes with the analytical solution and the results of ABAQUS. The chemical reaction and coupled thermo-mechanical parts are verified separately because there is no corresponding ablation experimental research. Hence, the calculation results of the ZrB2 oxidation kinetics model can only be compared with the relevant high-temperature oxidation experiments of ZrB2. Accordingly, the coupled thermo-mechanical part can be verified by commercial software or an analytical method. Finally, considering the growth and evolution of damage and oxide layer, the coupled thermo-chemo-mechanical model and fully-coupled thermo-mechanical model are compared and analyzed.

4.1 Verification of the Fully Coupled Thermo-Mechanical and Convergence Study

As shown in Fig. 6, the side length of the square plate is 1 m, the initial temperature of the plate is 0°C, and the top side is experiencing a thermal load of 1°C. Then, the constraints in the x-direction are restricted at the left and right boundaries of the plate and the constraints in the y-direction are restricted at the lower boundary. Table 1 presents the parameters of the homogeneous plate. The Poisson’s ratio is 0.333. Points a and b are located on the vertical central axis of the plate, respectively, at the top and the middle, as shown in Fig. 6.

4.1.1 Influence of Neighborhood Size

In peridynamic theory, the material point only interacts with other material points in the neighborhood of this point, so the size of horizon will affect the numerical calculation results of PD. In order to investigate the influence of horizon size on numerical calculation, this section fixes the distance between material points as

Figs. 7 and 8 respectively show the displacement change diagram of point a in the vertical direction and the temperature change diagram of point b in 0–1 s, and the detailed enlarged diagrams are given on the right. As can be seen from Fig. 7, under the fixed particle spacing, when m is taken as 3, the displacement numerical calculation results of peridynamic are closer to those of analytical solution and ABAQUS. At the same time, it can be seen from the results of Fig. 8, when m is taken as 3, the temperature numerical calculation results of peridynamic are also closer to the analytical solution and ABAQUS calculation results. Therefore, when m = 3 the model can accurately simulate the coupled thermo-mechanical problem.

4.1.2 Influence of Particle Spacing

In addition to the influence of neighborhood range on the calculation results of bond-based peridynamic, the BB-PD model is also highly sensitive to particle spacing, i.e., the grid size. In this section, three different grid spacing will be selected in turn, i.e.,

In summary, when the horizon size is 3 times the particle spacing, the calculation results of the BB-PD model are accurate enough; when the particle spacing is smaller, the displacement and temperature calculation results are closer to the analytical solution and the calculation results of commercial software. The above results also prove that the bond-based peridynamic coupled thermo-chemo-mechanical numerical model given in this paper can be used to deal with the thermo-mechanical coupling response.

4.2 Verification of Damage Model

In this section, the crack growth process of the ceramic plate under cold shock is simulated and compared with the ceramic quenching experiment of Li et al. [34] to verify the BB-PD coupled thermo-chemo-mechanical model considering damage mentioned in this paper. In Li et al.’s experiment, the ceramic plate was heated to a certain temperature and rapidly dropped into the normal-temperature water. The high-temperature ceramic plate exhibited obvious shrinkage during quenching, forming several cracks on the surface. These parallel periodic cracks expand and drive inward. Referring to Li et al.’s experimental model [35], a quarter of the ceramic plate model is established. As shown in Fig. 11, the upper and right sides of the ceramic plate suffered cold shocks, and the displacement constraints in the x and y directions were applied on the left and lower sides of the ceramic plate, respectively.

The initial temperature of the ceramic plate is 773 K. The water temperature is 293 K and the coefficient of convection heat exchange is 70,000 W/(m2·K). Table 2 presents the material parameters of ceramic plates. In PD calculations, the distance between material points is

Fig. 12 presents the crack path of the ceramic plate under cold shock. The comparison of Figs. 12a–12c reveal that the PD-simulated crack path is basically same as that of Li et al.’s experimental [34] and numerical results [35], demonstrating the reliability of the BB-PD coupled thermo-chemo-mechanical model.

4.3 Verification of Oxidation Kinetics Model

In this section, the calculation model in Section 2.1 is used to simulate the ablation experiment of ZrB2 at different temperatures, and predict the growth and evolution of oxide layer thickness at different temperatures. The oxide layer thickness can be calculated by Eqs. (33) and (34). The proposed model results are compared with the experimental results of Opeka et al. [36]. In addition, in the calculation of ZrB2 oxidation kinetics, the density of ZrO2 ceramics is taken as

When ZrB2 is oxidized in air for 5 h, the porosity (f) of oxide layer is 0.05 and the pore diameter (d) is 0.5 μm. Fig. 13 shows the relationship between oxidation temperature and oxide layer thickness. It can be found that the oxidation temperature and oxide layer thickness followed the parabolic trend. The effect of oxidation temperature on the oxide layer thickness is quite significant. When the oxidation temperature is 1273 K, the thickness of oxide layer is 14.6 μm, which increases to 75.3 and 321.0 μm at the oxidation temperature of 1473 and 1673 K, respectively. It can be seen that the growth of oxide layer thickness is obviously accelerated with the increase of oxidation temperature. In addition, it can be found that the predicted results of ZrB2 oxidation kinetics model are consistent with the experimental results, confirming that the proposed oxidation kinetics model can be used to simulate the growth in oxide thickness of ZrB2 ceramics at high temperatures.

4.4 Comparative Analysis of Coupled Thermo-Chemo-Mechanical and Coupled Thermo-Mechanical Calculations

In this section, the bond-based PD coupled thermo-chemo-mechanical model is used to simulate the dynamic response of ZrB2 ceramic plates under thermal shock. The boundary conditions of the ceramic plate are the same as those of Section 4.2. As shown in Fig. 14, the left and bottom boundaries are constrained in the x or y directions, and the top and right boundaries are subjected to thermal shock loads. The dimensions of ZrB2 ceramic plate are

The distance between material points is

Figs. 15–17 present the temperature distribution, damage distribution diagrams of the ZrB2 ceramic plate, and compare with the temperature distribution results on the vertical centerline, obtained from both calculation methods. As shown in Fig. 15, the rate of temperature transfer of ZrB2 ceramic plate under thermo-chemo-mechanic coupling is lower than the thermo-mechanical coupling. According to Fig. 16, the maximum temperature of the ceramic plate with oxide layer is also lower than pure ZrB2 ceramic plate. It can be found that the ZrO2 oxide layer effectively resists the heat flow and inhibits temperature transfer on the ceramic plate to a certain extent, playing a certain protective role for the ZrB2 ceramic matrix. Fig. 17 shows the evolution of damage crack under two coupling calculations. When t = 2 ms, the model boundary damage and crack distribution are basically the same; When t = 4 ms, cracks begin to form in the ceramic plate without considering the evolution of oxide layer, but there is no crack initiation and development in the ceramic plate considering the evolution of oxide layer, which only shows that the local damage at the boundary develops to the inside; When t = 6 ms, the internal crack of the ceramic plate without considering the evolution of oxide layer continues to expand, while a crack appears in the upper left corner of the ceramic plate with considering the evolution of oxide layer, and the damage at the top boundary continues to expand inward; When t = 8 ms, the damage cracks of the two ceramic plates did not expand obviously. From the damage results of ceramic plates in Fig. 17, it can be concluded that the existence of oxide layer does enhance the thermal protection ability of ceramic plates and delay the damage caused by structural deformation caused by temperature change.

From the above analysis, it can be seen that the evolution of oxide layer caused by chemical reaction during the ceramic ablation can indeed affect the thermal protection performance and damage development of ZrB2 ceramics, and the BB-PD coupled thermo-chemo-mechanical model based on ZrB2 ceramic ablation can deal with the problems of damage, fracture and temperature effect during the ZrB2 ceramic ablation to a certain extent.

Based on the bond-based peridynamic coupled thermo-mechanical theory and high-temperature oxidation kinetics model, a simple coupled thermo-chemo-mechanical model is established to simulate the ultrahigh-temperature ablation of ceramics. Compared with the coupled thermo-mechanical model, the coupled thermo-chemo-mechanical model considers the oxide layer generated by chemical reaction during ablation and the influence of chemical growth strain caused by the oxide layer on the bond basis force response function of PD. Finally, the influence of oxide layer on the ablation protection of ultra-high-temperature ceramics is analyzed in detail. The key results can be summarized as:

1. The BB-PD coupled thermo-chemo-mechanical model based on the ablation of ultra-high-temperature ceramics is proposed, which can describe the heat conduction and damage process of ZrB2 ceramics under high-temperature ablation conditions. It is proved that this numerical model can be used to describe the high-temperature ablation process of ZrB2 ceramics, as verified by the cold impact damage experiment and evolution of oxidation layer thickness.

2. The BB-PD coupled thermo-chemo-mechanical model is also used to simulate the high-temperature ablation process of ceramics. One should note that the results are significantly different from the calculation results of the BB-PD coupled thermo-mechanical, indicating that the presence of a protective oxide layer plays a key role in protecting the internal materials and preventing heat transfer.

In general, the bond-based peridynamic coupled thermo-chemo-mechanical model can describe the high-temperature ablation process of ZrB2 ceramics to a certain extent, revealing the changes in temperature distribution and damage behavior during the ceramic ablation.

Funding Statement: The authors acknowledge the support from the National Natural Science Foundation of China (11972267).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

## References

- Zhang, X., Hu, P., Du, S., Han, J., & Meng, S. (2015). Research progress on ultra-high temperature ceramic composites.
*Chinese Science Bulletin*,*60(3)*, 257-266. [Google Scholar] - Kuriakose, A. K., & Margrave, J. L. (1963). Kinetics of the reactions of elemental fluorine with zirconium carbide and zirconium diboride at high temperatures.
*Journal of Physical Chemistry*,*68(2)*, 290-295. [Google Scholar] - Fahrenholtz, W. G., Hilmas, G. E., Talmy, I. G., & Zaykoski, J. A. (2007). Refractory diborides of zirconium and hafnium.
*Journal of the American Ceramic Society*,*90(5)*, 1347-1364. [Google Scholar] - Han, J., Hu, P., Zhang, X., Meng, S., & Han, W. (2007). Oxidation-resistant ZrB-SiC composites at 2200°C.
*Composites Science and Technology*,*68(3)*, 799-806. [Google Scholar] - Wang, Y., Zhu, X., Zhang, L., & Cheng, L. (2011). Reaction kinetics and ablation properties of C/C-ZrC composites fabricated by reactive melt infiltration.
*Ceramics International*,*37(4)*, 1277-1283. [Google Scholar] - Qing, X., Sun, W., Tian, T., Xiong, X., & Zhang, H. (2020). Structural characteristics and ablative behavior of C/C-ZrC-SiC composites reinforced with Z-pins like Zr-Si-B-C multiphase ceramic rods.
*Ceramics International*,*46(11)*, 18895-18902. [Google Scholar] - Fahrenholtz, W. G. (2005). The ZrB volatility diagram.
*Journal of the American Ceramic Society*,*88(12)*, 3509-3512. [Google Scholar] - Talmy, I., Zaykoski, J., Opeka, M. (2008). Properties of ceramics in the ZrB2/ZrC/SiC system prepared by reactive processing. In: Ceramic engineering and science proceedings, vol. 19, no. 3, pp. 105–112. USA.
- Fahrenholtz, W. G. (2007). Thermodynamic analysis of ZrB-SiC oxidation: Formation of a SiC-depleted region.
*Journal of the American Ceramic Society*,*90(1)*, 143-148. [Google Scholar] - Parthasarathy, T., Rapp, R. A., Opeka, M., & Kerans, R. J. (2008). A model for transitions in oxidation regimes of ZrB.
*Materials Science Forum*,*595(1)*, 823-832. [Google Scholar] - Parthasarathy, T., Rapp, R., Opeka, M., & Kerans, R. (2007). A model for the oxidation of ZrB, HfB and TiB.
*Acta Materialia*,*55(17)*, 5999-6010. [Google Scholar] - Zhou, Z., Peng, X., & Wei, Z. (2014). A thermo-chemo-mechanical model for the oxidation of zirconium diboride.
*Journal of the American Ceramic Society*,*98(2)*, 629-636. [Google Scholar] - Wang, H., & Shen, S. (2017). A chemomechanical coupling model for oxidation and stress evolution in ZrB-SiC.
*Journal of Materials Research*,*32(7)*, 1-12. [Google Scholar] - Silling, S. (2000). Reformulation of elasticity theory for discontinuities and long-range forces.
*Journal of the Mechanics and Physics of Solids*,*48(1)*, 175-209. [Google Scholar] - Silling, S., & Askari, E. (2005). A meshfree method based on the peridynamic model of solid mechanics.
*Computers and Structures*,*83(17)*, 1526-1535. [Google Scholar] - Oterkus, S., Madenci, E., & Agwai, A. (2014). Fully coupled peridynamic thermomechanics.
*Journal of the Mechanics and Physics of Solids*,*64(1)*, 1-23. [Google Scholar] - Liao, Y., Liu, L., Liu, Q., Lai, X., & Assefa, M. (2017). Peridynamic simulation of transient heat conduction problems in functionally gradient materials with cracks.
*Journal of Thermal Stresses*,*40(12)*, 1-18. [Google Scholar] [CrossRef] - Wang, Y., & Zhou, X. (2019). Peridynamic simulation of thermal failure behaviors in rocks subjected to heating from boreholes.
*International Journal of Rock Mechanics and Mining Sciences*,*117(1–2)*, 31-48. [Google Scholar] [CrossRef] - Wang, Y., Zhou, X., & Zhang, T. (2019). Size effect of thermal shock crack patterns in ceramics: Insights from a nonlocal numerical approach.
*Mechanics of Materials*,*137(8)*, 103133. [Google Scholar] [CrossRef] - Chen, W., Gu, X., Zhang, Q., & Xia, X. (2021). A refined thermo-mechanical fully coupled peridynamics with application to concrete cracking.
*Engineering Fracture Mechanics*,*242(8)*, 107463. [Google Scholar] [CrossRef] - Cheng, Z., Wang, Z., & Luo, Z. (2019). Dynamic fracture analysis for shale material by peridynamic modelling.
*Computer Modeling in Engineering & Sciences*,*118(3)*, 509-527. [Google Scholar] [CrossRef] - Chu, B., Liu, Q., Liu, L., Lai, X., & Mei, H. (2020). A rate-dependent peridynamic model for the dynamic behavior of ceramic materials.
*Computer Modeling in Engineering & Sciences*,*124(1)*, 151-178. [Google Scholar] [CrossRef] - Wang, Y., Zhou, X., & Kou, M. (2019). An improved coupled thermo-mechanic bond-based peridynamic model for cracking behaviors in brittle solids subjected to thermal shocks.
*European Journal of Mechanics-A/solids*,*73(3)*, 282-305. [Google Scholar] [CrossRef] - Wang, Y., Zhou, X., & Kou, M. (2018). Peridynamic investigation on thermal fracturing behavior of ceramic nuclear fuel pellets under power cycles.
*Ceramics International*,*44(10)*, 11512-11542. [Google Scholar] [CrossRef] - Yu, H., & Li, S. (2020). On energy release rates in peridynamics.
*Journal of the Mechanics and Physics of Solids*,*142(3)*, 104024. [Google Scholar] [CrossRef] - Song, Y., Li, S., & Zhang, S. (2021). Peridynamic modeling and simulation of thermo-mechanical de-icing process with modified ice failure criterion.
*Defence Technology*,*17(1)*, 15-35. [Google Scholar] [CrossRef] - Li, J., Li, S., Lai, X., & Liu, L. (2022). Peridynamic stress is the static first piola-kirchhoff virial stress.
*International Journal of Solids and Structures*,*241(1–2)*, 111478. [Google Scholar] [CrossRef] - Loeffel, K., & Anand, L. (2011). A chemo-thermo-mechanically coupled theory for elastic-viscoplastic deformation, diffusion, and volumetric swelling due to a chemical reaction.
*International Journal of Plasticity*,*27(9)*, 1409-1431. [Google Scholar] [CrossRef] - Xu, G., Yang, L., Zhou, Y., Pi, Z., & Zhu, W. (2019). A chemo-thermo-mechanically constitutive theory for thermal barrier coatings under CMAS infiltration and corrosion.
*Journal of the Mechanics and Physics of Solids*,*133(1)*, 1-19. [Google Scholar] [CrossRef] - Clarke, D. (2003). The lateral growth strain accompanying the formation of a thermally grown oxide.
*Acta Materialia*,*51(5)*, 1393-1408. [Google Scholar] [CrossRef] - Lee, J., Oh, S. E., & Hong, J. (2017). Parallel programming of a peridynamics code coupled with finite element method.
*International Journal of Fracture*,*203(1)*, 99-114. [Google Scholar] [CrossRef] - Timoshenko, S. P., Goodier, J. N. (1970). Theory of elasticity. New York: Mcgraw-Hill.
- Carslaw, H. S., Jaeger, J. C. (1959). Conduction of heat in solids. Oxford: Clarendon Press.
- Li, J., Song, F., & Jiang, C. (2013). Direct numerical simulations on crack formation in ceramic materials under.
*Journal of the European Ceramic Society*,*33(13)*, 2677-2687. [Google Scholar] [CrossRef] - Li, J., Song, F., & Jiang, C. (2015). A non-local approach to crack process modeling in ceramic materials subjected to thermal shock.
*Engineering Fracture Mechanics*,*133(4)*, 85-98. [Google Scholar] [CrossRef] - Opeka, M. M., Talmy, I. G., Wuchina, E. J., Zaykoski, J. A., & Causey, S. J. (1999). Mechanical, thermal, and oxidation properties of refractory hafnium and zirconium compounds.
*Journal of the European Ceramic Society*,*19(13)*, 2405-2414. [Google Scholar] [CrossRef] - Justin, J. F., & Jankowiak, A. (2011). Ultra high temperature ceramics: Densification, properties and thermal stability.
*High Temperature Materials*,*1(3)*, 1-12. [Google Scholar] - Tan, Y. (2020). Peridynamic simulation of functionally graded materials under thermal shock loading (Master Thesis). Wuhan, China: Wuhan University of Technology.

## Cite This Article

Li, Y., Liu, Q., Liu, L., Mei, H. (2023). The Coupled Thermo-Chemo-Mechanical Peridynamics for ZrB_{2}Ceramics Ablation Behavior.

*CMES-Computer Modeling in Engineering & Sciences, 135(1)*, 417–439.