Discrete Element Modeling of Thermally Damaged Sandstone Containing Two Pre-Existing Flaws at High Conﬁning Pressure

: An underground coal gasiﬁcation (UCG) process is strongly exothermic, which will cause thermal damage on rock cap. We proposed a new thermal damage numerical model based on a two dimension particle ﬂow code (PFC2D) to analyze the inception and extension of cracks on pre-cracked red sandstone, which were thermally treated at a temperature of 25~1000 ◦ C. The results indicated that: (1) a thermal damage value D T obtained by extracting the thermal crack area of scanning electron microscope (SEM), which can be used as an indicator of the degree of thermal damage of the sandstone; (2) a thermal damage numerical model established by replacing the ﬂat-joint model with the smooth-joint model based on the thermal damage value D T , this approach can properly simulate the mechanical behavior and failure patterns of sandstone; (3) the critical temperature for strength reduction was 750 ◦ C. The peak strength increased as pre-treatment temperature increased from 25 to 750 ◦ C and then decreased. The elastic modulus E 1 decreased with the increasing thermal treatment temperature; (4) micro-scale cracks initiate from the tip of the prefabricated ﬁssure, and expand along the direction of prefabricated ﬁssure, ﬁnally developing into macroscopic fracture. This approach has the potential to enhance the predictive capability of modeling and presents a reliable model to simulate the mechanical behavior of thermally damaged sandstones, thereby offering a sound scientiﬁc basis for the utilization of space after UCG.


Introduction
Underground coal gasification (UCG) was identified as an environmentally friendly technique for in situ gasification of deep unmineable coal seams, which can expand our energy resources in a sustainable way; additionally, the captured carbon dioxide can be stored in the space after UCG [1,2].However, the UCG process generates a significant amount of heat, causing temperatures in the burn zone to occasionally exceed 900 • C [3,4], which significantly surpasses the critical point (T c ~400-500 • C) that triggers alterations in the mechanical properties of common sandstones [5].On the other hand, the mechanical response of rock is more complex when thermal induced cracks combined with pre-existing fissures [6].
Laboratory-scale specimens were extensively used to investigate how temperature affects the properties of rocks [7][8][9][10][11][12][13][14].These studies mainly focused on physical structure and mechanical response of thermal damaged rock.The temperature ranges from ambient to 400 • C and leads to the evaporation and release of water that is adhered to, combined with, or structurally bound to the sandstone.Between 400 and 600 • C, the minerals in the sandstone undergo thermal reactions and transformation, leading to an increase in porosity, a decrease in permeability and diffusivity, and changes in heat capacity [15].For red sandstones with pre-existing fissures, the elastic modulus and strength increase at 300 • C but decrease as the temperature rises to 900 • C, while the Poisson's ratio remains unaffected by temperature changes.As the temperature increases, the crack closure stress, crack initiation stress, and damage threshold first increase and then decrease [5].In addition, high-temperature damage (HTD) in sandstone induces considerably more acoustic emission (AE) signals [11].
The use of two-or three-dimensional discontinuum models in distinct element methods is a powerful method to simulate various aspects of rock deformation and failure [16][17][18][19], This approach allows for predicting mechanical behavior and tracking fracture evolution under different stress conditions [20].To account for thermal strain caused by temperature changes, disk radii are adjusted and bond forces are modified accordingly [21].In order to accurately capture the thermal behavior of rock specimens, different thermal expansion coefficients were assigned to various minerals [22,23].According to this method, a mesostructured-based numerical model adopted to analyze rock thermal cracking based on elastic damage mechanics and thermal-elastic theory [24,25].The mechanisms responsible for temperature-dependent mechanical properties of thermal damaged modeling granite [22].A cluster model in two-dimension particle flow code (PFC 2D ) was used to simulate a crystalline granite specimen, and crack number variation and crack evolution process was investigated [14].However, this approach had some shortcomings.Firstly, it is usually applicable to crystalline rocks (such as granite); however, the mineral grains of crystalline rocks are difficult to distinguish.Secondly, the heat temperature should be below 500 • C, because temperatures exceeding this threshold may trigger intricate mineral decomposition and other chemical reactions, such as α-β quartz phase transition, which cannot be accurately simulated by particle mechanics method [22].In addition, it is difficult to precisely determine the thermal expansion coefficient of different minerals [23].
The mechanical response of fractured rock masses is usually different from that of intact rock masses [26][27][28].Studies showed that stress concentrations are more likely to occur at the tips of fissures, thus promoting the initiating of new cracks and reducing the strength of the rock mass [29].Previous studies fully investigated intact rock masses as well as rock masses containing single fissures [30], but rock masses containing two pre-existing fissures, especially the interaction mechanism between fissures and fissures and between fissures and thermal cracks, need further research.Therefore, in this paper, the double-fissure specimens with 45 • fissure inclination were subjected to heat treatment at different temperatures in order to determine the influence of thermal damage on the mechanical behavior of fissured red sandstone.Furthermore, we proposed a new approach to establish two-dimension PFC thermal damage numerical model by replacing the flatjoint model with the smooth-joint model according to the thermal damage value D T .The main mechanism of failure and damage of prefabricated fractured sandstone specimens after heat treatment was explored by comparing and analyzing the numerical simulation results and the physical experimental results.

Preparation for Dual Fissured Sandstone
We used red sandstone from Rizhao City, Shandong Province, China, for our experiments.The sandstone had a dry density of 2402 kg/m 3 and an effective porosity of approximately 6.26 to 6.55% [14].The sandstone was primarily composed of quartz, feldspar, dolomite, hematite, and clay minerals [20] and the microstructure of red sandstone is shown in Figure 1.The samples were cut and polished to form cylinders with diameters of 50 mm and heights of 100 mm in according with ISRM [31].Subsequently, two nearly identical fissures with lengths (2a = 12 mm) and widths (d = 2 mm) were cut into the samples using a diamond wire saw.The two coplanar fissures had an inclination angle of 45º and were separated by an intervening rock bridge of 2b = 12 mm (Figure 2).The samples were dried for approximately 24 h and were then heated to 250, 500, 750, and 1000 • C at a rate of 5 • C per minute using a high-temperature furnace (Figure 3a), with a sample at room temperature 25 • C serving as the control group.The specimens were held in furnace at their target temperature for two hours before being allowed to cool naturally to room temperature (Figure 3b).The color of the samples changed during the heat treatment process, progressing from a drab-red color at 250 • C to increasingly redder hues at 500 and 750 • C, ultimately becoming brown at 1000 • C (as depicted in Figure 2).This was due to the gradual transformation of iron-containing minerals into hematite as a result of exposure to high temperature [32].The basic physical parameters of the sandstone specimens are presented in Table 1.
shown in Figure 1.The samples were cut and polished to form cylinders with dia of 50 mm and heights of 100 mm in according with ISRM [31].Subsequently, two identical fissures with lengths (2a = 12 mm) and widths (d = 2 mm) were cut into th ples using a diamond wire saw.The two coplanar fissures had an inclination angle and were separated by an intervening rock bridge of 2b = 12 mm (Figure 2).The sa were dried for approximately 24 h and were then heated to 250, 500, 750, and 1000 rate of 5 °C per minute using a high-temperature furnace (Figure 3a), with a sam room temperature 25 °C serving as the control group.The specimens were held in f at their target temperature for two hours before being allowed to cool naturally to temperature (Figure 3b).The color of the samples changed during the heat treatme cess, progressing from a drab-red color at 250 °C to increasingly redder hues at 5 750 °C, ultimately becoming brown at 1000 °C (as depicted in Figure 2).This was the gradual transformation of iron-containing minerals into hematite as a result of sure to high temperature [32].The basic physical parameters of the sandstone spec are presented in Table 1.shown in Figure 1.The samples were cut and polished to form cylinders with dia of 50 mm and heights of 100 mm in according with ISRM [31].Subsequently, two identical fissures with lengths (2a = 12 mm) and widths (d = 2 mm) were cut into th ples using a diamond wire saw.The two coplanar fissures had an inclination angle and were separated by an intervening rock bridge of 2b = 12 mm (Figure 2).The sa were dried for approximately 24 h and were then heated to 250, 500, 750, and 1000 rate of 5 °C per minute using a high-temperature furnace (Figure 3a), with a sam room temperature 25 °C serving as the control group.The specimens were held in fu at their target temperature for two hours before being allowed to cool naturally to temperature (Figure 3b).The color of the samples changed during the heat treatmen cess, progressing from a drab-red color at 250 °C to increasingly redder hues at 5 750 °C, ultimately becoming brown at 1000 °C (as depicted in Figure 2).This was the gradual transformation of iron-containing minerals into hematite as a result of sure to high temperature [32].The basic physical parameters of the sandstone spec are presented in Table 1.

Experiment Apparatus and Method
To account for the high ground stress experienced during the UCG operation, t confining pressure was set to 30 MPa in triaxial compression experiments.This cor sponded to a depth of 1270 m based on the calculated ground stress using the bulk weig (approximately 23.5 kN/m 3 ) of the red sandstone.The conventional triaxial compressi experiments were performed using a rock triaxial rheology testing system equipped w three servo-hydraulic control units for applying: axial stress, confining pressure, and po pressure loads, respectively.The system had a maximum axial loading capacity of 800 k and a maximum confining pressure and pore pressure of 60 MPa.The axial deformati of the specimen was monitored using a pair of linear variable displacement transduc (LVDTs), while the circumferential deformation was monitored using a circumferent extensometer.
The testing procedures were as follows: (ⅰ) the samples were subjected to confini pressure of 30 MPa at a rate of 4 MPa/min; (ⅱ) the samples were subjected to axial loadi at a rate of 0.04 mm/min until they reached their residual strength.In this manner, t stress-strain curves and triaxial compressive strengths (TCS) of the samples could be o tained.It should be noted that the specimens used in the test contained two prefabricat fissures, the confining pressure will break the rubber sleeve at the location of the fissu and cause the hydraulic oil to enter the specimens, which will affect the mechanical pro erties of the specimen.Therefore, before installing the rubber sleeve, a piece of iron w affixed to cover the fissures with hot-melt adhesive to prevent the hydraulic oil breaki the rubber sleeve into the specimen.However, placing the iron piece cover on the fissu caused a slight bulge in the middle of the specimen; to prevent this bulge, the circumf ential displacement sensor can only be installed in the lower and middle of the specim which resulted in a minor lateral deformation being measured.

Thermal Damage Analysis Based on Image Recognition
Scanning electron microscopy (SEM) is usually used for observing the microstructu of rocks, and there are few quantitative studies on the distribution and proportion

Experiment Apparatus and Method
To account for the high ground stress experienced during the UCG operation, the confining pressure was set to 30 MPa in triaxial compression experiments.This corresponded to a depth of 1270 m based on the calculated ground stress using the bulk weight (approximately 23.5 kN/m 3 ) of the red sandstone.The conventional triaxial compression experiments were performed using a rock triaxial rheology testing system equipped with three servo-hydraulic control units for applying: axial stress, confining pressure, and pore pressure loads, respectively.The system had a maximum axial loading capacity of 800 kN, and a maximum confining pressure and pore pressure of 60 MPa.The axial deformation of the specimen was monitored using a pair of linear variable displacement transducers (LVDTs), while the circumferential deformation was monitored using a circumferential extensometer.
The testing procedures were as follows: (I) the samples were subjected to confining pressure of 30 MPa at a rate of 4 MPa/min; (ii) the samples were subjected to axial loading at a rate of 0.04 mm/min until they reached their residual strength.In this manner, the stressstrain curves and triaxial compressive strengths (TCS) of the samples could be obtained.It should be noted that the specimens used in the test contained two prefabricated fissures, the confining pressure will break the rubber sleeve at the location of the fissure and cause the hydraulic oil to enter the specimens, which will affect the mechanical properties of the specimen.Therefore, before installing the rubber sleeve, a piece of iron was affixed to cover the fissures with hot-melt adhesive to prevent the hydraulic oil breaking the rubber sleeve into the specimen.However, placing the iron piece cover on the fissures caused a slight bulge in the middle of the specimen; to prevent this bulge, the circumferential displacement sensor can only be installed in the lower and middle of the specimen, which resulted in a minor lateral deformation being measured.

Thermal Damage Analysis Based on Image Recognition
Scanning electron microscopy (SEM) is usually used for observing the microstructure of rocks, and there are few quantitative studies on the distribution and proportion of micro-cracks, especially micro-cracks caused by heat treatment.Therefore, based on image recognition methods, we extracted micro-cracks from sandstone with different thermal damage in SEM, and conducted a quantitative analysis to determine the influence of heat temperature on thermal cracks distribution.The image threshold segmentation method is a common technical method in the image recognition process, which is simple to operate and implement [34].The image was composed of two areas with different gray levels, bright and dark, denoted as f (x, y), such as cracks and matrix of rock.The bright and dark parts of this kind of image can be distinguished by the gray value histogram, so a threshold can be selected to distinguish bright and dark peaks.Threshold segmentation can be regarded as a function operation.
where x and y is coordinates of pixel, p(x, y) is local features of pixels, f (x, y) is gray value of pixels, the thresholder image is defined as follows: where the pixel grayscale value of 1 corresponds to the target object (cracks), and the pixel grayscale value of 0 corresponds to the background (rock block).
According to this method, the sandstone samples heat-treated at different temperatures were scanned by electron microscope at a magnification of 500 times.Then, a histogram was applied thresholding to the SEM images based to their grayscale value.Finally, an appropriate gray threshold was set to identify and extract the crack area.Due to the gray scale difference between the cracks and the skeleton in the SEM image of rock, according to the visual inspection of the histogram, the gray threshold of the cracks was set to 40~60, and then, the cracks were identified and extracted, as shown in Figure 4 (red area).
SEM images and extracted crack results of the red sandstone after different pretreatment temperatures are illustrated in Figure 4.In the SEM images of pre-treatment of 25 • C and 250 • C, the grain boundaries and micro voids were clearly visible, and some boundary cracks could also be observed (Figure 4a,b).From 500 • C, the mineral particles inside the rock began to expand due to thermal expansion (Figure 4c).Different mineral particles have different thermal expansion properties, which leads to the dislocation cracks produced [6].The thermal decomposition of calcite results in the formation of certain transgranular cracks (Figure 4d), especially when the heat treatment temperature continues to increase to 1000 • C, a large number of transgranular cracks and boundary cracks appeared (Figure 4e), which indicated that the thermal damage increased significantly.Furthermore, according to the number of pixels of the identified cracks, the crack areas under different pre-treatment temperature can be calculated; they were 1486 mm 2 , 1611 mm 2 , 2338 mm 2 , 3114 mm 2 , and 4507 mm 2 , respectively.The sample without thermal treatment (room temperature) is considered a comparative sample and has no thermal damage, its crack area is regarded as the baseline crack area of the sandstone.For thermal treatment samples, the thermal crack area can be obtained by the total crack area subtracting the original crack area.The damage value D T can be calculated by dividing the thermal crack area by the SEM image area, as shown in the following: where S T represents the thermal crack area, S S represents the SEM image area, measuring 56,100 mm 2 in total.D T can be characterizing the degree of thermal damage of the heat treatment rock samples.The D T values for varying treatment temperatures are calculated according to Equation (3), which correspond to 0%, 0.22%, 1.52%, 2.90%, and 5.39%, respectively.Figure 5 shows the relationship between D T and E against heat treatment temperature (the elastic modulus E presented in Yang et al. 2021).As the temperature increases, the thermal damage value D T gradually increases, and resulting in a decrease in the elastic modulus E. This provides further evidence that the thermal cracks caused the rock sample to become softened.
The DT values for varying treatment temperatures are calculated according to Equation (3), which correspond to 0%, 0.22%, 1.52%, 2.90%, and 5.39%, respectively.Figure 5 shows the relationship between DT and E against heat treatment temperature (the elastic modulus E presented in Yang et al. 2021).As the temperature increases, the thermal damage value DT gradually increases, and resulting in a decrease in the elastic modulus E. This provides further evidence that the thermal cracks caused the rock sample to become softened.

Particle Flow Code
The mineral properties, structures, deformation, and strength parameters of rock are sensitive to high temperature treatment.The triaxial compressive strength showed an increase as the temperature was raised from 250 to 750 • C, but experienced a significant Sustainability 2023, 15, 6318 8 of 17 decrease once the temperature reached 1000 • C. The main reason for this phenomenon was attributed to the vaporization and escape of the adhered water, combined water and structural water, which occurs between the temperatures of 250 and 750 • C. As a result, this leads to an increase in the coefficient of internal friction, thereby increasing the strength of the rock specimen [5].When the temperature reaches 1000 • C, the progressive accumulation of thermal damage leads to the development of both boundary and transgranular cracks, in addition to the widening of intergranular cracks; consequently, this contributes to a subsequent reduction in strength [14].It is difficult to simulate the aforementioned complex mechanical behaviour of rock.The PFC software is a discontinuum programs to simulate rock mechanical behavior by employing the distinct element method [17].In this study, we used two-dimension PFC numerical mode to reproduce the laboratory triaxial compression experiment on heat treatment sandstone specimen.The flat-joint model was adopted to bond the non-uniform-sized circular particles, because this model can represent the mechanical behavior of rock and obtained reasonable ratio of direct-tension strength to unconfined-compressive strength than the linear parallel bond model [35].The smoothjoint contact model can convert the rotation between particles on two sides of the joint into relative sliding (Figure 6) [36]; hence, the smooth-joint contact model was employed to replace the flat-joint model bonds proportionally according to D T value, in order to simulate the micro thermally induced cracks.
The mineral properties, structures, deformation, and strength parameters of rock sensitive to high temperature treatment.The triaxial compressive strength showed an crease as the temperature was raised from 250 to 750 °C, but experienced a signifi decrease once the temperature reached 1000 °C.The main reason for this phenome was attributed to the vaporization and escape of the adhered water, combined water structural water, which occurs between the temperatures of 250 and 750 °C.As a re this leads to an increase in the coefficient of internal friction, thereby increasing strength of the rock specimen [5].When the temperature reaches 1000 °C, the progres accumulation of thermal damage leads to the development of both boundary and tr granular cracks, in addition to the widening of intergranular cracks; consequently, contributes to a subsequent reduction in strength [14].It is difficult to simulate the af mentioned complex mechanical behaviour of rock.The PFC software is a discontinu programs to simulate rock mechanical behavior by employing the distinct elem method [17].In this study, we used two-dimension PFC numerical mode to reproduce laboratory triaxial compression experiment on heat treatment sandstone specimen.flat-joint model was adopted to bond the non-uniform-sized circular particles, beca this model can represent the mechanical behavior of rock and obtained reasonable r of direct-tension strength to unconfined-compressive strength than the linear par bond model [35].The smooth-joint contact model can convert the rotation between p cles on two sides of the joint into relative sliding (Figure 6) [36]; hence, the smooth-j contact model was employed to replace the flat-joint model bonds proportionally acc ing to DT value, in order to simulate the micro thermally induced cracks.

Set-Up of Thermal Damage Numerical Model
First, circular particles with a radius of 0.3~0.498mm were randomly distributed 50 × 100 mm square area, which produced an intact sample with a porosity of 6.4% a density of 2400 kg•m −3 .To improve the calculation efficiency and minimize the influe of floating particles, a function was compiled to identify and delete particles with fe than two contact bonds.Next, particles within a designated area were removed to cr numerical samples of pre-cracked sandstone with an inclination angle of 45°.The existing fissure had a length of 12 mm and a width of 2 mm, and the length of the bridge was 12 mm, identical to the sandstone specimen used in laboratory testing.final numerical sample contained 8820 particles with 23,756 flat-joint contact bonds nally, the flat-joint bonds were randomly and proportionally replaced with smooth-j bonds according to the thermal damage value DT to simulate the thermal cracking, res ing in the numerical thermal damage sample shown in Figure 7 (smooth-joint bonds shown in blue).

Set-Up of Thermal Damage Numerical Model
First, circular particles with a radius of 0.3~0.498mm were randomly distributed in a 50 × 100 mm square area, which produced an intact sample with a porosity of 6.4% and a density of 2400 kg•m −3 .To improve the calculation efficiency and minimize the influence of floating particles, a function was compiled to identify and delete particles with fewer than two contact bonds.Next, particles within a designated area were removed to create numerical samples of pre-cracked sandstone with an inclination angle of 45 • .The pre-existing fissure had a length of 12 mm and a width of 2 mm, and the length of the rock bridge was 12 mm, identical to the sandstone specimen used in laboratory testing.The final numerical sample contained 8820 particles with 23,756 flat-joint contact bonds.Finally, the flat-joint bonds were randomly and proportionally replaced with smooth-joint bonds according to the thermal damage value D T to simulate the thermal cracking, resulting in the numerical thermal damage sample shown in Figure 7 (smooth-joint bonds are shown in blue).

Calibrating Micro-Mechanical Parameters
In the process of parameter calibration, the numerical simulation results w pared with the experimental results, and then, the micro parameters were grad justed to make the simulation results approach the laboratory test results; this p called trial and error matching.Previous studies determined the correlation betw macro-and micro-parameters; for example, fj_ten and fj_coh mainly affect streng has a linear relationship with elastic modulus, krat affects lateral deformation, a affects the ductility and brittleness of post-peak failure [37].
When calibrating the micro-parameters of the numerical simulation for therm ment samples, two aspects need to be considered.Firstly, the strengthening effe sample strength after undergoing thermal treatment.The laboratory results sho the strength increased linearly with the increase in treatment temperature in the 25~750 °C.The previous studies indicated that the parameter fj_coh in the flat-joi has a linear relationship with the strength of the rock.Therefore, we fixed other ters and increased fj_coh accordingly to obtain results that were consistent strengthening effect of heat treatment (25~750 °C) on mechanical properties.S when the heat treatment exceeded 750 °C and reached 1000 °C, there was a si drop in rock strength.The decrease in strength was not solely caused by the in thermal cracking.We also assumeed that the cohesion will decrease.Therefore, w cohesion to be the same as at 250 °C.Additionally, to simulate thermal cracks, the smooth-joint model, which randomly and proportionally replaced the flat-joi based on the thermal damage value DT.We calibrated the micro-properties of the joint model in the numerical model against the laboratory experiments results.U method, the micro-scale-parameters of flat-joint model and smooth-joint model in damaged modeling specimens are listed in Table 2 and Table 3, respectively.

Calibrating Micro-Mechanical Parameters
In the process of parameter calibration, the numerical simulation results were compared with the experimental results, and then, the micro parameters were gradually adjusted to make the simulation results approach the laboratory test results; this process is called trial and error matching.Previous studies determined the correlation between the macro-and micro-parameters; for example, fj_ten and fj_coh mainly affect strength, emod has a linear relationship with elastic modulus, krat affects lateral deformation, and fj_fric affects the ductility and brittleness of post-peak failure [37].
When calibrating the micro-parameters of the numerical simulation for thermal treatment samples, two aspects need to be considered.Firstly, the strengthening effect on the sample strength after undergoing thermal treatment.The laboratory results showed that the strength increased linearly with the increase in treatment temperature in the range of 25~750 • C. The previous studies indicated that the parameter fj_coh in the flat-joint model has a linear relationship with the strength of the rock.Therefore, we fixed other parameters and increased fj_coh accordingly to obtain results that were consistent with the strengthening effect of heat treatment (25~750 • C) on mechanical properties.Secondly, when the heat treatment exceeded 750 • C and reached 1000 • C, there was a significant drop in rock strength.The decrease in strength was not solely caused by the increase in thermal cracking.We also assumeed that the cohesion will decrease.Therefore, we set the cohesion to be the same as at 250 • C. Additionally, to simulate thermal cracks, we used the smooth-joint model, which randomly and proportionally replaced the flat-joint model based on the thermal damage value D T .We calibrated the micro-properties of the smooth joint model in the numerical model against the laboratory experiments results.Using this method, the micro-scale-parameters of flat-joint model and smooth-joint model in thermal damaged modeling specimens are listed in Tables 2 and 3, respectively.

Validation Tests
Before conducting numerical simulation tests, we validated the previously proposed method for simulating thermal cracking and confirmed its effectiveness in simulating thermal damage.Keeping the parameters of the flat-joint model uniform, we set up two sets of numerical tests.The first set consisted of a numerical model containing thermal cracks by replacing a part of the flat-joint model with smooth-joint model based on the thermal damage value D T .The second set consisted of a numerical model without thermal cracks.The compression results of the two groups of tests are presented in Figure 8.The numerical specimen without thermal cracks exhibited an almost linear increase in strength with thermal treatment temperature in the range of 0~750 • C, while the numerical specimen with thermal cracks showed no significant effect on strength at 250 • C due to the small number of thermal cracks.However, as the thermal treatment temperature increased, resulting in a difference in strength between the numerical model without thermal cracks and the one with thermal cracks.The difference in strength was 15, 31, and 37 MPa at thermal treatment temperature of 500 • C, 750 • C, and 1000 • C, respectively.Based on the above analysis, we concluded that the smooth-joint model replaced the flat-joint contact model, enabling simulation of the phenomenon of increasing thermal damage caused by the rise in thermal treatment temperature.

Validation Tests
Before conducting numerical simulation tests, we validated the previously proposed method for simulating thermal cracking and confirmed its effectiveness in simulating thermal damage.Keeping the parameters of the flat-joint model uniform, we set up two sets of numerical tests.The first set consisted of a numerical model containing thermal cracks by replacing a part of the flat-joint model with smooth-joint model based on the thermal damage value DT.The second set consisted of a numerical model without thermal cracks.The compression results of the two groups of tests are presented in Figure 8.The numerical specimen without thermal cracks exhibited an almost linear increase in strength with thermal treatment temperature in the range of 0~750 °C, while the numerical specimen with thermal cracks showed no significant effect on strength at 250 °C due to the small number of thermal cracks.However, as the thermal treatment temperature increased, resulting in a difference in strength between the numerical model without thermal cracks and the one with thermal cracks.The difference in strength was 15, 31, and 37 MPa at thermal treatment temperature of 500 °C, 750 °C, and 1000 °C, respectively.Based on the above analysis, we concluded that the smooth-joint model replaced the flat-joint contact model, enabling simulation of the phenomenon of increasing thermal damage caused by the rise in thermal treatment temperature.

Mechanical Behavior
Numerical simulations were performed on pre-cracked specimens containing varied number of thermal cracks and calibrated fine-scale parameters.The deviatoric stressstrain curves were compared with the curves obtained from laboratory tests (Figure 9), where the solid line represents the numerical simulation results and the dashed line represents the laboratory tests results.In the initial stage, the dashed line did not exhibit the typical upper concave shape due to the high confining pressure that compressed the initial microcracks in the rock specimen.Thus, linear elastic deformation occurred at the early stage of loading.As the loading continued and approached the first peak strength, the

Mechanical Behavior
Numerical simulations were performed on pre-cracked specimens containing varied number of thermal cracks and calibrated fine-scale parameters.The deviatoric stress-strain curves were compared with the curves obtained from laboratory tests (Figure 9), where the solid line represents the numerical simulation results and the dashed line represents the laboratory tests results.In the initial stage, the dashed line did not exhibit the typical upper concave shape due to the high confining pressure that compressed the initial microcracks in the rock specimen.Thus, linear elastic deformation occurred at the early stage of loading.As the loading continued and approached the first peak strength, the curve growth rate slowed down and transitioned into the plastic deformation stage.A stress drop occurred after reaching the first peak strength due to the closure of the pre-existing cracks, resulting in an increase in the load capacity of the specimen.This closure of prefabricated cracks lead to a second peak strength that was consistently larger than the first peak strength.At a thermal treatment temperature of 1000 • C, after multiple fluctuations near the peak, the experimental curve did not exhibit a significant stress drop (Figure 9e), mainly due to the presence of numerous thermally induced cracks.The gradual closure and expansion of these cracks during compression released some energy, and ultimately, the energy released by the fracture of the pre-existing cracks at the peak was relatively small.We conducted a further investigation on how thermal treatment temperature affects the elastic moudulus and peak deviatoric stress, as shown in Figure 10.Before reaching 750 °C, there was only a small difference (about 1~2 GPa) between the elastic modulus of physical experiment and numerical simulations.However, this difference increased significantly to 9 GPa at the 1000 °C.Although both physical experiments and numerical simulations showed a decrease in elastic modulus with an increase in thermal treatment Due to the mechanical properties of the specimen after the first peak strength having a large gap with those before the peak strength, it was unrealistic to simulate two different mechanical behaviors with the same set of fine-scale parameters.Therefore, this study mainly focused on simulating the mechanical properties of the rock at the first peak strength.The simulated curve basically overlaps with the laboratory test curves in the range of 0~500 • C for both deviatoric stress-axial strain curves and the deviatoric stress-lateral strain curves (Figure 9a-c).However, the deviatoric stress-axial strain curves obtained from simulations thermal treatment temperatures of 750 • C and 1000 • C lie above the laboratory test curves (Figure 9d,e), indicating a higher elastic modulus than that obtained from laboratory tests.
We conducted a further investigation on how thermal treatment temperature affects the elastic moudulus and peak deviatoric stress, as shown in Figure 10.Before reaching 750 • C, there was only a small difference (about 1~2 GPa) between the elastic modulus of physical experiment and numerical simulations.However, this difference increased significantly to 9 GPa at the 1000 • C.Although both physical experiments and numerical simulations showed a decrease in elastic modulus with an increase in thermal treatment temperature, the reduction in numerical simulation after 750 • C was smaller than that in the physical experiment, especially at 1000 • C (Figure 10a).The main reason for this difference is that the physical specimen produced many thermal cracks after high temperature thermal treatment, which closed during the compression process.However, the numerical simulation, where the flat-joint contact was replaced by smooth-joint contact, did not produce any gap, resulting in a smaller elastic modulus than that observed in the physical experiments after 750 • C. Peak deviatoric stress in the numerical experiment was similar to that in the physical experiment at the same temperature (Figure 10b).Both of them increased as the thermal treatment temperature increased from 250 to 750 • C and then, decreased substantially as temperatures reached 1000 • C.These results demonstrate that our proposed thermal damage numerical model can properly simulate the mechanical behavior of sandstone specimens after thermal treatment.We conducted a further investigation on how thermal treatment temperature affects the elastic moudulus and peak deviatoric stress, as shown in Figure 10.Before reaching 750 °C, there was only a small difference (about 1~2 GPa) between the elastic modulus of physical experiment and numerical simulations.However, this difference increased significantly to 9 GPa at the 1000 °C.Although both physical experiments and numerical simulations showed a decrease in elastic modulus with an increase in thermal treatment temperature, the reduction in numerical simulation after 750 °C was smaller than that in the physical experiment, especially at 1000 °C (Figure 10a).The main reason for this difference is that the physical specimen produced many thermal cracks after high temperature thermal treatment, which closed during the compression process.However, the numerical simulation, where the flat-joint contact was replaced by smooth-joint contact, did not produce any gap, resulting in a smaller elastic modulus than that observed in the physical experiments after 750 °C.Peak deviatoric stress in the numerical experiment was similar to that in the physical experiment at the same temperature (Figure 10b).Both of them increased as the thermal treatment temperature increased from 250 to 750 °C and then, decreased substantially as temperatures reached 1000 °C.These results demonstrate that our proposed thermal damage numerical model can properly simulate the mechanical behavior of sandstone specimens after thermal treatment.

Failure Pattern
We also analyzed the damage patterns of numerical specimens subjected to different thermal treatment temperatures from a microscopic perspective and compared them with the physical specimens, as shown in Figure 11.In the numerical specimens, the black microcracks represented tensile cracks, while the red micro-cracks represented shear cracks.These micro-cracks gathered, expanded, and gradually formed macro-cracks.Macroscopic shear crack bands were produced at each end of the two prefabricated fissures, and the combination of macroscopic shear cracks and prefabricated fissures formed a crack band through the upper and lower ends of the specimen, similar to the damage pattern observed Sustainability 2023, 15, 6318 13 of 17 in the physical specimen.A large number of tensile cracks gathered on two sides of the prefabricated fissures, indicating that the prefabricated fissures were mainly dominated by tensile damage.As the thermal treatment temperature increased, the number of microcracks gradually decreased, and the width of the macroscopic shear fracture tended to narrow.The damage patterns of numerical simulations were in good agreement with results of the physical experiment.
We also analyzed the damage patterns of numerical specimens subjected to differen thermal treatment temperatures from a microscopic perspective and compared them with the physical specimens, as shown in Figure 11.In the numerical specimens, the black mi cro-cracks represented tensile cracks, while the red micro-cracks represented shear cracks These micro-cracks gathered, expanded, and gradually formed macro-cracks.Macro scopic shear crack bands were produced at each end of the two prefabricated fissures, and the combination of macroscopic shear cracks and prefabricated fissures formed a crack band through the upper and lower ends of the specimen, similar to the damage pattern observed in the physical specimen.A large number of tensile cracks gathered on two side of the prefabricated fissures, indicating that the prefabricated fissures were mainly domi nated by tensile damage.As the thermal treatment temperature increased, the number o microcracks gradually decreased, and the width of the macroscopic shear fracture tended to narrow.The damage patterns of numerical simulations were in good agreement with results of the physical experiment.Figure 12 shows the displacement of particles in ultimate failure numerical modeling the displacement of the particles was essentially similar and the particles located on two sides of the pre-existing fissures experienced the largest displacement due to the compres sion-induced closure.Following these particles, those particles located at the upper and lower ends of the specimen experienced some displacement as well, while the particle a the ends of the prefabricated fracture (indicated by the blue arrow) underwent minima displacement.Overall, the particles experienced relative displacement along the prefabri cated fissure direction.At the outer end of the prefabricated fissures, the direction of dis placement changed from parallel to fissure to toward the free surface, which caused sur face spalling in physical specimens.On the other hand, at the two sides of th Figure 12 shows the displacement of particles in ultimate failure numerical modeling; the displacement of the particles was essentially similar and the particles located on two sides of the pre-existing fissures experienced the largest displacement due to the compression-induced closure.Following these particles, those particles located at the upper and lower ends of the specimen experienced some displacement as well, while the particle at the ends of the prefabricated fracture (indicated by the blue arrow) underwent minimal displacement.Overall, the particles experienced relative displacement along the prefabricated fissure direction.At the outer end of the prefabricated fissures, the direction of displacement changed from parallel to fissure to toward the free surface, which caused surface spalling in physical specimens.On the other hand, at the two sides of the prefabricated fissures, the displacement direction shifted from being parallel to the fissure to vertical orientation, exhibiting a closure of fissure.

Crack Evolution
There were small differences in damage patterns of specimens due to consistent location, number, and inclination angle of prefabricated fissures.Hence, we selected the thermal damage numerical specimens at 25 • C and 750 • C to illustrate the evolution of deviatoric stress and the number of cracks versus axial strain.We chose five representative points from the deviatoric stress-axial strain curve and presented a microcrack corresponding to each point, as shown in Figure 13.At point a, a small number of tensile cracks appeared at the ends of the prefabricated fissures due to stress concentration.The cracks had a wing-shaped distribution, and the specimen without heat treatment exhibited more tensile cracks than the specimen after being subjected to 750 • C heat treatment.When the axial deviatoric stress reached point b, the number of micro tensile cracks on two sides of prefabricated fissures significantly increased and gradually expanded along the direction of the maximum principal stress, while a small number of micro shear cracks appeared, especially for a specimen that was not heat treated.From point b onward, the number of micro tensile cracks increased significantly, and the number of micro shear cracks increased gradually; both specimens exhibited this phenomenon.At point c, the axial deviatoric stress reached its peak, and macroscopic shear cracks formed and connected the two pre-existing cracks.At the outboard of the pre-existing cracks, the shear cracks extended to the upper and lower ends of the specimen.With the increase in axial deviatoric stress to d and e, the shape of macroscopic cracks was basically fixed, but the width of macroscopic cracks gradually increased with the increase microcracks.It should be noted that the specimen without thermal treatment exhibited more microcracks than the specimens after being subjected to 750 • C thermal treatment, and these extra cracks were mainly shear cracks.

Crack Evolution
There were small differences in damage patterns of specimens due to consistent cation, number, and inclination angle of prefabricated fissures.Hence, we selected thermal damage numerical specimens at 25°C and 750 °C to illustrate the evolution deviatoric stress and the number of cracks versus axial strain.We chose five representat points from the deviatoric stress-axial strain curve and presented a microcrack cor sponding to each point, as shown in Figure 13.At point a, a small number of tensile cra appeared at the ends of the prefabricated fissures due to stress concentration.The cra had a wing-shaped distribution, and the specimen without heat treatment exhibited m tensile cracks than the specimen after being subjected to 750 °C heat treatment.When axial deviatoric stress reached point b, the number of micro tensile cracks on two side prefabricated fissures significantly increased and gradually expanded along the direct of the maximum principal stress, while a small number of micro shear cracks appear especially for a specimen that was not heat treated.From point b onward, the numbe micro tensile cracks increased significantly, and the number of micro shear cracks creased gradually; both specimens exhibited this phenomenon.At point c, the axial de atoric stress reached its peak, and macroscopic shear cracks formed and connected two pre-existing cracks.At the outboard of the pre-existing cracks, the shear cracks tended to the upper and lower ends of the specimen.With the increase in axial deviato stress to d and e, the shape of macroscopic cracks was basically fixed, but the width macroscopic cracks gradually increased with the increase microcracks.It should be no

Conclusions
We proposed a numerical model to simulate the mechanical behavior of thermal damaged rock, based on a series of physical experiments and 2D granular mechanics simulations on pre-cracked sandstone specimens subjected to different temperature thermal treatments.By comparing the results of the numerical and physical experiments, we were able to define the mechanisms of crack evolution at fine scale.The main conclusions of this study are as follows: (1) Based on the SEM of the rock samples after thermal treatment at different temperatures, the thermal damage value D T , obtained by extracting the thermal crack area, can be used as an indicator of the degree of thermal damage of the sandstone.(2) The thermal damage numerical model, established by replacing the flat-joint model with the smooth-joint model according to the thermal damage value D T , can simulate the mechanical behavior and failure patterns of properly.(3) The maximum strength rises with increasing pre-treatment temperature from 25 to 750 • C, but subsequently declines as the pre-treatment temperature approaches 1000 • C. Therefore, the temperature at which the strength begins to decrease signifi-cantly is 750 • C, indicating a critical threshold.The elastic modulus E 1 decreases with the increasing thermal treatment temperature.(4) The evolution of microcracks starts with tensile cracks distributed in a wing shape at the ends of the prefabricated fissures to two sides of the prefabricated fissure, accompanied by a small number of shear cracks, and finally, develops into macroscopic fracture that extend through the entire specimen along the prefabricated fissures.The microcracks of the specimen without thermal treatment more than those of the specimens after 750 • C thermal treatment, and these extra cracks are mainly shear cracks.

Figure 1 .
Figure 1.The sampling location and the microstructure of fissured sandstone sample.The fi microstructure of sandstone adapted with permission from Ref. [14].2018, Rock Mech.Rock

Figure 1 .
Figure 1.The sampling location and the microstructure of fissured sandstone sample.The fi microstructure of sandstone adapted with permission from Ref. [14].2018, Rock Mech.Rock

Figure 3 .
Figure 3. High-temperature furnace and heating scheme.(a) High temperature heating furnace.(b) Heating scheme of sandstone specimens.

Figure 3 .
Figure 3. High-temperature furnace and heating scheme.(a) High temperature heating furnace.(b) Heating scheme of sandstone specimens.

Figure 4 .
Figure 4. SEM observation results and crack extraction results of different thermally damaged red sandstones (part of the SEM results were drawn from literature [33]).

Figure 5 .Figure 4 .
Figure 5.Effect of thermal treatment on E and DT.

Figure 4 .
Figure 4. SEM observation results and crack extraction results of different thermally damaged red sandstones (part of the SEM results were drawn from literature [33]).

Figure 5 .Figure 5 .
Figure 5.Effect of thermal treatment on E and DT.

Figure 6 .
Figure 6.Notation used to define joint and smooth-joint contact [36].

Figure 6 .
Figure 6.Notation used to define joint and smooth-joint contact [36].

Figure 8 .
Figure 8.The effects of thermal cracks on the deviatoric stress in simulations.

Figure 8 .
Figure 8.The effects of thermal cracks on the deviatoric stress in simulations.

Sustainability 2023 , 19 Figure 9 .
Figure 9. Comparisons of the stress-strain curves in tests and simulations of the fissured red sandstones after different temperature thermal treatment.

Figure 9 .
Figure 9. Comparisons of the stress-strain curves in tests and simulations of the fissured red sandstones after different temperature thermal treatment.

Figure 9 .
Figure 9. Comparisons of the stress-strain curves in tests and simulations of the fissured red sandstones after different temperature thermal treatment.

Figure 10 .Figure 10 .
Figure 10.Comparisons of the elastic modulus and peak deviatoric stress against heat temperature in experiments and simulations.(a) The effects of the temperature on elastic modulus.(b) The effects of the temperature on peak deviatoric stress.

Figure 11 .
Figure 11.Comparisons of the failure pattern for the fissured red sandstones with different therma temperature treatments.

Figure 11 .
Figure 11.Comparisons of the failure pattern for the fissured red sandstones with different thermal temperature treatments.

Sustainability 2023 ,
15,  x FOR PEER REVIEW 14 o prefabricated fissures, the displacement direction shifted from being parallel to the fiss to vertical orientation, exhibiting a closure of fissure.

Figure 12 .
Figure 12.Displacement of particles of the fissured red sandstones with different thermal temp ture treatments.

Figure 12 .Figure 13 .
Figure 12.Displacement of particles of the fissured red sandstones with different thermal temperature treatments.Sustainability 2023, 15, x FOR PEER REVIEW 15 of 18

Table 1 .
Basic physical parameters and experimental design scheme of fissured red sandstones.

Table 1 .
Basic physical parameters and experimental design scheme of fissured red sandstones.

Table 2 .
Modeling parameters of Flat-joint model in PFC2D simulation for the specimens

Table 2 .
Modeling parameters of Flat-joint model in PFC2D simulation for the specimens.

Table 3 .
Modeling parameters of the smooth-joint model in PFC2D simulation for the specimens.

Table 3 .
Modeling parameters of the smooth-joint model in PFC2D simulation for the specimens.