Thermo-Mechanical Coupling Model of Bond-Based Peridynamics for Quasi-Brittle Materials

The mechanical properties of quasi-brittle materials, which are widely used in engineering applications, are often affected by the thermal condition of their service environment. Moreover, the materials appear brittle when subjected to tensile loading and show plastic characteristics under high pressure. These two phenomena manifest under different circumstances as completely different mechanical behaviors in the material. To accurately describe the mechanical response, the material behavior, and the failure mechanism of quasi-brittle materials with the thermo-mechanical coupling effect, the influence of the thermal condition is considered in calculating bond forces in the stretching and compression stages, based on a new bond-based Peridynamic (BB-PD) model. In this study, a novel bond-based Peridynamic, fully coupled, thermo-mechanical model is proposed for quasi-brittle materials, with a heat conduction component to account for the effect of the thermo-mechanical coupling. Numerical simulations are carried out to demonstrate the validity and capability of the proposed model. The results reveal that agreement could be found between our model and the experimental data, which show good reliability and promise in the proposed approach.


Introduction
Thanks to their unique brittle characteristics, materials such as ceramics, rocks, and concrete are widely used in various engineering fields, including the aerospace, armor protection, and construction industries. However, temperature, as an essential environmental factor, can impact the mechanical properties and service performance of the above materials during operation [1]. It may even cause serious accidents associated with the thermal cracking of surrounding rocks upon the underground storage of nuclear waste [2] and heat-induced spalling of building concrete [3]. Therefore, to ensure the safety and load-bearing capacity of quasi-brittle engineering materials, their mechanical properties and damage behavior must be understood under the thermos-mechanical coupling.
To study the behavior of quasi-brittle materials undergoing thermal damage, numerous experimental facilities have been developed and utilized, such as X-ray computed tomography (X-CT) [4], acoustic emission (AE) [5], scanning electron microscopy (SEM) [6], and thermal stress devices (TSD) [7]. These approaches have enabled researchers to significantly improve their knowledge about the hermos-mechanical properties of the materials. However, experimental studies are not only time-consuming and laborious, but it is also challenging to observe the sprouting and extension of micro-cracks inside the material in situ. Thus, by providing a detailed and cost-effective prediction, numerical methods shed insight into the mechanical failure processes in quasi-brittle materials.
Extensive numerical investigations have been devoted to understanding the thermal damage behavior of quasi-brittle materials, which are generally based on continuum to cold, uniaxial compression experiments after heat treatment were conducted, and the coincidence between the simulated and experimental results were analyzed.

Thermo-Mechanical Coupling Model
In this section, the classic fully coupled thermo-mechanical BB-PD model, in which the bonds remain elastically deformed during deformation until broken and are not applicable to quasi-brittle materials, is first introduced. Then, the mechanical behavior of quasi-brittle materials in tension and compression is presented. The Peridynamics model for quasi-brittle materials is presented in detail in Section 2.3, and the substance of this paper, i.e., the study of thermal effects in the tensile and compression phases, is presented. Finally, the numerical discretization and time integration methods of the proposed model are presented.

Fully Coupled Thermo-Mechanical Equation
Unlike the partial derivative of deformation with respect to the spatial coordinate in continuum mechanics, the BB-PD theory adopts the spatial integral equation, which can be applied to discontinuous bodies [21]. As shown in Figure 1, each material point x in a region R of an object interacts with all other points within its neighborhood radius H x . For the domain x , the material point is called a neighborhood particle of the material point x. When a solid is deformed under an external load, the matter points x and x arrive at the post-deformation positions y and y through displacements u and u , respectively, and |ξ| is the distance between the two particles before the deformation and |ξ + η| is the distance between the two particles after the deformation.
In this study, a model suitable for the failure of quasi-brittle materials under thermomechanical coupling is proposed, and the effect of temperature within the elastic and plastic stages is considered. To verify the reliability and validity of this model, numerical simulations of ceramics under a heating load and the pre-cracked Brazilian disk under uniaxial compression were conducted. Finally, two-dimensional granite plates were exposed to cold, uniaxial compression experiments after heat treatment were conducted, and the coincidence between the simulated and experimental results were analyzed.

Thermo-Mechanical Coupling Model
In this section, the classic fully coupled thermo-mechanical BB-PD model, in which the bonds remain elastically deformed during deformation until broken and are not applicable to quasi-brittle materials, is first introduced. Then, the mechanical behavior of quasi-brittle materials in tension and compression is presented. The Peridynamics model for quasi-brittle materials is presented in detail in Section 2.3, and the substance of this paper, i.e., the study of thermal effects in the tensile and compression phases, is presented. Finally, the numerical discretization and time integration methods of the proposed model are presented.

Fully Coupled Thermo-Mechanical Equation
Unlike the partial derivative of deformation with respect to the spatial coordinate in continuum mechanics, the BB-PD theory adopts the spatial integral equation, which can be applied to discontinuous bodies [21]. As shown in Figure 1, each material point x in a region R of an object interacts with all other points within its neighborhood radius x H . For the domain x′ , the material point is called a neighborhood particle of the material point x . When a solid is deformed under an external load, the matter points x and x′ arrive at the post-deformation positions y and y′ through displacements u and u′ , respectively, and ξ is the distance between the two particles before the deformation and ξ η + is the distance between the two particles after the deformation. In contrast to the classical derivation of the heat equation [22], PD laws are derived based on irreversible thermo-mechanicals, i.e., energy conservation and free energy density functions. The fully coupled thermo-mechanical equations for BB-PD are as follows [23][24][25]: In contrast to the classical derivation of the heat equation [22], PD laws are derived based on irreversible thermo-mechanicals, i.e., energy conservation and free energy density functions. The fully coupled thermo-mechanical equations for BB-PD are as follows [23][24][25]: Here, Equation (1) is the PD thermal diffusion equation with a structural coupling term, where c v is the specific heat capacity, κ is the thermal conductivity of the bond in the PD system, h s (x, t) is the rate of heat production per unit volume at the point of matter x at time t. The temperature difference between the substance points x and x can be expressed as τ = T(x, t) − T(x , t), where T 0 (cα/2) . e are the deformation terms caused by heating and cooling, T 0 is the reference temperature, and . e is the rate of change of bond lengths, denoted as: Equation (2) is the equation of motion for the BB-PD system with a temperature coupling term. In this equation, ρ is the mass density; ..
u(x, t) is the acceleration of the matter point x at time t; H is the neighborhood range of the matter point x; and V is the volume of the matter point x ; c is the bond constant; α is the coefficient of thermal expansion; b(x, t) is the force density of the matter point x at time t; s is the stretch rate between x and x ; and T avg is the average temperature of the matter points x and x , respectively, expressed as:

The Characterization of the Mechanical Behavior of Quasi-Static Brittle Materials
For quasi-brittle materials such as ceramics and concrete, the damage models under compression and tension are distinct. Under tension, the behavior is brittle, while a more ductile behavior can be observed under compression. The brittle behavior of materials under tensile loading is attributed to macro-crack formation. On the other hand, the ductile behavior of ceramics under compression can be explained by micro-crack formation and plasticity [26][27][28].
Under tensile loading, the quasi-brittle materials undergo direct, brittle damage at the end of elastic deformation, with essentially no plastic deformation. Figure 2 depicts a typical stress-strain diagram of the quasi-brittle material model in the tensile phase, showing a linear relationship between stress and strain, where the material breaks down and loses its load-bearing capacity after reaching the tensile-strength limit of the material. x′ can be expressed as , , where 0 ( /2) T c e α  are the deformation terms caused by heating and cooling, 0 T is the reference temperature, and e  is the rate of change of bond lengths, denoted as: e ξ η η ξ η Equation (2)  x and x′ , respectively, expressed as:

The Characterization of the Mechanical Behavior of Quasi-Static Brittle Materials
For quasi-brittle materials such as ceramics and concrete, the damage models under compression and tension are distinct. Under tension, the behavior is brittle, while a more ductile behavior can be observed under compression. The brittle behavior of materials under tensile loading is attributed to macro-crack formation. On the other hand, the ductile behavior of ceramics under compression can be explained by micro-crack formation and plasticity [26][27][28].
Under tensile loading, the quasi-brittle materials undergo direct, brittle damage at the end of elastic deformation, with essentially no plastic deformation. Figure 2 depicts a typical stress-strain diagram of the quasi-brittle material model in the tensile phase, showing a linear relationship between stress and strain, where the material breaks down and loses its load-bearing capacity after reaching the tensile-strength limit of the material.  When a quasi-brittle material is subjected to a compressive load, the model remains intact at the initial stage, and as the load continues to increase, structural crushing occurs internally due to the generation of micro-cracks and lattice plasticity, and the interaction of these cracks leads to a decrease in the compressive strength of the quasi-brittle material, a phenomenon described as plastic-softening behavior. The quasi-brittle material is treated as elastic material before damage occurs, and after damage occurs, it is treated as material that remains intact but whose strength decreases with the accumulation of damage. Figure 3 shows a typical stress-strain diagram for the ideal quasi-brittle material which is typical of ultra-high-strength concrete [29]. Before reaching the elastic limit of the material, for a sufficiently small segment, the relation is close to linear. After exceeding the elastic limit, micro-cracks appear in the concrete, resulting in an inelastic response that differs from the plastic-flow behavior of ductile materials, as shown by the fact that the strength of the material decreases with the accumulation of plastic damage, i.e., the Materials 2022, 15, 7401 5 of 20 material begins to soften. Finally, the concrete is completely broken and loses its loadbearing capacity. age. Figure 3 shows a typical stress-strain diagram for the ideal quasi-brittle material which is typical of ultra-high-strength concrete [29]. Before reaching the elastic limit of the material, for a sufficiently small segment, the relation is close to linear. After exceeding the elastic limit, micro-cracks appear in the concrete, resulting in an inelastic response that differs from the plastic-flow behavior of ductile materials, as shown by the fact that the strength of the material decreases with the accumulation of plastic damage, i.e., the material begins to soften. Finally, the concrete is completely broken and loses its load-bearing capacity.

Description of the Stretching Stage
The mechanical behavior of quasi-brittle materials in the tensile phase, as described in Section 2.2, has been known to exhibit mainly brittle characteristics. When using PD to describe the behavior of brittle materials in the tensile phase, the bond forces are considered to be related only to the relative elongation of the bonds and can be described as: where c is the micro-elastic modulus, obtained from the consistency between the strain energy density of the isotropic material and the theoretical strain energy density of the continuum mechanics. The micro-elastic modulus of the isotropic material in the plane stress state is ( ) , where E is the elastic modulus and δ is the neighborhood radius.

Description of the Compression Stage
In the compression stage, the quasi-brittle material can be divided into two phases: the elastic phase before reaching the elastic limit and the plastic phase after exceeding the elastic limit. For the elastic stage, the bond force is calculated in the same way as for the tensile stage. For the plastic stage, the quasi-brittle material exhibits plastic softening, and the influence of plastic deformation must be considered in the calculation of the bond force, which can be expressed as:  The mechanical behavior of quasi-brittle materials in the tensile phase, as described in Section 2.2, has been known to exhibit mainly brittle characteristics. When using PD to describe the behavior of brittle materials in the tensile phase, the bond forces are considered to be related only to the relative elongation of the bonds and can be described as: where c is the micro-elastic modulus, obtained from the consistency between the strain energy density of the isotropic material and the theoretical strain energy density of the continuum mechanics. The micro-elastic modulus of the isotropic material in the plane stress state is c = 9E/ πhδ 3 , where E is the elastic modulus and δ is the neighborhood radius.

Description of the Compression Stage
In the compression stage, the quasi-brittle material can be divided into two phases: the elastic phase before reaching the elastic limit and the plastic phase after exceeding the elastic limit. For the elastic stage, the bond force is calculated in the same way as for the tensile stage. For the plastic stage, the quasi-brittle material exhibits plastic softening, and the influence of plastic deformation must be considered in the calculation of the bond force, which can be expressed as: where s p indicates the amount of plastic deformation after the bond enters the plastic phase.

Yield Criteria
In compression, once the bond force exceeds the critical force of the elastic limit of the bond, plastic deformation occurs in the bond, resulting in the accumulation of damage. The bond force decreases due to the accumulation of plastic deformation, and the critical force of the bond decreases due to the accumulation of damage until the bond force is less than the critical force of the bond.
In order to accurately calculate the true bond force for quasi-brittle materials undergoing plastic softening behavior in the compression phase, the bond strengths need to be known. Chu et al. [18], in their work, expressed the bond strengths as follows by emulating the strength expressions in the JH-2 model in the classical framework [30]: where p i is the critical force when the bond is undamaged, p f is the critical force when the bond is fully damaged, expressed as: and D represents the cumulative damage from plastic deformation of the bond: where ∆s p is the plastic deformation in one time step, ∑ ∆s p is the accumulated plastic deformation, s e is the elastic compression limit, and s 1 is the plastic-compression deformation limit.
In order to characterize the relationship of the bond in the compressive plastic softening phase, it is necessary to define a function to determine whether the bond force has reached the maximum allowed value; referring to the theory of continuum-media mechanics, the yield criterion of the bond is expressed as: where . s is the compressive deformation rate of the bond, which is the time derivative of the relative stretch of the bond, expressed as follows: . s(ξ, η) = . η(η + ξ) |η||η + ξ| (12)

Flow Rule
Similar to the consistency condition in elastoplastic mechanics, the consistency condition for the bond is defined by deriving Equation (11) as follows: . ϕ s, s p , where the second-order derivative of the relative deformation of the bond is neglected. According to Equation (13), the bond force remains equal to the elastic limit value of the bond. The derivation of Equation (7) is as follows: Based on Equations (13) and (14), the relative compression plastic deformation rate of the bond can be solved as: where H = (1−β)p i s 1 −s e .

Consideration of Thermal Effects
Since the mechanical behavior of the material in the elastic phase is not affected by the loading rate, the thermo-mechanical coupling bond-force function in this phase can be simply expressed in a linear form as follows: where α(T) is the coefficient of thermal expansion under different temperatures; and c(T) is the bond constant considering the effect of temperature, which is a material-dependent constant. This can be obtained by making the elastic strain energy density in the theory of elastic mechanics equal to the deformation energy density at the material point. The relation between c(T) and the elastic modulus of the material under a plane-stress state can be expressed as: where the E(T) is Young's modulus of the material at different temperatures. When the bond is in the plastic phase, the bond-force function can be expressed as follows due to the plastic softening behavior and temperature effects: Meanwhile, the yield function, Equation (13), which determines whether the bond force reaches the maximum permissible value, can be expressed as: where D(T) is the damage that gradually accumulates, considering the thermal effect, could be expressed as: Here, ∆s p (T) is the plastic compression deformation, considering the thermal effect, at each time step ∆t, which can be expressed as: To obtain the change rate of bond stretching . s ξ, η, . η, T , we derived Equation (12): It should be noted that the consistency condition also needs to be satisfied, i.e., the time derivative of Equation (13): The relationship between η, T can be similarly obtained from the above equation and Equation (15) and is expressed as: . s p ξ, η, Therefore, the BB-PD thermo-mechanical coupling equations applicable to quasi-brittle materials are as follows:

Numerical Discretization and Time Integration
For Equation (25), the motion equation and heat conduction can be replaced by a discretized form, as given below: where f ij is the bond force function between the matter points i and j. This can be calculated using Equation (16) in the elastic phase and Equation (18) in the plastic phase.
In this study, in order to describe the thermo-mechanical coupling in the framework of the PD model of quasi-brittle materials, an interleaved scheme is used to approximate the solution. This means that the coupled system equations are solved separately, and different time-explicit algorithms are employed to solve the heat conduction equation and the kinetic equation. In particular, the explicit integration algorithm with the first-order forward difference is used to solve the heat conduction equation, and the temperature of the next time step is obtained as: where ∆t T is the thermo-mechanical time step. In addition, similar to the quasi-static problem, virtual inertia and damping terms are introduced to solve the dynamics equations, which can be expressed as: where D is the virtual diagonal density matrix and c is the damping coefficient, obtained from Reschgorin law and Rayleigh quotient [31], respectively. Using the central difference algorithm, the displacement and velocity for the next time step are defined as: .
where ∆t M is the kinetic time step.
In thermo-mechanical coupled problems, the kinetic characteristic time scale depends on the propagation velocity of the stress wave in the material, and the heat conduction characteristic time scale depends on the thermal diffusivity of the material. In general media, the time scale of heat transfer characteristics is usually much larger than that of kinetic characteristics. Therefore, the whole thermal coupling solution can be divided into the following three steps.
Step 1: The heat conduction equation is solved and the temperature field distribution of the whole model is calculated.
Step 2: The motion equations are solved until the whole model reaches a steady state.
Step 3: The heat conduction equation is solved for the next thermo-mechanical time step. The above steps are repeated to obtain the entire thermo-mechanical coupling solution. It should be noted that the convergence criterion is also needed to determine the steady state of the kinetic iterative solution. When the whole model reaches the steady state, the displacement increment of each kinetic iteration step tends to 0. Chen et al. [14] have previously defined the parametric number, Re, as shown in Equation (32), and provided a minimal value, Ω. When Re ≤ Ω, the system reaches the steady state and can be solved in the next thermo-mechanical time step. Otherwise, the iteration continues until it converges with the formula below: where M is the entire number of model particles.

Model Verification and Convergence Analysis
In this section, the BB-PD model proposed in Section 2 is implemented in Fortran code, and two typical cases are applied to verify the efficiency of the model. The convergence of the numerical model is also analyzed in the following section. Figure 4 shows the computational model of the ceramic plates subjected to heating loads. A flat directional plate with a side length L = W = 1 m is adiabatically constrained to the normal directional displacement with respect to three sides, except for the top. The initial temperature of the whole plate is 0 • C and the temperature of T = 1.0 • C is applied to the top boundary. The peridynamic, mechanical, and thermo-mechanical parameters used in the numerical simulation are listed in Table 1. Three points, referred to as A, B, and C on the vertical symmetry axis at the center of the plate, and located on the top boundary, positive center, and bottom boundary of the plate, respectively, were selected as reference points. media, the time scale of heat transfer characteristics is usually much larger than that of kinetic characteristics. Therefore, the whole thermal coupling solution can be divided into the following three steps.

Ceramic Plates Subjected to Heating Loads
Step 1: The heat conduction equation is solved and the temperature field distribution of the whole model is calculated.
Step 2: The motion equations are solved until the whole model reaches a steady state.
Step 3: The heat conduction equation is solved for the next thermo-mechanical time step.
The above steps are repeated to obtain the entire thermo-mechanical coupling solution. It should be noted that the convergence criterion is also needed to determine the steady state of the kinetic iterative solution. When the whole model reaches the steady state, the displacement increment of each kinetic iteration step tends to 0. Chen et al. [14] have previously defined the parametric number, Re , as shown in Equation (32), and provided a minimal value, Ω . When ≤ Ω Re , the system reaches the steady state and can be solved in the next thermo-mechanical time step. Otherwise, the iteration continues until it converges with the formula below: where M is the entire number of model particles.

Model Verification and Convergence Analysis
In this section, the BB-PD model proposed in Section 2 is implemented in Fortran code, and two typical cases are applied to verify the efficiency of the model. The convergence of the numerical model is also analyzed in the following section. Figure 4 shows the computational model of the ceramic plates subjected to heating loads. A flat directional plate with a side length L = W = 1 m is adiabatically constrained to the normal directional displacement with respect to three sides, except for the top. The initial temperature of the whole plate is 0 °C and the temperature of T = 1.0 °C is applied to the top boundary. The peridynamic, mechanical, and thermo-mechanical parameters used in the numerical simulation are listed in Table 1. Three points, referred to as A, B, and C on the vertical symmetry axis at the center of the plate, and located on the top boundary, positive center, and bottom boundary of the plate, respectively, were selected as reference points.   Thermal parameters

Ceramic Plates Subjected to Heating Loads
An identical finite element model with the same material properties was built using the commercial ABAQUS finite element software. The model was discreted into 200 × 200 grids with a time step of 1 × 10 −5 s, and a direct thermo-mechanical coupling method was used. Additionally, the theoretical calculation formulas for the temperature and vertical displacement of three reference points were provided by Timoshenko [32] and Carslaw [33]: In both simulations, in spite of the temperature loads applied, the temperature change causes deformation inside the plate due to the thermo-mechanical coupling effect, and the reference points are displaced at the same time. The simulation results of different reference points according to both models are shown in Figure 5, the calculation results of PD and ABAQUS are consistent with those of analytical results, thus, demonstrating the reliability of the proposed approach in solving the thermo-mechanical coupling problem.

Thermal parameters
An identical finite element model with the same material properties was built using the commercial ABAQUS finite element software. The model was discreted into 200 200 × grids with a time step of , and a direct thermo-mechanical coupling method was used. Additionally, the theoretical calculation formulas for the temperature and vertical displacement of three reference points were provided by Timoshenko [32] and Carslaw [33]: In both simulations, in spite of the temperature loads applied, the temperature change causes deformation inside the plate due to the thermo-mechanical coupling effect, and the reference points are displaced at the same time. The simulation results of different reference points according to both models are shown in Figure 5, the calculation results of PD and ABAQUS are consistent with those of analytical results, thus, demonstrating the reliability of the proposed approach in solving the thermo-mechanical coupling problem.

Pre-Cracked Brazilian Disk under Uniaxial Compression
The Ayatollahi's experiments on the brittle fracture of polycrystalline graphite were simulated [34]. A modified version of the cracked Brazilian disk (CBD) specimen called the V-notched Brazilian disk (VBD) specimen was used in this experiment. As shown in Figure 6, the specimen is a circular disk of diameter D containing a central rhombic hole with an opening angle 2α and length d for the VBD specimen. The disk diameter and the notch depth were 60 mm and 15 mm, respectively, and the angles used in the experiment were 2α = 30 • , β = 15 • . The basic material properties of polycrystalline graphite are as follows: density of 1710 kg/m 3 , Young's modulus of 8.05 GPa, Poisson's ratio of 0.33, and the fracture toughness of 1.0 MPa m 0.5 . In the experiment, the fracture test was performed by using a universal tension-compression test machine under displacement conditions with a loading rate of 0.05 mm/min. Using the proposed model to simulate the crack extension of the VDB specimen, the splitting damage process is shown in Figure 7. At 70 s, the tips of both sides of the pre-existing crack begin to accumulate damage due to stress concentration exceeding the strength limit, and crack initiation occurs here. Then at 80 s, cracks appear at the tips of both sides of the pre-existing crack and begin to propagate outward. Next, at 110 s, the cracks on both sides remain symmetrical and propagate to both sides of the loading. Finally, at 170 s, cracks penetrate the tips of both ends of the precast crack and both sides of the loading point, at which point the Brazilian disc specimen fails. the notch depth were 60 mm and 15 mm, respectively, and the angles used in the experiment were 2 30 α =°, . In the experiment, the fracture test was performed by using a universal tension-compression test machine under displacement conditions with a loading rate of 0.05 mm/min. Using the proposed model to simulate the crack extension of the VDB specimen, the splitting damage process is shown in Figure 7. At 70 s, the tips of both sides of the preexisting crack begin to accumulate damage due to stress concentration exceeding the strength limit, and crack initiation occurs here. Then at 80 s, cracks appear at the tips of both sides of the pre-existing crack and begin to propagate outward. Next, at 110 s, the cracks on both sides remain symmetrical and propagate to both sides of the loading. Finally, at 170 s, cracks penetrate the tips of both ends of the precast crack and both sides of the loading point, at which point the Brazilian disc specimen fails.   Using the proposed model to simulate the crack extension of the VDB specimen, the splitting damage process is shown in Figure 7. At 70 s, the tips of both sides of the preexisting crack begin to accumulate damage due to stress concentration exceeding the strength limit, and crack initiation occurs here. Then at 80 s, cracks appear at the tips of both sides of the pre-existing crack and begin to propagate outward. Next, at 110 s, the cracks on both sides remain symmetrical and propagate to both sides of the loading. Finally, at 170 s, cracks penetrate the tips of both ends of the precast crack and both sides of the loading point, at which point the Brazilian disc specimen fails.  The comparison of the prefabricated cracked Brazilian disc splitting damage before and after the experiment with the PD simulation results is shown in Figure 8. It can be observed that the simulated results are in high agreement with the experimental results. The comparison of the prefabricated cracked Brazilian disc splitting damage before and after the experiment with the PD simulation results is shown in Figure 8. It can be observed that the simulated results are in high agreement with the experimental results.

Convergence Analysis
Furthermore, the numerical convergence of the model was also analyzed using the case study of Section 3.1. Currently, the convergence analysis for PD consists of two main types: m-convergence and δ -convergence [35].
In the m-convergence, horizon δ is kept constant as , respectively, as shown in Figure 9. The vertical displacements of the reference point obtained usinig the proposed model, the analytical solution, and the FEM simulation are shown in Figure 10. For a fixed horizon, as the value of m increases, the error rate between the simulation results and the analytical results becomes smaller. Although

Convergence Analysis
Furthermore, the numerical convergence of the model was also analyzed using the case study of Section 3.1. Currently, the convergence analysis for PD consists of two main types: m-convergence and δ-convergence [35].
In the m-convergence, horizon δ is kept constant as δ = 5.0 × 10 −3 m throughout the computation, while m and ∆x are taken as 2, 3, and 4, and 2.5 × 10 −3 m, 1.66 × 10 −3 m, and 1.25 × 10 −3 m, respectively, as shown in Figure 9. The vertical displacements of the reference point obtained usinig the proposed model, the analytical solution, and the FEM simulation are shown in Figure 10. For a fixed horizon, as the value of m increases, the error rate between the simulation results and the analytical results becomes smaller. Although the error can be captured when m = 3 or 4, m = 3 is preferred, considering the effect of computational efficiency.

Convergence Analysis
Furthermore, the numerical convergence of the model was also analyzed using the case study of Section 3.1. Currently, the convergence analysis for PD consists of two main types: m-convergence and δ -convergence [35].
In the m-convergence, horizon δ is kept constant as , respectively, as shown in Figure 9. The vertical displacements of the reference point obtained usinig the proposed model, the analytical solution, and the FEM simulation are shown in Figure 10. For a fixed horizon, as the value of m increases, the error rate between the simulation results and the analytical results becomes smaller. Although the error can be captured when 3 or 4 m = , 3 m = is preferred, considering the effect of computational efficiency.

Convergence Analysis
Furthermore, the numerical convergence of the model was also analyzed using the case study of Section 3.1. Currently, the convergence analysis for PD consists of two main types: m-convergence and δ -convergence [35].
In the m-convergence, horizon δ is kept constant as , respectively, as shown in Figure 9. The vertical displacements of the reference point obtained usinig the proposed model, the analytical solution, and the FEM simulation are shown in Figure 10. For a fixed horizon, as the value of m increases, the error rate between the simulation results and the analytical results becomes smaller. Although the error can be captured when 3 or 4 m = , 3 m = is preferred, considering the effect of computational efficiency.  In the δ-convergence, the non-locality parameter m is kept as a constant, i.e., m = 3. Three different horizon sizes are chosen as δ = 1.5 × 10 −2 m, 3 × 10 −2 m, 6 × 10 −2 m, and ∆x = 5 × 10 −3 m, 1 × 10 −2 m, 2 × 10 −2 m (see Figure 11). The vertical displacements calculated at reference points using the different δ are compared to those obtained from the finite element method, as shown in Figure 12. For a fixed non-locality parameter, with a decrease in the horizon sizes, the error calculation will also decrease, but it will lead to an increase in the calculation efficiency. Therefore, choosing the right value of horizon sizes requires special consideration. In the δ -convergence, the non-locality parameter m is kept as a constant, i.e., m = 3.
Three different horizon sizes are chosen as  Figure 11). The vertical displacements calculated at reference points using the different δ are compared to those obtained from the finite element method, as shown in Figure 12. For a fixed non-locality parameter, with a decrease in the horizon sizes, the error calculation will also decrease, but it will lead to an increase in the calculation efficiency. Therefore, choosing the right value of horizon sizes requires special consideration.

Numerical Applications
The previous studies have shown that the proposed model is able to accurately simulate the mechanical behavior of quasi-brittle materials under thermal loading (in Section 3.1) and static loading (in Section 3.2). For the purpose of clarifying the applicability of the model, this section applies it to two more complex coupled thermal-force processes: the ceramic quenching process (in Section 4.1) and the process of compressing the rock after heat treatment (in Section 4.2).

Ceramic under Cold Shock
Referring to the quenching experiments of ceramic plates at different temperatures conducted by Jiang et al. [36], an alumina ceramic plate with the dimensions of 50mm 10mm × is heated to 873 K and subsequently allowed to freefall into the water at 293 K. Considering the symmetry of the load and boundary conditions, a 1/4 model,

Numerical Applications
The previous studies have shown that the proposed model is able to accurately simulate the mechanical behavior of quasi-brittle materials under thermal loading (in Section 3.1) and static loading (in Section 3.2). For the purpose of clarifying the applicability of the model, this section applies it to two more complex coupled thermal-force processes: the ceramic quenching process (in Section 4.1) and the process of compressing the rock after heat treatment (in Section 4.2).

Ceramic under Cold Shock
Referring to the quenching experiments of ceramic plates at different temperatures conducted by Jiang et al. [36], an alumina ceramic plate with the dimensions of 50 mm × 10 mm is heated to 873 K and subsequently allowed to freefall into the water at 293 K. Considering the symmetry of the load and boundary conditions, a 1/4 model, shown in Figure 13, could be established to perform the calculation. In such a model, the left and lower boundaries are constrained during the normal displacements, and a uniform and constant-cold impact load is applied to the upper and right boundaries. The convective heat transfer coefficient h = 70, 000 W/ m 2 · K is taken when the ceramic plate is dropped into the water. The PD and thermo-mechanical parameters involved in the numerical simulations are listed in Table 2.  Table 2.   Figure 13. Schematic diagram of the geometry and boundary condition of the ceramic subjected to cold shock. Thermal parameters [36] Thermal conductivity k T (W · m −1 K −1 ) 31 Coefficient of thermal expansion α (1/K) 7.5 × 10 −6 Specific heat capacity c v (J · Kg −1 K −1 ) 880 When the high-temperature ceramic plate (873 K) enters the room-temperature water (293 K), the heat energy of the plate spreads rapidly to the surrounding ambient medium. The surface temperature of the ceramic plate decreases sharply, as shown in Table 3, forming a vast temperature gradient with the interior. This temperature gradient from the inside to the outside (i.e., hot inside and cold outside) causes the ceramic surface to undergo tensile stress while the interior is subjected to compressive stress. When the tensile stress on the ceramic surface exceeds that of the interior of the ceramics, damage occurs, and cracks propagate throughout the material. Since the temperature bond also breaks due to the thermal effect, heat conduction through the crack is blocked and the temperature on both sides of the bending crack exhibits a significant temperature jump.  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at  For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at For cracks due to thermal shock, all cracks are distributed in a parallel manner at equal distances on the outer surface of the ceramic (upper and right side) and extend from the outside to the inside. As the cold shock continues, some of the initial cold shock cracks stop growing, while other thermal shock cracks continue to grow.
During the first 10 ms, the cracks are uniformly distributed at intervals of roughly 0.001 mm, and the length of each crack remains consistent. As the cold shock continues, some initial cold shock cracks stop growing at 50 ms, while other heat shock cracks continue to grow. Thereafter, the ceramic plate temperature gradient decreases, the thermal stress becomes smaller, the crack expansion slows down, and the crack stops growing at 600 ms, reaching the maximum length.
In the work of Jiang et al. [36], the area within 10 mm of the ends of the specimen was excluded in order to eliminate the effect of the end boundaries. The average dimensionless crack spacing s and dimensionless crack length p were proposed, denoted as s = s/L C and p = p/L C , respectively, where s is the crack spacing, p is the crack length, and L C is the specimen width. The average dimensionless crack spacing in the simulation results is 0.112, compared with 0.12 in experiments, and the dimensionless crack length in the simulation results is 0.715, compared with 0.79 in experiments [36]. The thermal impact cracks show a clear spacing distribution, i.e., there are short cracks in the middle of long cracks. The comparison between the simulated and experimental results is shown in Figure 14, where the thermal impact cracks remain similar in terms of spacing, length, length hierarchy, and periodicity. However, since the model used in the experiments is not an ideal model, the ceramic plate is a non-homogeneous material and there are small gaps in the structure, which cannot be consistent with the simulated results, as evidenced by the asymmetry of the thermal cracks in the experimental results.

Granite under Uniaxial Compression after Heat Treatment
The granite specimen with prefabricated cracks was first subjected to thermocycling and then compressed uniaxially, as performed by Yang et al. [37]. The dimensions of the granite specimen were 80 mm 160 mm × (see Figure 15a), and there was a prefabricated crack with a length of 20 mm , width of 1.5 mm , and inclination angle of 30° in the center of the specimen. At the thermal cycling stage, Yang et al. first heated the granite specimen to 573 K and then kept the temperature constant to make the inside and outside of the sample converge to the same temperature. Subsequently, the sample was placed in the open air and cooled down naturally to room temperature (293 K). At the uniaxial compression stage, the top and bottom ends of the specimens were loaded in compression using a loading speed of 0.1 m s , and the crack nucleation and expansion were observed.
The mechanical and thermo-mechanical parameters in the PD simulation are listed in Table 4.

Granite under Uniaxial Compression after Heat Treatment
The granite specimen with prefabricated cracks was first subjected to thermocycling and then compressed uniaxially, as performed by Yang et al. [37]. The dimensions of the granite specimen were 80 mm × 160 mm (see Figure 15a), and there was a prefabricated crack with a length of 20 mm, width of 1.5 mm, and inclination angle of 30 • in the center of the specimen. At the thermal cycling stage, Yang et al. first heated the granite specimen to 573 K and then kept the temperature constant to make the inside and outside of the sample converge to the same temperature. Subsequently, the sample was placed in the open air and cooled down naturally to room temperature (293 K). At the uniaxial compression stage, the top and bottom ends of the specimens were loaded in compression using a loading speed of 0.1m/s, and the crack nucleation and expansion were observed. The mechanical and thermo-mechanical parameters in the PD simulation are listed in Table 4.
to 573 K and then kept the temperature constant to make the inside and outside of the sample converge to the same temperature. Subsequently, the sample was placed in the open air and cooled down naturally to room temperature (293 K). At the uniaxial compression stage, the top and bottom ends of the specimens were loaded in compression using a loading speed of 0.1 m s , and the crack nucleation and expansion were observed.
The mechanical and thermo-mechanical parameters in the PD simulation are listed in Table 4.   Thermal parameters [37] Thermal conductivity k T (W · m −1 K −1 ) 3.5 Specific heat capacity c v (J · Kg −1 K −1 ) 900 Before the numerical simulation, the same geometric model as that for the granite specimen was established, as shown in Figure 15b. Due to the non-homogeneous properties of rock materials, the Weibull distribution is often introduced to describe the statistical distribution of the characteristic parameters [38], such as elastic modulus, Poisson's ratio, and thermal expansion coefficient in the PD simulation. However, the Weibull distribution does not accurately reflect the properties of granite due to the variety of mineral components and contents of rocks and their vastly different material properties. Therefore, this study consists of reconstructing the PD calculation model of the non-homogeneous granite with the non-uniform and discontinuous thermal expansion coefficients using the Knuth-Durstenfeld stochastic algorithm proposed by Yang et al. [17] (see Figure 15c). The proportions of mineral compositions and thermal expansion coefficients of granite materials are listed in Table 5. After 1000 s of thermal loading, the temperature of the granite specimen increased from 293 K to 573 K. Subsequently, after 3200 s, the sample naturally cooled down to room temperature (293 K). The simulation of the whole heat treatment process is shown in Figure 16a. Due to the slow temperature rise at the warming stage, the temperature difference between the inside and outside of the granite is small, and the non-uniform thermal stress caused by the temperature gradient is minor. Therefore, only a tiny amount of discontinuous thermal cracks is generated inside the granite during the entire heating-up stage. In addition, the higher compression strength of the granite also suppresses the crack generation at the warming stage. Unlike the warming stage, the outer surface of the granite decays sharply to room temperature during the natural cooling stage. This drastic heat transfer behavior leads to the formation of a vast temperature gradient inside and outside the granite, which provokes a rapid increase of the tensile (thermal) stress on the surface of the specimen under tensile strength, causing more discontinuous cracks to occur on both sides of the sample. These cracks continue to expand in the course of the cooling process and then gradually penetrate and fall off. Moreover, due to the inconsistency between the thermal expansion coefficients of different mineral components inside the specimen, there are more and more cracks induced by uneven thermal expansion. The crack initiation and propagation process of the granite specimens containing prefabricated cracks under uniaxial compression simulated through the PD model is shown in Figure 16b. In a future study, the crack types will be analyzed according to the classification of crack types proposed by Yang et al. [38]. At the initial stage of loading, the main strain concentrations are found from the tips of the pre-existing fissure; the granite specimens had no macroscopic crack generation except for a large number of thermal microcracks on both sides, due to the thermal cycling process. With the increase of load, when The crack initiation and propagation process of the granite specimens containing prefabricated cracks under uniaxial compression simulated through the PD model is shown in Figure 16b. In a future study, the crack types will be analyzed according to the classification of crack types proposed by Yang et al. [38]. At the initial stage of loading, the main strain concentrations are found from the tips of the pre-existing fissure; the granite specimens had no macroscopic crack generation except for a large number of thermal micro-cracks on both sides, due to the thermal cycling process. With the increase of load, when the time reached 1 ms, the main strain concentrations develop obviously, the secondary tensile crack appeared at both the upper and lower ends of the prefabricated crack of the specimen, but the development of the secondary tensile crack was not symmetrical due to the uneven distribution of thermal micro-cracks inside the model. Subsequently, at 1.35 ms, a downward expanding tensile wing crack appeared at the upper end of the precast crack, while an upward expanding anti-shear crack appeared at the lower end of the precast crack. At the same time, thermal micro-cracks can be observed developing into secondary tensile cracks on both sides of the prefabricated cracks that afterward form web-like cracks. The macroscopic cracks in the final granite specimens were classified as secondary tensile cracks and anti-shear cracks.
The comparison between the experimental and PD simulation results of uniaxial compression after thermal cycling of precast cracked granite is shown in Figure 17. Secondary tensile cracks and anti-shear cracks existed at the tips of both sides of the precast cracks and were approximately the same in the form of extension. However, both tensile wing cracks and anti-shear cracks exist on the right side of the precast crack in the experimental results [37], while only anti-shear cracks exist in the simulation process. This is caused by the existence of tiny voids inside the granite, which is a non-homogeneous material, and the presence of a large number of non-uniformly distributed thermal micro-cracks during the thermal cycling process, which prevent the anti-tensile cracks that should appear along the axial stress direction; hence, only the anti-shear cracks mainly caused by shear damage appeared. The inconsistency of crack forms on both sides in the experimental results also indicates the non-homogeneous nature of granite. crack, while an upward expanding anti-shear crack appeared at the lower end of the precast crack. At the same time, thermal micro-cracks can be observed developing into secondary tensile cracks on both sides of the prefabricated cracks that afterward form weblike cracks. The macroscopic cracks in the final granite specimens were classified as secondary tensile cracks and anti-shear cracks.
The comparison between the experimental and PD simulation results of uniaxial compression after thermal cycling of precast cracked granite is shown in Figure 17. Secondary tensile cracks and anti-shear cracks existed at the tips of both sides of the precast cracks and were approximately the same in the form of extension. However, both tensile wing cracks and anti-shear cracks exist on the right side of the precast crack in the experimental results [37], while only anti-shear cracks exist in the simulation process. This is caused by the existence of tiny voids inside the granite, which is a non-homogeneous material, and the presence of a large number of non-uniformly distributed thermal microcracks during the thermal cycling process, which prevent the anti-tensile cracks that should appear along the axial stress direction; hence, only the anti-shear cracks mainly caused by shear damage appeared. The inconsistency of crack forms on both sides in the experimental results also indicates the non-homogeneous nature of granite.  Existing studies have shown that the deformation process of concrete is very complex due to its heterogeneity and involves progressive damage, such as the generation, propagation, and coalescence of microcracks [39]. While we have considered the parameters of the different components of the granite in an effort to construct a macroscopic heterogeneous model subjected to thermo-mechanical coupling loads, we neglected the mechanical characteristics at the microscopic scale. Through the references [40][41][42][43][44], it may be observed that the study of thermo-mechanical coupling under microstructures tends to focus on the mechanical properties of microstructures using the theory of nonlocal elasticity, by combining intermolecular or interatomic bonds into their specific intrinsic structural relationships. Further study in this area would increase the validity of the model.

Conclusions
In this paper, a coupled model capable of simulating the thermal-force damage behavior of quasi-brittle materials was proposed based on the bond-based PD theory using Existing studies have shown that the deformation process of concrete is very complex due to its heterogeneity and involves progressive damage, such as the generation, propagation, and coalescence of microcracks [39]. While we have considered the parameters of the different components of the granite in an effort to construct a macroscopic heterogeneous model subjected to thermo-mechanical coupling loads, we neglected the mechanical characteristics at the microscopic scale. Through the references [40][41][42][43][44], it may be observed that the study of thermo-mechanical coupling under microstructures tends to focus on the mechanical properties of microstructures using the theory of nonlocal elasticity, by combining intermolecular or interatomic bonds into their specific intrinsic structural relationships. Further study in this area would increase the validity of the model.

Conclusions
In this paper, a coupled model capable of simulating the thermal-force damage behavior of quasi-brittle materials was proposed based on the bond-based PD theory using the fully thermodynamic coupling equation. The model consists in describing different mechanical properties of quasi-brittle materials in tensile and compressive states, constructing bond force functions in the tensile and compressive phases, and introducing the role of temperature terms in the bond-based peridynamic model. The simulations of the thermal expansion process of ceramics and the static compression damage of polycrystalline graphite were applied to the proposed model, and the numerical model results showed agreement with the experimental results in the references. In addition, the model was also used to simulate thermal damage processes in ceramics and in homogeneous rocks, revealing the potential capacity of the model in analyzing the post-thermal damage behavior of quasi-brittle materials.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the findings of this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.