As the depth of underground engineering increases (e.g., resource extraction, space development, and energy utilization), high geo-temperature has become one of the inevitable geological disasters (Breede et al. 2013; Scott et al. 2015; Vilarrasa et al. 2017). Especially for the Sichuan-Tibet Railway, there are about 15 tunnels along the line encounter high geo-temperature problems (Xue et al. 2020; Tian et al. 2021). Such as the Sangzhuling tunnel that has been completed, the measured maximum temperature of the rock during construction exceeds 89 °C, and the maximum buried depth reaches 1500 m (Yan et al. 2018); for the Layue tunnel under construction, the maximum borehole temperature measurement data reaches 93.5 °C, and the maximum buried depth reaches 2080 m (Wang et al. 2021). In addition, the high geo-temperature tunnel of Sichuan-Tibet Railway also faces ultra-low ambient temperature conditions, and the lowest temperature can reach below −30 °C (Hu et al. 2021). In summary, rocks are generally subjected to the combined effects of high temperature difference and high crustal stress during the excavation. Therefore, studying the mechanical behavior and meso-structure of rock under high temperature difference and axial loading is important for revealing the failure mechanism of rock and preventing the occurrence of engineering disasters.

Up to now, many scholars have studied the physical and mechanical properties of rock materials after high temperature effect (e.g., density, porosity, P-wave velocity, thermal conductivity, Poisson's ratio, uniaxial compressive strength (UCS), elastic modulus) (David et al. 1999; Glamheden and Lindblom 2002; Funatsu et al. 2004; Chaki et al. 2008; Zhang et al. 2012; Sengun 2013; Kim et al. 2014; Wang et al. 2016; Yüksek 2019; Guo et al. 2020a). González-Gómez et al. (2015) studied the effect of temperature on the physical properties of limestone with different porosity. The results showed that the limestone with lower porosity showed brittle behavior and had higher compressive strength and elastic modulus. Zhao et al. (2018) treated Beishan granite at high temperature and studied the effect of treatment temperature on the thermal conductivity of granite. The results show that the thermal conductivity of granite decreases nonlinearly with the increase of temperature. Jiang et al. (2018) shows that the P-wave velocity and S-wave velocity of granite almost decrease linearly with the increase of temperature, and the permeability of granite almost increases exponentially with the increase of temperature. Wang et al. (2020) studied the effects of treatment temperatures of 400, 600 800 and 1000 °C on the physical, mechanical and thermal properties of Eibenstock granite. The results show that the treatment temperature can significantly affect basic parameters such as mineral composition, P-wave velocity, UCS, Young's modulus and thermal conductivity, and the thermal expansion coefficient increases sharply at about 573 °C and 870 °C.

Microscopic damage mechanism and macroscopic failure characteristics of rock materials after high temperature have also been widely studied (Li and Sh. 1993; Hale and Shakoor 2003; Yavuz et al. 2006; Kim and Kemeny 2009; Enayatpour and Patzek 2013; Zeng et al. 2019; Guo et al. 2020b). Yu et al. (2014) used digital image modeling to study the thermal cracking process caused by different expansibility of granite at high temperature. The results show that the elastic modulus of rock is only affected by thermal damage, and the UCS is affected not only by high temperature thermal stress, but also by cumulative thermal damage. Zhang et al. (2019) mainly studied the influence of temperature and prefabricated cracks on the failure mode of granite, and explained the formation mechanism of different failure modes of granite under the combined action of temperature effect and cracks. Shen et al. (2019) prepared antifreeze by adjusting the concentration of calcium chloride solution, and conducted cooling shock (20 °C, 0 °C and − 30 °C) on granite samples with target temperatures of 150, 350, 550 and 750 °C. The results showed that with the decrease of refrigerant temperature, the local distribution characteristics of granite macro cracking became more obvious.

In recent years, some scholars have carried out relevant research on the mechanical behavior of rocks under thermo-mechanical (TM) coupling effect (Li et al. 2010; Xia et al. 2014; Wang et al. 2015; Xue et al. 2016; Tao et al. 2021). Zuo et al. (2012) conducted a TM coupling failure test on Beishan granite with double precast notches. The results show that the mineral composition and particle size of granite seriously affect its failure mechanism and fracture toughness. The granite is mainly intragrain fracture at low temperature, and intragrain and trans granular coupling fracture mechanism at high temperature. Yin et al. (2016) respectively carried out static mechanical tests of granite after pre-heating treatment and real-time high temperature treatment. The results show that there are significant differences in elastic modulus, damage characteristics, peak stress and failure mode under the two conditions. Liu et al. (2020) according to the variation law of mechanical properties of rock mass with temperature, a TM coupling simulation method considering the heterogeneity of rock mass at high temperature is put forward. The change of mechanical properties of rock mass after heating is the main reason for the change and failure of underground gasification strata in coalfield. Yin et al. (2021) studied the triaxial stress of granite at high temperature (100–400 °C), revealing the difference of thermal and mechanical characteristics between deep granite and shallow granite at high temperature.

To sum up, although many scholars have done a lot of research on thermal treatment and TM coupling treatment of rock. However, most studies focus on: (1) The influence of temperature above 200 °C on the rock; (2) The rapid cooling of rock is concentrated in water cooling (20 °C) and liquid nitrogen cooling (−196 °C), which is inconsistent with the working conditions of surrounding rock of high geo-temperature tunnel affected by ultra-low temperature (−30 °C). In addition, there are few studies on the failure mechanism of rock under the coupling action of high temperature difference and axial loading. Therefore, the purpose of this study is to investigate the mechanical behavior and failure mechanism of sandstone under high temperature difference (−30 to 150 °C) and axial loading. The cooling rate is used as a variable to simulate the effect of different intensity of high temperature difference on rock. First, the TM coupling experimental system is established to study the mechanical properties and failure modes. Then the present study establishes a sandstone TM coupling simulation system to further analyze the influence of cooling rate on the mechanical behavior and failure mechanism of sandstone. It includes the change of mechanical properties, the crack propagation path and macroscopic failure mode in the process of uniaxial compression, the damage characteristics of cooling rate on meso-structure and so on. The research results are intended to provide reference for the construction of high geo-temperature tunnel and other deep engineering encountering high geo-temperature.

2.1 Sample preparation

In the present study, the sandstone widely distributed in high geo-temperature tunnel was tested (Xie et al. 2022). All samples were taken from the same block of each rock, and selected from density determinations and longitudinal wave velocity measurements. The samples were prepared as cylindrical, with diameter of 50 mm and height of 100 mm. Both the samples’ size and processing accuracy satisfied the International Society of Rock Mechanics (ISRM) requirements. The X-ray diffraction (XRD) show that the sample mainly consists of quartz, albite and clay minerals (kaolinite, chlorite and illite), with a content of 28.3%, 39.5% and 26.8%, as well as a small amount of calcite and micro plagioclase, which are 2.3% and 3.1%, respectively. In addition, the natural density of the sample is 1.88 ± 0.3 g/cm³.

2.2 Experimental apparatus

AE monitoring system, triaxial test system and high & low temperature test system were combined in this study (Fig. 1a). The PCI-2 system of PAC acoustic emission instrument was adopted (Fig. 1a(1–1)). The triaxial test system is a 2000 kN electric servo triaxial tester (Fig. 1a(2)), the measuring accuracy is better than 1%, and the loading rate is 0.01–20 kN/s. The temperature test system is a GD8/10 high & low temperature test chamber (Fig. 1a(3)). The test chamber adopts a split structure, includes test box (Fig. 1b(3–1)) and refrigeration unit (Fig. 1b(3–2)), which are connected by ventilation pipe. The controllable temperature range is −80 °C–150 °C, and the effective volume of box is 0.01 m³.

Experimental platform a Schematic of experimental system b Main test equipment

Firstly, sensors such as extender, temperature sensor and AE sensor are connected from the reserved pipe (on the back of the test box), and relate to the rock sample respectively (Fig. 1b(1–2)). Then the test box (Fig. 1b(3–1)) is placed on the triaxial tester (Fig. 1b(2–2)), and the pressure bar exerted by triaxial tester is closely connected with the test box to ensure that all the axial pressure exerted on rock sample. After connecting the ventilation pipe, adjust the air flow temperature to achieve the purpose of heating and cooling the test box. After that, the AE damage location of sandstone is carried out by the time difference location method, and four AE probes are arranged symmetrically. Vaseline is used as a coupling agent between the AE sensor and the sample surface to achieve good acoustic coupling. In addition, in order to reduce the effect of background noise, the gain value of the pre-amplifier, and the threshold for AE detection were set to 40 dB and 50 dB, respectively.

2.3 Experimental process

This study sets the axial stress loading level as about 60% of the UCS, and the specific value is 22 MPa. Chen et al. (2022) shows that after the action of this stress level, the damage value of sandstone can reach more than 24% with porosity as the index. The target temperature of thermal treatment is 150 °C and −30 °C, respectively. The axial loading and thermal treatment are realized in the test box (Fig. 1b(3–1)). For better results, two samples are repeated under each condition. All the samples were dried before test. The specific experimental process is as follows:

1. (1)

   The axial stress is applied to the sample, the loading rate is 2 MPa/min, the loading value is 22 MPa, and the axial stress is stable for 5–10 min (Li et al. 2018).

2. (2)

   Heat the sample to 150 °C at the heating rate of 5 °C/min, the heating rate is the same as that of the previous experiments on sandstone (Ding et al. 2016; Zhu et al. 2016). After reaching 150 °C, keep the temperature for 2 h.

3. (3)

   The temperature of the sample was cooled to −30 °C at a cooling rate of 5 °C/min. After reaching − 30 °C, keep the temperature for 2 h. AE positioning was carried out in this process.

4. (4)

   The sample is subjected to uniaxial compression experiment until failure at −30 °C, and the loading rate is 0.1 mm/min. At the same time, the force–displacement data are recorded once per 500 ms.

5. (5)

   Change the cooling rate in step (3) and repeat steps (1)–(4).

Due to the limitation of the test equipment, the cooling rate is only set at 5 °C/min and 10 °C/min in the laboratory test, Fig. 2 shows the experimental flow chart.

Experimental flow chart

Figure 3a presents the axial stress–strain curves of the samples, where room temperature represents the samples without TM coupling treatment, and the other two groups represent the samples treated at different cooling rates. It is clear from Fig. 3a that the curve can be roughly divided into four stages: compaction, elasticity, yield and failure. Compared with the curves of room temperature and 5 °C/min, the TM coupling effect has significantly affected the mechanical properties of sandstone. The UCS decreases by about 16% and the peak strain also decreased. The UCS difference between 5 and 10 °C/min curves is small, with a decrease of only about 3%, and the peak strain still tends to decrease.

Experimental results a Stress–strain curves b AE event distribution c Failure mode

During the cooling stage, the AE event distribution with a height of 20 mm in the middle of the sample was selected. Figure 3b shows the results of four samples. Under the cooling rate of 5 °C/min, the AE event distribution of the two samples is relatively random, and there is no obvious rule in spatial distribution. For the two samples treated at a cooling rate of 10 °C/min, the distribution density of AE events in the exterior region of the sample increases slightly, and the total number of AE events increased. It indicating that the cooling rate of 10 °C/min intensifies the damage of sandstone.

Figure 3c shows the macroscopic failure mode of sandstone samples. The samples show a typical oblique shear failure at room temperature, and the integrity after failure is high. Except for a main crack that runs through the sample, there are only a few secondary cracks derived from the main crack. After TM coupling treatment, the mode changed significantly. The number of macroscopic main cracks increases, and the secondary cracks derived from the main cracks also increase obviously. The degree of fragmentation is higher than that at room temperature. In the samples with a cooling rate of 10 °C/min, it is also found that the surface falls off (the place indicated by the circle), and the integrity of the samples after failure is the lowest.

4.1 Numerical model

A general particle-flow model simulates the mechanical behavior of a system made up of a collection of arbitrarily shaped particles. The model is composed of distinct particles that displace independent of one another and interact at pair-wise contacts. Newton’s laws of motion provide the fundamental relationship between particle motion and the forces causing that motion. The force system may be in static equilibrium, or it may be such as to cause the particles to flow.

The rock model can be composed of many rigid particles and bond models. Several rigid particles are combined to represent minerals, and the physical properties (e.g., density, stiffness) of minerals are reflected on particles. Rigid particles cannot be destroyed, and the overall strength of rock is controlled by the bond model between particles. The stiffness and strength of the bond model need to be adjusted according to different mineral components, so that the model can better reproduce many characteristics of rock behavior, including elasticity, fracture, expansion, damage accumulation and so on. It is worth noting that when the bond between two particles is broken, microcracks will be formed in the middle of the fractured bond through the built-in function. The direction of the microcrack is perpendicular to the fractured bond, and the length is equal to the minimum diameter of the particles at both ends of the fractured bond.

Based on the Grain-based model (GBM) in PFC2D, an equivalent model of heterogeneous sandstone is established, as shown in Fig. 4a. The height of the model is 100 mm, the width is 50 mm. There are 17,393 particles, the radius is 0.15–0.40 mm, the particle size obeys uniform distribution. Compared with the Particle-based model (PBM), the advantages of the GBM model are as follows: (1) Heterogeneity in the mechanical properties of mineral grains can be reproduced; (2) The geometric heterogeneity induced by grain boundaries can be revealed; and (3) Grains can be crushed under certain stress levels (Li et al. 2019). In the GBM model, mineral crystals are composed of multiple particles. The linear parallel bond model (PB model) is used inside the mineral crystals, and the smooth joint model (SJ model) is adopted between the mineral crystals in order to better reflect the sliding behavior, as shown in Fig. 4b. The process of modeling is as follows: (1) Firstly, the mineral composition and content of the model are determined according to the XRD results of sandstone, and the PBM model is generated according to basic physical and mechanical parameters of each mineral; (2) Convert PBM model to GBM model according to the rblock module built in PFC; (3) Before the assignment of microscopic parameters, the confining pressure is applied by the servo principle to simulate the underground occurrence environment of rock; (4) Assign a set of microscopic parameters to the bond model, unloading confining pressure to simulate the process of taking the rock from ground; (5) Repeat step (4) to complete the calibration of microscopic parameters. Then carry out the simulation steps such as axial loading and thermal treatment.

Numerical model a Mineral distribution of the model b Mechanical bond model c Thermal resistance distribution of thermal pipe

The simulation steps are the same as those in laboratory test. For axial stress loading, the stress change may cause large disturbance in the model. Therefore, based on servo principle, a multi-gradient graded loading model is established, as presented in Eq. (1). Where σ₀ represents the target stress of the current loading stage, m represents the current number of cycles, n represents the total number of cycles, β represents the attenuation speed, σ represents the current stress of the model, σ_t represents the total target stress (that is 22 MPa). In the process of multi-stage loading, the next stage loading is carried out when the unbalanced force ratio of the model reaches 0.00001.

The thermal calculation begins after the axial stress loading is completed, and the heating rate is still 5 °C/min. When the heating boundary of the model reaches 150 °C, the temperature shall be kept constant until the temperature field is uniform. Different from laboratory tests, the model simulates the effects of various cooling rates on sandstone, including 5, 10, 20, 50, 100 and 180 °C/min. Figure 5 shows the change curve of the model temperature boundary.

Temperature boundary variation curve of numerical model

4.2 Parameter calibration

In order to make the mechanical properties of the model consistent with the sandstone samples, the micro parameters need to be calibrated according to static experiments. The "trial and error method" proposed by Peng et al. (2018) is used for calibration. In the present study, the stress–strain curve of the sandstone samples is obtained by the uniaxial compression experiment. After that, the curve obtained by the model is consistent with the experiment results (Fig. 6a). Table 1 shows the microscopic parameters after calibration. Due to the model does not simulate the closure behavior of original fractures, the stress–strain curve obtained does not reflect the compaction stage. According to the post-peak stage of stress–strain curve, the model shows obvious brittle fracture characteristics, which is consistent with the experiment results.

Calibration results of numerical model a Room temperature b 5 °C/min c 10 °C/min

PB model plays an irreplaceable role in PFC thermal module, it is the only bond model with thermal expansion phenomenon. The change of temperature leads to the change of particle radius and bond force, which produces thermal strain and thermal stress respectively.

where ΔR represents the change of particle radius, α represents the linear thermal expansion coefficient of the particle, R represents the particle radius, and ΔT represents the temperature change of the particle. \(\Delta \bar{F}_{{\text{n}}}\) is the normal component of the force vector carried by the bond, \(\bar{k}_{{\text{n}}}\) is the bond normal stiffness, A is the area of the bond cross-section, ᾱ is the expansion coefficient of the bond material (taken equal to the average value of the expansion coefficients of the particles at the two ends of the pipe associated with the bond), \(\overline{L}\) is the bond length (taken equal to the distance between the centroids of the two particles at the ends of the pipe associated with the bond), and ΔT is the temperature increment of the bond (taken equal to the average temperature change of the two particles at the ends of the pipe associated with the bond).

The thermal parameters of different minerals need to be assigned respectively, including thermal conductivity, linear thermal expansion coefficient, specific heat capacity and so on. It is worth noting that the difference of thermal conductivity between different minerals has not been considered in most studies. However, the thermal conductivity is assigned respectively by Fish function in this study, which increases the accuracy of thermal calculation. Table 1 shows the main thermal parameters of the model. Figure 4c shows the thermal resistance distribution of each thermal pipe. The specific formula is as follows:

where η denote the thermal resistance of the thermal pipe, k denote the thermal conductivity, n is the porosity within the measurement sphere, V^(b) is the volume of ball (b), l^(p) is the length of thermal pipe (p) associated with balls in the measurement sphere, N_b and N_p is the number of balls and thermal pipes, η₁ and η₂ denote the thermal resistance of piece 1 and 2, respectively.

As a thermal boundary, wall is only a one-dimensional line model and does not participate in mechanical calculation in PFC2D, so it can only be assigned a specific temperature. In order to better reproduce the impact of continuous change of ambient temperature on rock samples, the present study realizes the continuous change of wall temperature according to callback function in Fish. For reduce the influence of unbalanced stress on the model, each step is calculated until the unbalanced force ratio reaches 0.0001. Finally, the thermal calculation of the model is also verified (Fig. 6b, c). The uniaxial compression simulation is carried out for the sandstone treated with different cooling rates. The results show that the numerical simulation results are consistent with the experiment results, and the error is kept in a reasonable range.

5.1 Mechanical properties of sandstone at different cooling rates

The uniaxial compression simulation of cooled sandstone was carried out to analyze the changes of its mechanical properties. The loading rate was 0.1 mm/min, and acoustic emission counts (AE counts) monitoring was realized by PFC2D. Figure 7 shows the changes of axial stress, AE counts and AE cumulative counts with axial strain under different cooling rates, and shows the corresponding macroscopic failure mode (the type of microcracks represented by different colors is the same as Fig. 12). AE counts have obvious quiet period and active period. In the elastic deformation stage and before, the sample accumulates a lot of energy, almost no AE events occur, and the AE counts is in a quiet period. The accumulated energy is released rapidly in the plastic deformation stage, the AE counts increases rapidly and begins to enter the active stage. Obviously, the AE cumulative counts curve is related to cooling rate. When the cooling rate is less than 20 °C/min, the early change of the curve is relatively smooth, and a sharp increase occurs near the peak stress. As the increase of cooling rate, the growth rate of the curve increases and the overall slope becomes larger, which is consistent with the distribution characteristics of AE counts, that is, more AE events were detected at the beginning of active period.

Evolutions of the axial stress, AE counts and cumulative AE counts of sandstone as axial strain

Figure 8 shows the changes of peak stress, peak strain and elastic modulus of sandstone under different cooling rates. All of them showed a decreasing trend with the increase of cooling rate. Compared with untreated sandstone, the peak stress of sandstone decreased by 15%, 23%, 25%, 37%, 41% and 54%, respectively, the peak strain decreased by 20%, 26%, 28%, 37%, 38% and 44% respectively, and the elastic modulus decreased by 5%, 7%, 9%, 10%, 15% and 21%, respectively. It can be seen that the peak stress and peak strain change obviously. However, the variation range of elastic modulus is relatively small. The elastic modulus at room temperature is 9.7 GPa, and the corresponding elastic modulus at 180 °C/min is 7.95 GPa.

Changes of mechanical properties

5.2 Uniaxial failure characteristics of sandstone at different cooling rates

The present study realizes the microcrack hotspot map which is used to characterize microcrack density by means of Fish function and Scalar module in PFC2D, as shown in Fig. 9. The microcrack hotspot map is composed of many square measuring units of equal size, all of which spread out sequentially and occupy the whole sample space. Each square measuring unit can retrieve the number of microcracks within its search radius. Finally, the size of visualization area is based on the number of microcracks it searches. The microcrack density in different regions can be shown qualitatively through microcrack hotspot map.

Calculation flow of microcrack hotspot map

The cooled sandstone shows a different crack propagation path from the untreated sandstone (namely room temperature) during uniaxial compression. Limited by the length of paper, this study only takes the sandstone with a cooling rate of 50 °C/min as the representative, and shows the crack propagation path between it and untreated sandstone in the form of microcrack hotspot map (Fig. 10). In room temperature, before the sandstone reaches UCS, the microcracks randomly occur in the sample, and the distribution is relatively uniform. Near the UCS, the connectivity of microcracks is strengthened, and macroscopic cracks tend to be formed uniformly. After failure, the macroscopic crack is relatively simple, and the overall fragmentation degree is low. Before the UCS of sandstone treated by TM coupling, the distribution of microcracks is uneven in space, and the density of microcracks in the exterior region of the sandstone is obviously higher. Near the UCS, the macroscopic crack is more likely to be formed from outside to inside. The macroscopic mode after failure is also different from room temperature, showing X-type conjugate failure, accompanied by the formation of obvious secondary cracks, and the overall fragmentation degree increases.

Microcrack hotspot map of crack propagation path under uniaxial compression simulation

According to the influence of cooling rate on macroscopic failure mode (refer to Fig. 7), it is found that the fragmentation tends to increase with the increase of cooling rate. However, the failure mode does not change much, there is X-type conjugate failure in most cases. To further analyze the effect of cooling rate on it, this study quantifies the type and number of microcracks according to the preset function. Figure 11 shows the proportion of various intragrain cracks and grain boundary cracks produced during uniaxial compression. The ratio of various crack types fluctuates within a certain range under different cooling rates, so the cooling rate does not have a significant effect on the proportion. The preliminary analysis of the cause of this phenomenon is related to the damage characteristics of meso-structure of sandstone in the cooling stage, which will be discussed in more detail in Sect. 5.3.

Percentage of microcracks under uniaxial compression simulation

5.3 Mesoscopic damage analysis during the cooling stage

The present study selects three critical moments to analysis the distribution of microcrack and the change of microcrack density during the cooling stage, namely Moment A, Moment B and Moment C, as shown in Fig. 12. In Moment A, microcracks mainly occur in the exterior region of the sample, and this phenomenon becomes more obvious with the increase of cooling rate. The occurrence of microcracks in this stage is mainly due to the temperature change in the exterior region, which leads to the thermal stress in this region exceeds its strength. In Moment B, except for the increasing number of microcracks in the exterior region, some randomly distributed microcracks begin to appear in the interior. However, according to the hotspot map, the distribution of microcracks is still mainly in the exterior region. Moment C means that the sample is cooled to −30 °C. It can be seen that the spatial distribution of microcrack tends to be uniform, and the regions with high microcrack density are randomly distributed while the cooling rate below 20 °C/min. As the increase of cooling rate, although the number of microcracks increases compared with the previous two moment, the damage in the exterior region is still dominant.

Microcrack distribution and crack density distribution of sandstone

The number of microcracks produced in the cooling stage is also quantified according to preset function, as shown in Fig. 13. Along the radial direction, the sample is divided into five regions equidistant, marked as I, II, III, IV and V, and the number of microcracks is counted by region. The number of microcracks in each region tends to be uniform while the cooling rate below 20 °C/min. As the increase of cooling rate, the number of microcracks shows a distribution characteristic of more exterior and less interior. The faster the cooling rate, the more significant the trend. For instance, when the cooling rate is 50 °C/min, the region with the largest number of microcracks is I, a total of 149, followed by 111 in V, and the region with the smallest number is III, only 65. When the cooling rate is 100 °C/min, the region with the largest number of microcracks is I, 324 in total, followed by 264 in V, and the region with the smallest number is III, only 113. Affected by the inhomogeneity of mineral distribution and thermal conduction in numerical model, the microcracks are not symmetrically distributed along the axial direction, and the number of I is always the largest. Finally, the damage of sandstone in the cooling stage is mainly concentrated in the exterior region, and the destroy of interior structure is not serious. This is one of the reasons why there is no significant change about failure morphology in uniaxial compression simulation.

Number of microcracks in different region

6.1 Effect of cooling rate on internal heat conduction of sandstone

Wu and Liu (2003) proposed that the heating mode of rock can be divided into two ways: uniform heating and sharp heating. Due to the low thermal conductivity of rock, the hysteresis effect of heat conduction in time and space will be occurred while heated rapidly. There will be a great difference in temperature distribution between the surface and interior of rock, which will lead to the stress concentration and difference of stress distribution (Tan et al. 2006). In addition, the temperature gradient produced by the temperature change rate of 5 °C/min is about 0.27 °C/mm. When the temperature gradient exceeds 0.3 °C/mm, the effect of thermal stress caused by the temperature gradient on microcracks cannot be ignored (Fredrich and Wong 1986; Avanthi Isaka et al. 2019). Thus, the change rate of rock ambient temperature will have a significant impact on its internal temperature conduction, which will further affect the distribution of thermal stress. When the thermal stress exceeds its strength, thermal damage will occur in the rock.

Figure 14 presents the temperature changes at different positions while the temperature boundary of sample decreases from 150 to − 30 °C. The position of the temperature probe is marked in Fig. 14a, and the temperature field distribution of the model is also given when the temperature boundary reaches −30 °C. As shown in Fig. 14a, the slope of the three curves tends to be the same, and the internal temperature change of the sample is relatively uniform. With the increase of cooling rate, the difference between the curve’s slope under each cooling rate gradually increases, and the temperature difference inside and outside the rock gradually increases. It indicates that the faster the cooling rate is, the more obvious the difference of temperature change rate at different positions is, that is, the temperature change rate gradually decreases from exterior region to interior region of the sandstone. This further shows that the temperature gradient in the exterior region is relatively larger.

Variation of temperature probe with loading step

Table 2 shows the immediate temperature of each probe when the temperature boundary reaches −30 °C. The internal maximum temperature (data monitored by probe 3) of rock is only −14.89 °C while the cooling rate is 5 °C/min. With the increase of cooling rate, the temperature difference inside and outside the rock gradually increases. The maximum temperature corresponding to 180 °C/min is 142.53 °C, and the temperature difference inside and outside reaches 172.53 °C. Due to the increase of temperature difference, the temperature gradient in sandstone cannot be ignored, and the thermal failure is the result of the joint action of temperature gradient and thermal stress mismatch of adjacent minerals. In addition, this study is assumed that the thermal conductivity and thermal expansion coefficient of minerals in sandstone are not affected by the cooling rate. This further shows that the excessive temperature gradient in the exterior region of the sandstone is the main reason for the damage concentration in this region (refer to Fig. 12).

6.2 Contribution of intragrain cracks and grain boundary cracks to sandstone damage during the cooling stage

Figure 15 shows the number of intragrain and grain boundary cracks produced in the cooling stage. In the process of cooling, intragrain failure is the main form, and the proportion of grain boundary failure is small. It is worth noting that when the cooling rate is less than 100 °C/min, the damage of clay minerals is the largest, followed by albite. With the increase of cooling rate, the number of microcracks in albite gradually exceeds that of clay minerals, and the damage of albite crystal is the most serious. It is clear that the severity of crystal damage is not only related to cooling rate, but also related to mineral composition, the most seriously damaged crystal types are not invariable. The causes may include: (1) There are many grain boundary cracks in the pre-treatment process, so the grain boundary damage is not serious in the cooling stage; (2) The cold load acts directly on the surface of the sample and has a large contact area with mineral crystal, which promotes the destruction of crystal structure; (3) The thermal crack will be affected by the interaction of axial stress, mineral distribution characteristics and so on in the process of gradual penetration, which may lead to different cumulative damage rates of minerals.

Number of microcracks in the cooling stage

In addition, Fig. 15 shows the change of total number of microcracks produced in the cooling stage. The curve increased exponentially with the increase of cooling rate. The total number of cracks between different cooling rates increased by 13%, 29%, 60%, 71% and 47%, respectively. It is found that 50 °C/min can be used as a demarcation point. When the cooling rate is less than 50 °C/min, the cumulative damage rate of sandstone is in a relatively slow growth stage, on the contrary, the cumulative damage rate is significantly accelerated.

6.3 Stress distribution of sandstone under axial pressure during rapid cooling

Kim et al. (2014) found that when the pre-heated rock is cooled, tensile stress is generated in the exterior region of the sample, and compressive stress is generated in the region near the center of the sample. Li et al. (2020) also describes this phenomenon, as shown in Fig. 16a. When the external region of the sample is under tension, it will promote the further expansion of the crack in this region. When the central region of the sample is under pressure, the crack in this region will be compacted. Additionally, the rock material has obvious brittle characteristics, which is more prone to failure under tensile stress. In the present study, the cooling process of pre-heated rock is also simulated by numerical model, and the stress distribution is shown in the form of force chain diagram (Fig. 16b), the results are consistent with the previous research conclusions. However, the premise of the above conclusion is that the boundary condition of rock is free and is not subject to axial loading. It is worth noting that the sandstone model in this study will be controlled by constant axial pressure during the cooling stage. This factor will change the stress distribution in sandstone, resulting in the stress distribution no longer appearing the phenomenon described in Figs. 16a, b. As shown in Fig. 16c, the stress distribution is compression along the axial stress loading direction, and the tensile stress is distributed along the radial direction or at a certain angle with the radial direction.

Stress distribution and crack rose diagram of pre-heated rock in cooling stage a Stress distribution without axial stress studied by Li et al. (2020) b Stress distribution without axial stress studied by the model in this study c Magnification of stress distribution under axial pressure d Inclination distribution of microcracks

The present study calculates the inclination distribution of cracks during the cooling stage (take 50 °C/min as an example, the inclination is obtained when the crack rotates counterclockwise along the horizontal line), as shown in Fig. 16d. The inclination of microcracks is mostly distributed around 90°. The number of microcracks with an inclination of 70°–90° accounts for 27% of the total number of cracks, and the number of cracks with an inclination of 90°–110° accounts for 24%. It is worth noting that most of the cracks in the cooling stage are tensile cracks, which indicates that the cracks in this stage are mostly caused by the fracture of tensile bond distributed along the radial direction. Therefore, the stress distribution characteristics affected by axial loading have an inevitable influence on its failure in the cooling stage.

In the present study, laboratory tests and numerical calculations have been carried out on sandstone. Firstly, the influence of cooling rate on mechanical properties and failure models was been studied using the TM coupling experimental system. Then, the sandstone TM coupling simulation system is established according to the GBM model in PFC, and the effects of cooling rate on macro-mechanical characteristics and meso-structure are further analyzed. The main conclusions are as follows:

1. (1)

   The peak stress, peak strain and elastic modulus decrease with an increase in cooling rate. The uniaxial crack propagation path also changes significantly. The exterior region failure occurs first, and macroscopic cracks gradually penetrate from outside to inside. The fragmentation degree increases with an increase in cooling rate.

2. (2)

   Sandstone is dominated by intragrain failure during the cooling stage, and the number of microcracks increases exponentially with an increase in cooling rate. When the cooling rate exceeds 20 °C/min, the cracks produced in the cooling stage are concentrated in the exterior region, mainly due to the fracture of radial bonds.

3. (3)

   The stress distribution along the axial direction is in the form of compressive stress, and most of the tensile stress is distributed along the radial direction in the cooling stage. Therefore, the cracks in the cooling stage result from axial loading, temperature gradient and thermal stress mismatch between adjacent minerals. Moreover, the excessive temperature gradient in the exterior region of the sandstone is the main reason for the damage concentration in this region.