The Coupled Thermo-Chemo-Mechanical Peridynamics for ZrB2 Ceramics Ablation Behavior

Yuanzhe Li,Qiwen Liu,Lisheng Liu and Hai Mei

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

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

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 ZrB2 ceramics oxidation kinetics model and coupled thermomechanical 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 ZrB2 ceramics;ablation;coupled thermo-chemo-mechanic;peridynamics model;oxide layer

1 Introduction

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 ZrO2layer 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 ZrO2layer when ZrC was exposed to a high-temperature oxygen environment.When the temperature is increased to 2800°C,the solid ZrO2transforms into a molten state and covers the material surface.The formation of a dense ZrO2layer 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 B2O3film in the low-temperature zone,SiO2film in the medium-temperature zone and ZrO2film 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 ZrB2ceramic 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, ZrO2and B2O3.The proposed model can also predict the mass fraction of liquid B2O3remaining 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, HfB2and 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 ZrB2and proposed a thermal force coupling model, which could describe the multi-field coupling behavior of ZrB2during 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 thermomechanical 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 twoand 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 ZrB2ceramic, 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 ZrB2ceramics in this paper.In Section 2, we will introduce the ablation mechanism, high-temperature oxidation kinetics theory of single-phase ZrB2ceramics, 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 thermochemo-mechanical model is verified by numerical simulations and examples.Based on ZrB2ablation,the calculation results of the fully coupled thermo-mechanical model and the coupled thermo-chemomechanical 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 ZrB2ceramics,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 ZrB2in the oxygen environment begins to become obvious and ZrB2oxidizes to produce solid particles of ZrO2and 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 ZrO2and B2O3in 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 B2O3on the surface becomes significant and,with the increase of temperature,B2O3gradually shrinks into the pores of ZrO2.When the temperature is higher than 2073 K,the volatilization rate of B2O3becomes greater than the generation rate and almost only porous ZrO2remains in the oxide layer.

Figure 1:Schematic diagram of ZrB2 oxidation process with respect to temperature

2.1.2 ZrB2 Ceramic Oxidation Kinetics Model

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

Figure 2:Schematic diagram of the oxidation model of ZrB2

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]:

whereRis the gas constant (8.314 J/(mol·K)).Forqin Eq.(4), it is a constant that varies with temperature, and related to the liquid B2O3layer thickness (hint) and the oxide layer thickness (L),which can be defined as:

For the calculation ofhint,DO2andDB2O3please refer to Parthasarathy et al.[11],which will not be described in detail here.The growth strain appears with the growth of oxide layer.Based on the phenomenon of oxide layer growth,Clarke[30]has proposed the following equation to describe the growth strain and oxide layer thickness growth rate,which be written as:

whereεgrepresents the chemical growth strain of the oxide layer;Doxdenotes the chemical reaction growth coefficient;hox(t)corresponds to the thickness of the oxide layer,which isLin Eq.(4).

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ρ(x)represents the mass density of the material pointx.(x,t)represents the second derivative of the displacementu(x,t)of the material pointxat timet.f (u′－u,x′－x,t)is the force response function,which is defined as the force loss per unit volume square applied by the material pointx′on the material pointx,andb(x,t)represents the physical density of the material pointxat timet.As shown in Fig.3,yandy′refer to the position vectors of the material pointsxandx′after displacementsuandu′,respectively.

Figure 3:The interactions of material points x and x′before and after deformation

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:

whereT(u′－u,x′－x,t)represents the average value of the temperature change of material pointxand

Combined with Eq.(10),s(u′－u,x′－x,t)represents the distance between two material points andαT(u′－u,x′－x,t)refers to the thermal expansion deformation caused by temperature change.In classical continuum mechanics,the elastic strain can be expressed as:

Herein,εerepresents elastic strain.εtotrefers to total strain andεTcorresponds to thermal strain.It can be found thatεtot－εTcorresponds tos－αTavgin Eq.(10).Combined with the oxidation kinetics in Section 2.1,the total strain of the oxide layer can be expressed as:

Eq.(13) can be re-written asεe=εtot－εT－εgto replaces－αTavgof Eq.(10), which can be equivalent that the elongation of the bond minus the thermal expansion deformation and the chemical reaction deformation.At this time,the bond force response function can be written as:

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

wherecrefers 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:whereErepresents the elastic modulus.The definitions of initial relative position vector and relative displacement vector areξ=x′－xandη=u′－u,respectively.αrepresents the coefficient of linear thermal expansion andTavgrefers to the average temperature change between material pointsxandx′,which can be defined 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:

wherecvrefers to the specific heat capacity of the substance,T(x,t)represents the temperature of the material pointxat timet,hs(x,t)represents the heat source density,andfhcorresponds to the thermal response function,which is solely governed by the interactions between material pointsxandx′and can be expressed as:

Introducing the temperature change termdue to deformation and substituting Eq.(18)into Eq.(17),the bond-based PD thermal conduction equation containing the deformation effect term on temperature can be written as:

its time derivative can be written as:

whereκrepresents the micro-thermal conductivity of PD,which is related to the thermal conductivity(kT) of the classical thermodynamic theory and micro-thermal conductivity isκ= 6kT/(πhδ3)in two-dimensional problems.

In summary,Eqs.(4),(15)and(19)form the bond-based peridynamic thermo-chemo-mechanical coupling equation for ultra-high temperature ablation of ZrB2ceramics.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 distributionTBC(x,t)is known,the boundary temperature can be applied directly to the material points of the boundary layerRBC,which can be given as:

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 boundaryRBCisqBC,the heat flux at the boundaryRBCcan be given as:

If there is a temperature difference between the material surface(RBC)and surrounding environment, the temperature difference will cause the heat transfer between the material and environment medium.Hence, the temperature boundary belongs to the heat convection boundary, which can be given as:

wherehrefers to the convection heat transfer coefficient between boundary (RBC) and surrounding environment.Δis the thickness of the boundaryRBC.

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 vectorUBC, vectorVBCand vectoraBC, respectively.These constraints can be applied directly to the particles of the outer boundary(Rs):

When distributed pressurep(x,t)or concentrated forceP(t)are applied on the boundary(RBC)layer surface,the body force density vectors can be expressed as:

2.4.3 Chemical Initial Conditions

It is assumed that the high-temperature oxidation of ZrB2only 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 oxygenin the chemical reaction needs to be considered.The initial concentration of oxygenremains constant throughout the chemical reaction.

2.5 Damage Model

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μ(x,x′,t)[15] is introduced to deal with the damage caused by temperature changes.The criteria for this time scalar function can be given as:

wheres0refers to the critical elongation of the bond in the planar stress state and can be expressed as:

where the material’s critical energy release rateG0is related to the fracture toughnessKIC, and their relationship with the elastic modulusEcan be expressed as:G0=.Whenμij=1,it indicates that the bond between material pointsxandx′is intact.Whenμij= 0,it indicates that the bond between material pointsxandx′is broken.Herein, the temperature bond is not adherent to the force bond and,when the force bond is broken,the heat transfer efficiency of the temperature bond is weakened.The corresponding numerical expression can be given as:

whereφi(t)represents the local damage to the material pointxi,which is related to the time-dependent functionμijof the bond and can be defined as:

Whenφ= 1, it indicates that all bonds within the neighborhood of the substance points are broken.Whenφ= 0,it indicates that the bonds within the neighborhood of the material points are intact.

3 Numerical Implementation

3.1 Numerical Discretization

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 pointxiisViand the neighborhood containsNiinteracting material pointxj.

Eq.(32) shows thathoxneeds to be solved for solving this coupled formula further, while the thickness of the oxide layer att(n+1)can be obtained from the calculations,as given below:

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

Figure 4:The evolution of oxide layer during ablation

3.2 Solution Procedue

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-chemomechanical model.

Figure 5:The flow diagram of numerical simulations

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 stepof temperature field calculation can be calculated by the following formula:

where,is the critical thermal diffusion time step.In order to meet the stability conditions of numerical calculation of thermal diffusion, the thermal diffusion time stepΔtTHshould meet the following conditions:

ΔtTHis the time step of thermal diffusion calculation.After each thermal diffusion time step,the Adaptive Dynamic Relaxation (ADR) method [31] is used by introducing artificial damping to the solution system, and the static or quasi-static solution in the dynamic problem is solved by the explicit forward difference time integration method.The static or quasi-static solution is the steadystate part of the instantaneous response solution.When the static or quasi-static solution is obtained,the displacement and damage caused by the coupled thermo-mechanical material points will be solved as the initial parameters of the next thermal diffusion time step.In the ADR method,the time step of displacement calculation is=1 s.

4 Numerical Examples and Discussion

In this section,first,the correctness and accuracy of the bond-based peridynamic coupled thermochemo-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 ZrB2ceramic 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 ZrB2oxidation 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 thex-direction are restricted at the left and right boundaries of the plate and the constraints in they-direction are restricted at the lower boundary.Table 1 presents the parameters of the homogeneous plate.The Poisson’s ratio is 0.333.Pointsaandbare located on the vertical central axis of the plate,respectively,at the top and the middle,as shown in Fig.6.

Figure 6:The top boundary of the plate is under thermal load

Table 1: The parameters of the homogeneous plate

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Δx= 0.005 m and horizon size is set asδ=mΔx.The time step isΔt= 0.00001 s.The plate model is divided into 200×200 = 40000 particles,select different horizon sizeδ= 2Δx,δ= 3Δx,δ= 4Δx,respectively,i.e.,m=2,3,4.We compare the calculation results of PD with the results of ABAQUS and analytical method,and these three methods are used to calculate the vertical displacement of pointaand temperature of pointb, respectively, where the analytical solution can be obtained from Timoshenko et al.[32]and Carslaw et al.[33],as shown in Eqs.(37)and(38).

Figs.7 and 8 respectively show the displacement change diagram of pointain the vertical direction and the temperature change diagram of pointbin 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, whenmis 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,whenmis taken as 3,the temperature numerical calculation results of peridynamic are also closer to the analytical solution and ABAQUS calculation results.Therefore,whenm=3 the model can accurately simulate the coupled thermo-mechanical problem.

Figure 7:The vertical displacement of point a,(b)is an enlarged detail from(a)

Figure 8:The temperature change of point b,(b)is an enlarged detail from(a)

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., Δx= 0.005 m, Δx= 0.01 m,Δx= 0.02 m,and the neighborhood size is taken as 3 times of particle spacing in Section 4.1.1,i.e.,δ=3Δx.Figs.9 and 10 show the displacement change diagram of pointain the vertical direction and the temperature change diagram of pointb.It can be seen that when the neighborhood range is fixed asm= 3, the denser the grid distribution, the closer the displacement and temperature calculation results of PD are to the analytical solution and ABAQUS’s results.

Figure 9:The vertical displacement of point a,(b)is an enlarged detail from(a)

Figure 10:The temperature change of point b,(b)is an enlarged detail from(a)

In summary,when the horizon size is 3 times the particle spacing,the calculation results of the BBPD 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-chemomechanical 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 thermochemo-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 thexandydirections were applied on the left and lower sides of the ceramic plate,respectively.

Figure 11:Geometry and boundary conditions of ceramic plates

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Δ= 0.05 mm,the time step isΔt=0.0001 s and the neighborhood size isδ=3Δ.

Table 2: The parameters of ceramic plates

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.

Figure 12:The comparison of the final crack path

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 ZrB2at 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 ZrB2oxidation kinetics,the density of ZrO2ceramics is taken asρ= 5600 kg/m3.The ambient oxygen concentrationis set as 9.69 mol/m3.

When ZrB2is 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 ZrB2oxidation 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 ZrB2ceramics at high temperatures.

Figure 13:The comparison between theoretically calculated and experimentally measured oxide layer thickness

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 ZrB2ceramic 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 thexorydirections,and the top and right boundaries are subjected to thermal shock loads.The dimensions of ZrB2ceramic plate are 0.015 m×0.005 m.The initial temperature is 293 K and the ambient temperature is 2293 K.The convection heat transfer coefficient is 70,000 W/(m2·K).

Figure 14:The boundary conditions of ZrB2 ceramic plate

The distance between material points is Δ=0.00005 m,the horizon size isδ=3Δ,and the time step isΔt=0.00001 s.Table 3 presents the parameters of ZrB2matrix and ZrO2oxide layer.In order to analyze the influence of chemical reactions on the ablation process,the dynamic response of ZrB2ceramic plate under fully coupled thermo-mechanical conditions is simulated.

Table 3: Thermo-mechanical parameters of ZrB2 matrix and ZrO2 oxide layer

Table 3 (continued)Parameters ZrB2[37] ZrO2[38]Thermal conductivity coefficient kT (W/(m·K)) 60 2.063 Thermal expansion coefficient α(1/K) 5.9×10－6 7.11×10－6 Fracture toughness KIC (MPa·m1/2) 4.52×106 5×106

Figs.15–17 present the temperature distribution, damage distribution diagrams of the ZrB2ceramic 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 ZrB2ceramic 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 ZrB2ceramic plate.It can be found that the ZrO2oxide 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 ZrB2ceramic 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.

Figure 15: The comparison of temperature (K) nephogram under (a) Thermo-mechanical coupling calculations and(b)Thermo-chemo-mechanical coupling calculations

Figure 16:The comparison of temperature distribution on the vertical centerline of ZrB2 ceramic plate under different calculation methods

Figure 17:The comparison of damage nephogram under(a)Thermo-mechanical coupling calculations and(b)Thermo-chemo-mechanical coupling calculations

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 ZrB2ceramics, and the BB-PD coupled thermo-chemo-mechanical model based on ZrB2ceramic ablation can deal with the problems of damage,fracture and temperature effect during the ZrB2ceramic ablation to a certain extent.

5 Conclusions

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-hightemperature ceramics is proposed,which can describe the heat conduction and damage process of ZrB2ceramics under high-temperature ablation conditions.It is proved that this numerical model can be used to describe the high-temperature ablation process of ZrB2ceramics, 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 hightemperature 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 ZrB2ceramics 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.