Microdamage study of granite under thermomechanical coupling based on the particle flow code

Abstract

The thermomechanical coupling of rocks refers to the interaction between the mechanical and thermodynamic behaviors of rocks induced by temperature changes. The study of this coupling interaction is essential for understanding the mechanical and thermodynamic properties of the surrounding rocks in underground engineering. In this study, an improved temperature-dependent linear parallel bond model is introduced under the framework of a particle flow simulation. A series of numerical thermomechanical coupling tests are then conducted to calibrate the micro-parameters of the proposed model by considering the mechanical behavior of the rock under different thermomechanical loadings. Good agreement between the numerical results and experimental data are obtained, particularly in terms of the compression, tension, and elastic responses of granite. With this improved model, the thermodynamic response and underlying cracking behavior of a deep-buried tunnel under different thermal loading conditions are investigated and discussed in detail.

Graphical abstract

Keywords

thermomechanical coupling effect / granite / improved linear parallel bond model / thermal property / particle flow code

1 Introduction

The thermomechanical coupling interaction of rocks is a common problem encountered in underground engineering, particularly for granite in nuclear engineering projects. During engineering operations, the extra heat produced by high-level radioactive waste is continuously transmitted into the surrounding rock through engineering barriers such as waste tanks and buffer materials. As a result, various types of thermal damage and some microcracks may appear in the surrounding rock around the disposal repository, further degrading the thermal and mechanical properties of the surrounding rocks [1]. Therefore, understanding the thermodynamic properties of the surrounding rock is crucial for determining the volume, layout, and cost of disposal repositories, optimizing the spacing of disposal units, and ensuring the safety of engineering projects [2].

Therefore, a large number of studies have been carried out to investigate the deformation and rupture of surrounding rock induced by thermomechanical coupling action [3–6]. Here, we consider granite as an example, as granite is often encountered in deep-buried facilities and is generally considered a preferential geo-material for high-level radioactive waste. Fan et al. [7] first investigated the local fracturing of granite due to temperature changes by analyzing its internal microstructural evolution using computed tomography and 3D image reconstruction. The results indicated that the density of granite decreased sharply between 400 and 800 °C. Chen et al. [8] completed a series of unconfined compressive strength (UCS) tests on granite after heating and found that the elastic response sharply declined as the temperature increased from 400 to 800 °C. Sun et al. [9] obtained similar results and further noted that the failure of granite comprised four stages. Moreover, the crack width and density of granite changed abruptly above 400 °C. Using scanning electron microscopy, X-ray, and acoustic emission techniques, Xu et al. [10] further explored the microcrack propagation process of granite via UCS tests with different temperature loadings and found that the failure pattern exhibited a clear transition from brittle to flexible. Moreover, Yin et al. [11] and Yang et al. [12] also observed that a critical temperature exists at which the physical and mechanical behavior of granite changes significantly. More precisely, when the temperature is sufficiently high, a reduction in the chemical strength of minerals, such as the α−β quartz phase transition, leads to failure of the granite [11,13–16]. Therefore, most previous experimental studies have shown that the type, distribution, and content of minerals in granite directly affects its thermomechanical responses [17–19] in terms of its elastic properties, compression, tension, and Poisson’s ratio [20–23].

Numerical methods have made significant progress in the last decade. Different types of numerical thermomechanical coupling algorithms have been proposed and used to study the thermodynamic behavior of granite [24–27]. Jang and Yang [28] simulated the granite heating process using a thermal elastic finite element program and compared the results with those of indoor tests. Schrank et al. [29] proposed a granite model based on an implicit Lagrangian finite element code to capture the amplification of rock pores caused by cracks during the heating process. By adopting a combined approach of rock failure process analysis and digital image processing, Yu et al. [30] developed a numerical model incorporating meso-structure characteristics that could successfully reproduce the local cracking related to thermal stress. Wang and Konietzky [31] used a continuum approach to investigate the effects of different heating speeds on granite. Other numerical methods, including meshed [32–35] and mesh-free methods [36–38], have also been used to capture the cracking initiation and propagation behavior in the failure process of samples; the results can provide a deeper understanding of local cracking in rocks under different loading conditions [39,40].

The particle flow code (PFC), which was first proposed by Cundall and Strack [41], is a discrete element method that is widely used [42,43] to simulate the mechanical behavior of cohesiveless sand–soil media [44] and analyze the cohesiveness of rocks [45–47]. In terms of thermodynamic coupling, PFC exhibits a remarkable advantage in modeling the heat-induced initiation and propagation of cracks in rocks. In this method [48], the continuous rock is regarded as an assembly of discrete particles bonded by interfaces. During thermomechanical coupling, these particles are considered as thermal reservoirs and the bonds are simultaneously considered as thermal pipes to transfer heat. The thermodynamic process of heat flowing outward from a heat source can thus be simulated. Wanne and Young [49] simulated the crack propagation of granite under heating conditions using a bonded-particle model and compared the results with those of acoustic emission tests. Using a coupling model, Xia et al. [50] considered the variation in the elastic strength and modulus with temperature. Zhao [14] simulated the generation of microscopic and macroscopic cracks in granite under heating and illustrated the mechanical mechanism of granite cracking under different temperature conditions using particle flow simulations.

However, rocks such as granite have often been simulated as isotropic materials in previous studies, thus neglecting the internal microstructural features. In experiments, the internal micro-inclusions of granite, such as plagioclase, feldspar, quartz, and mica, have been shown to have an obvious effect on the thermodynamic behavior. Therefore, this structural characteristic should be considered during modeling of the thermomechanical coupling. A simple thermomechanical coupling model originally embedded in the PFC has frequently been used in simulations; this model defines the thermal expansion of the contacts as a linear function of the strain. Different thermal loadings thus result in different nonlinear thermal strain responses. Therefore, this issue remains controversial.

A numerical sample generation approach is first proposed in this study to describe the internal inclusions of granite by defining different types of contacts. Then, by including a thermal damage function at the contacts, an improved thermal-dependent linear parallel bond model is suggested and programmed based on the PFC2D algorithm. The micro-parameters involved in the model are calibrated using laboratory test data obtained for Beishan granite from China. The thermodynamic properties of the granite during heating are investigated. The proposed model is then further applied to study the thermodynamic response and potential cracking behavior of a typical deep-buried tunnel under different thermal loading conditions.

2 Methodology

In this study, the linear parallel bond model (LPBM) [51] in PFC2D is used; this model can effectively reproduce the mechanical behavior of cementation related to rock minerals [52,53]. For the thermomechanical coupling, a thermal expansion function is incorporated into the model to describe the expansion of particles as well as that at the contacts. In addition, owing to oxidation or deterioration with increasing temperature [6,50,54], the mechanical and thermal properties of minerals will change during the heating process; thus, an improved LPBM based on temperature is proposed and used to establish a thermomechanically coupled granite sample.

For the LPBM, to describe the local deformation and failure behavior at the contacts, as shown in Fig.1, stiffness parameters (k_n and k_s) are first introduced into the contact laws to define the linear relationships between the contact forces and relative displacements, respectively, with the normal force and displacement (F_n and u_n, respectively), as well as the shear force and displacement (ΔF_s and Δu_s, respectively) [55,56]. It should be noted that the normal stiffness, k_n, is the secant modulus, whereas the shear stiffness, k_s, is tangential. The corresponding mathematical equations for the contact laws in the LPBM are given as follows:

(1)F_{\text{n}} =k_{\text{n}}{u}_{\text{n}} ,

(2)\mathrm{\Delta }{F}_{\text{s}}=-k_{\text{s}} \mathrm{\Delta }{u}_{\text{s}}.

To model the thermomechanical coupling behavior, the micromechanical parameters in terms of stiffness and strength in the linear contact bond portion of the LPBM are temperature-dependent. As presented in Fig.2, pieces 1 and 2, which can be a facet and a particle or two particles, represent the two ends of the contact. In the thermal calculation, these two pieces serve as thermal sources that can transfer heat to each other through the thermal pipe.

During the temperature change, ΔT, the radius of each particle and bond length of the contact are expressed as follows:

(3)\mathrm{\Delta }R\text{=}\alpha R\mathrm{\Delta }T ,

(4)\mathrm{\Delta }L\text{=}\overline{\alpha }L\mathrm{\Delta }T,

where ΔR and R refer to the radius change after heating and radius before heating of the particle, respectively, α is the coefficient of the linear thermal expansion (CLTE) of the particle, ΔL is the length change of the contact after heating, \overline{\alpha } is the CLTE of the contact, and L is the length of the contact before heating, which is a heat-transfer medium in the thermal calculations.

Therefore, the normal contact force change, \mathrm{\Delta } {\overline{F}}_{\text{n}}, of the bond between two particles at each time step can be obtained as follows:

(5)\mathrm{\Delta }{\overline{F}}_{\text{n}} =2{\overline{k}}_{\text{n}} R_{\text{c}} \mathrm{\Delta }L,

where {\overline{k}}_{\text{n}} is the normal stiffness of the bond, and R_{\text{c}} is the width of the contact.

The heat transmission equation for a continuum is given in the thermal calculations as follows:

(6)\frac{\mathrm{\partial }{\mathit{q}}_i}{\mathrm{\partial }{x}_i}+ q_{\text{v}} =\rho C_{\text{v}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t},

(7){\mathit{q}}_i=-{\mathit{k}}_{ij}\frac{\mathrm{\partial }T}{\mathrm{\partial }{x}_j},

where {\mathit{q}}_i and q_{\mathrm{v}} denote the heat flux vector and volumetric heat source intensity, respectively, and \rho, C_{\text{v}}, and {\mathit{k}}_{ij} denote the mass density, specific heat, and thermal conductivity tensors, respectively.

The heat transfer equation for two particles can be expressed as follows:

(8)-\sum _{p=1}^N{Q}_p+Q_{\text{v}} =m{C}_{\text{v}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t},

(9)Q_p=-\frac{\mathrm{\Delta }{T}_p}{{\eta }_p{L}_p},

(10){\eta }_p=\frac{{\eta }_{p1}+{\eta }_{p2}}{L_p},

where Q_p and Q_v denote the thermal energy in the pth pipe flowing out of the heat source and the heat source intensity, respectively; N is the number of interparticle contacts, i.e., the number of thermal pipes; m is the thermal mass; ΔT_p is the temperature difference between the two ends of the pth pipe; η_p is the thermal resistance at the unit length of the pth pipe; η_p1 and η_p2 are the resistances of pieces 1 and 2, respectively; and L_p is the length of the pth pipe.

3 Calibration and verification of the thermomechanical coupling model

3.1 Calibration of the mechanical properties

A numerical sample of Beishan granite for thermomechanical testing is established based on experimental data for Beishan granite from China [57]. The granite samples in the experimental data are composed of plagioclase (33%–44%), feldspar (22%–26%), quartz (29%–38%), and mica (3%–5%).

Image recognition and processing [58,59] is a common numerical remodeling method. As shown in Fig.3, the numerical sample is generated in PFC2D with the appropriate mineral content using the random composition method for the minerals [60]. The numerical sample is composed of plagioclase (37.2%), feldspar (25.8%), quartz (32.9%), and mica (4.1%). The parameters of the numerical sample are summarized in Tab.1 and Tab.2.

The microscopic parameters are obtained by trial and error such that the stress–strain curves of the numerical sample match the experimental data [57] under different confining pressures to establish a reasonable granite sample in the numerical simulation. The strengths of the different minerals are defined in accordance with the strength relationship (quartz > plagioclase ≈ feldspar > mica) [61]. In this study, the modulus parameters of plagioclase, feldspar, and mica are set to 0.5, 0.45, and 0.25 times that of feldspar, respectively. The strength parameters of plagioclase, feldspar, and mica are set to 0.7, 0.8, and 0.4 times that of feldspar, respectively. The stress–strain curves of the numerical simulation and experimental data at confining pressures of 5, 10, and 15 MPa are shown in Fig.4. For the boundary conditions of the numerical tests, the contacts between the side walls and particles are defined as linear with an effective modulus of 0.52 GPa and a ratio of the normal stiffness to shear stiffness of 20.0. The contacts between the top wall and particles are defined as linear with an effective modulus of 10.4 GPa and ratio of the normal stiffness to shear stiffness of 0.0. The step strain rate is 4.02 × 10⁻⁹ step⁻¹, which is appropriate for quasi-static loading conditions [62].

The numerical simulation results and experimental data are in good agreement, and the parameters of the LPBM are listed in Tab.3. Each contact is named based on the mineral with the lower strength between the two ends of the contact.

3.2 Calibration of the thermal properties

In the thermomechanical coupling analysis of the PFC, the timescales associated with the thermal and mechanical processes are different, and the ratio of the thermal-to-mechanical timescales is extremely large; therefore, taking several thermal time steps without the mechanical state of the sample is feasible if the system approaches thermal equilibrium [48]. In the heating process in this study, the thermal-only mode is applied for a specific time with a temperature increment of 1 °C to ensure that the system approaches thermal equilibrium. The mechanical-only mode is then applied until the system approaches mechanical equilibrium. As shown in Fig.5, temperature is applied to the four walls, and heat is transferred from the walls to the particles and then from those particles to other particles through the thermal pipes. The temperature of the walls is initially set to 25 °C (room temperature), and the sample is solved to approach mechanical equilibrium. In the first stage of the heating process, the thermal-only mode is used, and the sample is cycled for 12 s after the wall temperature is set to 26 °C. The mechanical-only mode is then used, and the sample is run for a number of simulation cycles until the sample reaches a mechanical equilibrium state. This process is repeated until the temperature of the sample reaches 1000 °C. During this process, the heating rate is 5 °C/min.

The microscopic thermal parameters of the thermal resistance can be calibrated by measuring the macroscopic thermal conductivity of the numerical sample. In experiments, steady-state measurements are frequently employed to determine the thermal conductivity. The material is heated until the thermal equilibrium point is reached. When the heat transfer and heat dissipation rates are equal, the material is considered to have reached equilibrium. The thermal conductivity can be calculated using the heat flux density passing through the specimen, thickness, and temperature differential between the two sides in accordance with the Fourier 1D steady-state heat conduction model [63]. The heat conductivity measurements are shown in Fig.6(a).

In the numerical method, a constant thermal power of 100 W is applied to the particles at the bottom, and the sample is solved such that the heat from the bottom particles is transferred to the particles at the top of the sample. Two narrow areas are selected on the heated sample: the first area with an average temperature of T₁ is very close to the heat source, and the second area with an average temperature of T₂ is close to the top of the sample. The thermal conductivity, k, is calculated using Fourier’s law, which can be expressed as follows:

(11)k=\frac{Qh}{(T_1-T_2)d},

where Q is the thermal power of 100 W applied to the bottom particles, h is the distance between the centers of the two areas, and d is the width of the sample. As shown in Fig.6(b), the thermal conductivity curve converges to a relatively stable value during the heating process, and this value represents the thermal conductivity at a certain temperature. The thermal conductivity curves are shown in Fig.6(c). Under these conditions, the thermal conductivity obtained by the numerical simulation can be matched with the thermal conductivity obtained experimentally at various temperatures. The resulting mesoscopic thermal conductivity of the ith particle is negatively correlated with the temperature and satisfies the following exponential distribution:

(12){\eta }_i(T)=0.68{\eta }_{i0}\times {1.00045}^{\frac{T}{1.04}}.

According to Refs. [64,65], the CLTE, α, approximately satisfies an exponential distribution. The CLTE affects the strength of the rock after heating by influencing the degree of thermal damage and number of thermal cracks. The CLTE of the ith particle is obtained as follows:

(13){\alpha }_i(T )=(0.4{\text{e}}^{\frac{T}{350}}+ 0.36){\alpha }_{i0}.

The specific heat, C_v, increases linearly at low temperatures [66], and mutation occurs with the phase transition [67]. Therefore, the specific heat capacity of the ith particle is fitted as follows [68]:

(14)C_{\text{v}i}(T)= \{\begin{array}{rlr}& 0.16T+1011,& (T⩽573{}^{\circ }\text{C})\\ & a+b{T}^{-2}+c{T}^{-3}+d{T}^{-0.5}+e{T}^{-1},& (T>573{}^{\circ }\text{C})\end{array}

where a, b, c, d, and e are set to −228.24, 3.3801 × 10⁸, −3.4738 × 10¹⁰, 1.11 × 10⁵, and −2.38 × 10⁶, respectively [69].

Several experiments have demonstrated that quartz material undergoes an α–β phase transition at 573 °C [54,70,71], which causes a change in the properties of granite-containing quartz materials at high temperatures. Piecewise functions are used in the expressions for the elastic modulus, E, and strength, ψ (including the tensile strength, σ, and shear strength, τ), of the cith contact to reflect this change, which can be expressed as follows:

(15)E_i(T)=\{\begin{array}{l}(0.00059T+0.985)E_{i0},(T⩽573{}^{\circ }\text{C})\\ (0.000028{\text{e}}^{\frac{T}{85}}+1.27)E_{i0},(T>573{}^{\circ }\text{C})\end{array}

(16){\psi }_{ci}(T)=\{\begin{array}{l}{\psi }_{ci0},(T⩽573{}^{\circ }\text{C})\\ (0.000098{\text{e}}^{\frac{T}{114}}+0.992){\psi }_{ci0},(T>573{}^{\circ }\text{C})\end{array}

where i represents different types of minerals; α_i0, η_i0, E_i0, and ψ_ci0 are the CLTE, thermal resistance, effective modulus, and strength parameter of contact i at 25 °C, respectively.

The preheating process is conducted from 25 to 1000 °C using the model described above. The variation in deformation inside the granite sample leads to thermal cracking when the temperature increases because of the rapid production of cracks and the different CLTE values of the minerals. The distribution of microcracks formed by thermal cracking is relatively uniform. Owing to the expansion of the mineral particles after heating, which causes thermal stress, a total of 1757 thermal cracks form. The number of thermal cracks continues to increase with increasing temperature. Fig.7(a) shows the distribution of the thermally produced microcracks in the sample at various temperatures. The simulation results show that the growth rate of the crack number and order of crack generation differ for various minerals. Quartz experiences minimal damage during heating owing to its high strength. Cracking occurs earliest in mica, initiating at 429 °C. The plagioclase begins to crack at 536 °C, and the feldspar begins to crack at 717 °C. The crack production rate is slow during the heating process because of the lower mica content. When the sample is heated to 791 °C, the cracks induced by the failure of plagioclase dominate. The number of cracks in the plagioclase increase rapidly. More than 60% of the total cracks are induced above 800 °C.

As shown in Fig.7(b), the expansion of the numerical sample increases with the temperature and reaches approximately 3% at 1000 °C. The expansion trend with temperature agrees well with the results of the laboratory testing [64]. The thermal conductivity curves of the Beishan granite are obtained from laboratory experiments in Refs. [72,73]. The microscopic thermal parameters are obtained by trial and error to ensure that the thermal conductivity curves of the numerical sample match the experimental data. The thermal parameters of the numerical samples are presented in Tab.4.

3.3 Verification of the stress–strain curves

Following the compression tests using the parameters in Tab.3 and Tab.4, the stress–strain curves at different temperatures are obtained, as shown in Fig.8. The stress–strain curves generated by the numerical simulations (lines) and the results of the experiments (scatters) are in good agreement, as shown in Fig.8(a). The strength and elastic modulus of the granite sample decrease with temperature. However, when constant parameters are used, as shown in Fig.8(b), the compressive strength and elastic modulus obtained numerically are higher than the experimental values, particularly at high temperatures. This result is attributed to the fact that the thermal and mechanical properties of minerals vary with temperature, and the complex chemical reactions involving minerals occur at temperatures exceeding 573 °C.

Fig.9 shows the microcrack distributions at various temperature peak stress points; the results indicate that the macrocracks at the peak change with temperature. At lower temperatures (approximately 200 °C), the microcrack distribution exhibits a common brittle failure mechanism in Beishan granite, manifested as an “X”-shaped failure pattern [14]. At higher temperatures (> 400 °C), the main failure mode is shear failure with a steep angle of approximately 60°–70°. As a result of thermally generated microcracks and thermal stress, granite samples may become more ductile as the temperature increases.

A tensile test is conducted on the heated Beishan granite using the same sample as that used for the compression test. Fig.10 shows the resulting stress–strain curves and microcrack distributions. According to the results of the direct tensile testing, the temperature affects the tensile strength of Beishan granite. No obvious fluctuation is observed in the curves at low temperatures (< 400 °C), whereas the curves fluctuate sharply at high temperatures (> 400 °C). The curves fluctuate less with the generation of fewer thermal cracks during preheating and then sharply increase with an increase in the thermal cracks. This is because adjacent cracks open and propagate under tension when thermal cracks are present. The connection and transfixion with adjacent cracks result in a sudden decrease in the strength. In contrast, the cracks close under pressure in the compression tests. Therefore, fluctuations do not occur in the compression tests. The loading cracks are dominated by tensile cracks. When the temperature is lower than 600 °C, the microcracks at the peak tensile strength generated by loading appear as a horizontal line in the failure area, which is a typical tensile failure mode. More thermal cracks than loading cracks occur at temperatures above 600 °C, and because the thermal cracks prevent the loading cracks from growing, the loading cracks are dispersed randomly with intermittent shear cracks.

4 Influence of the temperature field of a tunnel on the excavation damage

4.1 Size effect

The size effect is investigated to ensure that the thermal and mechanical characteristics of the small-scale sample and large-scale model are consistent. The particle distribution remains the same, but the height and width of the sample, gaps between particles, and radii of the particles are increased, as shown in Fig.11(a). Six samples are thus obtained at magnifications of 1, 2, 5, 100, 500, and 2000×. The stress–strain curves obtained in the compression tests and the microcrack distributions are shown in Fig.11(b) and Fig.11(c). Under the condition that the particle distribution remains unchanged, the physical and mechanical characteristics of the sample remain unchanged regardless of the magnification (error within 5%). The main failure mode is shear failure with a steep angle of approximately 60°–70° for all samples. Therefore, a magnified sample can be used to perform the tunnel heating simulation.

4.2 Model tunnel construction

It is necessary to investigate the thermomechanical coupling effect while considering the failure mechanisms of tunnel excavations at high temperatures. To confirm the accuracy of the proposed model, it is necessary to analyze the impact of the temperature field on tunnel excavation, as illustrated in Fig.12.

Based on the size effect tests, 18 samples magnified at 500× are combined to form a large model, as shown in Fig.12. The length and width of the model are 150 m. A rectangular area of 125 m × 140 m is removed from the model as the initial tunnel model. The obtained model system is guaranteed to be consistent with the small-scale sample in the compression and tensile tests. Thus, it is possible to investigate the impact of the temperature field on tunnel excavation using the parameters listed in Tab.3 and Tab.4. A flexible membrane boundary [74] is used to apply geostress to the model and to simulate the geostress gradient at different depths. The burial depth of the model is assumed to be 80 m. The stress conditions are shown in Fig.12. The lateral stress increases linearly with depth. The bottom of the model is fixed. A circular area with a radius of 10 m is excavated. The total number of discs is approximately 70000.

To investigate the failure mechanism of tunnel excavation under thermomechanical coupling, the temperature of the heating area is increased from 25 to 1000 °C at a rate of 5 °C/min. The temperature of the boundary is fixed at 25 °C. During this process, the temperature of the surrounding rock around the tunnel gradually increases, and heat is gradually transferred to the inside of the rock mass. The crack propagation in the surrounding rock is analyzed through the simulated heating process.

4.3 Crack propagation behavior of a tunnel under temperature changes

The tunnel is heated using the method described above. To examine the impact of temperature-dependent thermal factors, Fig.13 shows the development of the number of cracks with temperature for the tunnel model.

When the temperature is below 200 °C, the tunnel model contains almost no cracks, regardless of whether the temperature-dependent parameters are considered. When the temperature is increased to 200 °C or higher, scattered cracks appear throughout the tunnel. The frequency of cracks increases dramatically between 300 and 400 °C, and a through-fissure develops in the middle of the model. The number of cracks increases linearly between 400 and 1000 °C. When the temperature reaches 1000 °C, heat can be transmitted to a distance of approximately 45 m from the tunnel. The deformation caused by the cracks prevents heat transfer in a standard circular form. The temperature gradient in the undamaged area is relatively uniform, whereas the temperature of the crack area is discontinuous. A crack can cause a significant temperature differential between nearby rock masses, and this temperature difference-induced deformation continues to increasing, causing the fissure to widen. When the temperature around the tunnel drops to 25 °C, the temperature near the through-crack decreases rapidly, returning to close to the normal temperature. More cracks are observed in the heating area around the hole, and the rock at the boundary of the tunnel is mostly detached. The heat transfer process decelerates, and the surrounding rock becomes more dispersed because of the cracks. However, points of contact are still observed between the rock masses, and the maximum temperature remains at 400 °C and decreases evenly with the depth.

As shown in Fig.13(a), for the tunnel modeled using the conventional thermomechanical model, cracks start to occur at 325 °C, and the number of cracks grows evenly above 350 °C. The number of cracks at 1000 °C is 3361. Fig.13(b) shows that the cracks are induced at 300 °C and begin to grow evenly after 375 °C when the temperature-dependent thermomechanical model is adopted. The number of cracks is 3100 at 1000 °C. Compared with the model without temperature-dependent thermal parameters, the rock mass begins to break earlier, and fewer cracks are present in the initial damage stage. The elastic modulus decreases with temperature; as a result, the rock mass is soft at high temperatures, and less deformation can be caused by thermal expansion. Therefore, the crack growth rate is reduced significantly at high temperatures.

5 Discussion

The heating rate has a significant impact on the thermal properties of granite [75]. The Beishan granite model is heated at different hating rates from room temperature (25 °C) to 1000 °C to elucidate the growth and distribution behavior of cracks after heating under different heating rates. The heating rates are 1, 5, 10, 30, and 60 °C/min. The numbers of tensile, shear, and total cracks are obtained, as shown in Fig.14. At lower (< 30 °C/min) or higher (> 30 °C/min) heating rates, the number of cracks is greater after heating to 1000 °C. When the heating gradient is 30 °C/min, the final number of cracks is the minimum, only 77% of that observed at a heating rate of 1 °C/min. The numerical results reveal that the heating gradient affects the number of thermal damage cracks, which in turn affect the mechanical characteristics of granite.

6 Conclusions

Based on the thermal calculation framework in PFC2D, an improved LPBM that characterizes the mechanical and thermal characteristics of rocks with varying temperature is proposed, and a temperature-dependent thermomechanical model is established. The mechanical properties of granite under thermomechanical coupling are analyzed based on laboratory tests of Beishan granite. The thermodynamic response of crack propagation in deep-buried tunnels is investigated. The main conclusions are as follows.

1) The results of thermal experimental tests of granite indicate that the mechanical behavior of granite changes with temperature owing to thermal damage. This thermal damage can be divided into two parts: the first part is caused by thermal expansion effects, and the second part is caused by the reduction of the mineral mechanical and thermal properties with temperature.

2) To simulate the damage induced by the thermal effect in granite at high temperatures, a temperature-dependent LPBM is implemented using the PFC2D algorithm. The particles and contacts are assigned the relevant thermal and mechanical properties. Temperature decreases the compressive strength, tensile strength, and elastic modulus of the Beishan granite sample. The tensile curves fluctuate at high temperatures. The simulation results closely reflect the laboratory test results. Thus, the models are capable of recreating the microscopic and macroscopic cracks, stress–strain behavior, and failure mechanisms of Beishan granite.

3) In the Beishan granite tunnel model, heat is continuously transferred to the interior of the rock layer under heating of the tunnel. Random microcracks appear around the hole, and a through-crack is observed in the middle of the model.

4) A number of sensitivity investigations of the heating rate are carried out. The number of thermal damage cracks in granite is significantly influenced by the heating gradient. The specimen exhibits the fewest thermal cracks at a heating rate of 30 °C/min.

References

[1] Chaki S, Takarli M, Agbodjan W P. Influence of thermal damage on physical properties of a granite rock: Porosity, permeability and ultrasonic wave evolutions. Construction & Building Materials, 2008, 22(7): 1456–1461

[2] Cai Y Y, Luo C H, Yu J, Zhang L M. Experimental study on mechanical properties of thermal-damage granite rock under triaxial unloading confining pressure. Journal of Geotechnical Engineering, 2015, 37: 1173–1180

[3] Zhao F, Shi Z M, Sun Q. Fracture mechanics behavior of jointed granite exposed to high temperatures. Rock Mechanics and Rock Engineering, 2021, 54(5): 2183–2196

[4] Li D Y, Su X L, Gao F H, Liu Z D. Experimental studies on physical and mechanical behaviors of heated rocks with pre-fabricated hole exposed to different cooling rates. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2022, 8(4): 125

[5] Huang Y H, Yang S Q, Tian W L, Zhao J, Ma D, Zhang C S. Physical and mechanical behavior of granite containing pre-existing holes after high temperature treatment. Archives of Civil and Mechanical Engineering, 2017, 17(4): 912–925

[6] Sun Q, Zhang W Q, Zhu Y M, Huang Z. Effect of high temperatures on the thermal properties of granite. Rock Mechanics and Rock Engineering, 2019, 52(8): 2691–2699

[7] Fan L F, Gao J W, Wu Z J, Yang S Q, Ma G W. An investigation of thermal effects on micro-properties of granite by X-ray CT technique. Applied Thermal Engineering, 2018, 140: 505–519

[8] Chen Y L, Ni J, Shao W, Azzam R. Experimental study on the influence of temperature on the mechanical properties of granite under uni-axial compression and fatigue loading. International Journal of Rock Mechanics and Mining Sciences, 2012, 56: 62–66

[9] Sun Q, Zhang W Q, Xue L, Zhang Z Z, Su T M. Thermal damage pattern and thresholds of granite. Environmental Earth Sciences, 2015, 74(3): 2341–2349

[10] XuX LKang Z XJiMGeW XChenJ. Research of microcosmic mechanism of brittle-plastic transition for granite under high temperature. Procedia Earth and Planetary Science. 2009, 1(1): 432–437

[11] Yin T B, Li X B, Cao W Z, Xia K W. Effects of thermal treatment on tensile strength of Laurentian granite using Brazilian test. Rock Mechanics and Rock Engineering, 2015, 48(6): 2213–2223

[12] Yang S Q, Ranjith P G, Jing H W, Tian W L, Ju Y. An experimental investigation on thermal damage and failure mechanical behavior of granite after exposure to different high temperature treatments. Geothermics, 2017, 65: 180–197

[13] Glover P W J, Baud P, Darot M, Meredith P G, Boon S A, Leravalec M, Zoussi S, Reuschlé T. α/β phase transition in quartz monitored using acoustic emissions. Geophysical Journal International, 1995, 120(3): 775–782

[14] Zhao Z H. Thermal influence on mechanical properties of granite: A microcracking perspective. Rock Mechanics and Rock Engineering, 2016, 49(3): 747–762

[15] Van der Molen I. The shift of the α−β transition temperature of quartz associated with the thermal expansion of granite at high pressure. Tectonophysics, 1981, 73(4): 323–342

[16] Branlund J M, Hofmeister A M. Thermal diffusivity of quartz to 1000 °C: Effects of impurities and the α−β phase transition. Physics and Chemistry of Minerals, 2007, 34(8): 581–595

[17] Staněk M, Geraud Y. Granite microporosity changes due to fracturing and alteration: Secondary mineral phases as proxies for porosity and permeability estimation. Solid Earth, 2019, 10(1): 251–274

[18] Ghasemi S, Khamehchiyan M, Taheri A, Nikudel M R, Zalooli A. Crack evolution in damage stress thresholds in different minerals of granite rock. Rock Mechanics and Rock Engineering, 2020, 53(3): 1163–1178

[19] ZhangYZhao Y S. Thermal cracking meso-characteristic of LuHui granite. In: Proceedings of the International Conference on Mechanical Engineering and Green Manufacturing (MEGM) 2010. Xiangtan: Mechanical Engineering and Green Manufacturing, 2010

[20] ZhaoY SMeng Q RKangT HZhangNXiB P. Micro-CT experimental technology and meso-investigation on thermal fracturing characteristics of granite. Journal of Geotechnical Engineering, 2008, 27: 28−34 (in Chinese)

[21] Chen S W, Yang C H, Wang G B. Evolution of thermal damage and permeability of Beishan granite. Applied Thermal Engineering, 2017, 110: 1533–1542

[22] Lin W R. Permanent strain of thermal expansion and thermally induced microcracking in Inada granite. Journal of Geophysical Research. Solid Earth, 2002, 107(B10): 107

[23] Wang H F, Bonner B P, Carlson S R, Kowallis B J, Heard H C. Thermal stress cracking in granite. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1989, 26(5): 234

[24] Hu X D, Song X Z, Liu Y, Li G S, Shen Z H, Lyu Z H, Liu Q L. Lowest required surface temperature for thermal spallation in granite and sandstone specimens: Experiments and simulations. Rock Mechanics and Rock Engineering, 2019, 52(6): 1689–1703

[25] Liu H, Zhang K, Liu T, Cao H, Wang Y. Experimental and numerical investigations on tensile mechanical properties and fracture mechanism of granite after cyclic thermal shock. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2022, 8(1): 18

[26] Wang F, Konietzky H, Frühwirt T, Dai Y J. Laboratory testing and numerical simulation of properties and thermal-induced cracking of Eibenstock granite at elevated temperatures. Acta Geotechnica, 2020, 15(8): 2259–2275

[27] Zhang Y M, Mang H A. Global cracking elements: A novel tool for Galerkin-based approaches simulating quasi-brittle fracture. International Journal for Numerical Methods in Engineering, 2020, 121(11): 2462–2480

[28] Jang M, Yang H S. Experiment and numerical analysis for thermo-mechanical behavior of granite block. Geosystem Engineering, 2005, 8(1): 9–14

[29] SchrankCFusseis FKarrechARevetsSRegenauer Lieb KJieL. Thermal cracking of Westerly granite: From physical to numerical experiment. In: Proceedings of the EGU General Assembly Conference Abstracts. Vienna: EGU General Assembly, 2010

[30] Yu Q L, Ranjith P G, Liu H Y, Yang T H, Tang S B, Tang C A, Yang S Q. A mesostructure-based damage model for thermal cracking analysis and application in granite at elevated temperatures. Rock Mechanics and Rock Engineering, 2015, 48(6): 2263–2282

[31] Wang F, Konietzky H. Thermal damage evolution of granite under slow and high-speed heating conditions. Computers and Geotechnics, 2020, 123: 103590

[32] Zhang Y M, Zhuang X Y. Cracking elements: A self-propagating strong discontinuity embedded approach for quasi-brittle fracture. Finite Elements in Analysis and Design, 2018, 144: 84–100

[33] Zhang Y M, Zhuang X Y. Cracking elements method for dynamic brittle fracture. Theoretical and Applied Fracture Mechanics, 2019, 102: 1–9

[34] Zhang Y M, Huang J U, Yuan Y, Mang H A. Cracking elements method with a dissipation-based arc-length approach. Finite Elements in Analysis and Design, 2021, 195: 103573

[35] Zhang Y M, Gao Z R, Li Y Y, Zhuang X Y. On the crack opening and energy dissipation in a continuum based disconnected crack model. Finite Elements in Analysis and Design, 2020, 170: 103333

[36] Rabczuk T, Zi G, Bordas S, Nguyen Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[37] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[38] Yu S Y, Ren X H, Zhang J X, Wang H J, Sun Z H. An improved form of smoothed particle hydrodynamics method for crack propagation simulation applied in rock mechanics. International Journal of Mining Science and Technology, 2021, 31(3): 421–428

[39] Zhang Y M, Yang X Q, Wang X Y, Zhuang X Y. A micropolar peridynamic model with non-uniform horizon for static damage of solids considering different nonlocal enhancements. Theoretical and Applied Fracture Mechanics, 2021, 113: 102930

[40] Zhang Y M, Lackner R, Zeiml M, Mang H A. Strong discontinuity embedded approach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations. Computer Methods in Applied Mechanics and Engineering, 2015, 287: 335–366

[41] Cundall P A, Strack O D L. A discrete numerical model for granular assembilies. Geotechnique, 1979, 29(1): 47–65

[42] Zhang Y L, Shao J F, Liu Z B, Shi C, de Saxce G. Effects of confining pressure and loading path on deformation and strength of cohesive granular materials: A three-dimensional DEM analysis. Acta Geotechnica, 2019, 14(2): 443–460

[43] Zhang Y L, Shao J F, de Saxce G, Shi C, Liu Z B. Study of deformation and failure in an anisotropic rock with a three-dimensional discrete element model. International Journal of Rock Mechanics and Mining Sciences, 2019, 120: 17–28

[44] Cundall P A, Hart R D. Numerical modelling of discontinua. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1993, 30: 101–113

[45] Zhang Y L, Liu Z B, Han B, Zhu S, Zhang X. Numerical study of hydraulic fracture propagation in inherently laminated rocks accounting for bedding plane properties. Journal of Petroleum Science Engineering, 2022, 210: 109798

[46] Zhang Y L, Zhang Y P, Han B, Zhang X, Jia Y. Parameter studies on hydraulic fracturing in brittle rocks based on a modified hydromechanical coupling model. Energies, 2022, 15(7): 2687

[47] Zhang Y L, Shao J F, Zhu S, Liu Z B, Shi C. Effect of rock anisotropy on initiation and propagation of fractures due to fluid pressurization. Acta Geotechnica, 2022, 18(4): 2039–2058

[48] ItascaConsulting Group. User’s Manual: PFC2D-particle Flow Code in 2 Dimensions. Minneapolis: Itasca Consulting Group, 2004

[49] Wanne T S, Young R P. Bonded-particle modeling of thermally fractured granite. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(5): 789–799

[50] Xia M, Zhao C B, Hobbs B E. Particle simulation of thermally-induced rock damage with consideration of temperature-dependent elastic modulus and strength. Computers and Geotechnics, 2014, 55: 461–473

[51] Potyondy D O, Cundall P A. A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329–1364

[52] Zhang Y L, Liu Z B, Shi C, Shao J F. Three-dimensional reconstruction of block shape irregularity and its effects on block impacts using an energy-based approach. Rock Mechanics and Rock Engineering, 2018, 51(4): 1173–1191

[53] Zhang Y P, Shi C, Zhang Y L, Yang J X, Chen X. Numerical analysis of the brittle-ductile transition of deeply buried marble using a discrete approach. Computational Particle Mechanics, 2021, 8(4): 893–904

[54] Sun Q, Chen S E, Gao Q, Zhang W Q, Geng J S, Zhang Y L. Analyses of the factors influencing sandstone thermal conductivity. Acta Geodynamica et Geomaterialia, 2017, 14: 173–180

[55] Zhang Y L, Shao J F, Liu Z B, Shi C. An improved hydromechanical model for particle flow simulation of fractures in fluid-saturated rocks. International Journal of Rock Mechanics and Mining Sciences, 2021, 147: 104870

[56] ZhangY LHan BZhangXJiaYZhuC. Study of interactions between induced and natural fracture effects on hydraulic fracture propagation using a discrete approach. Lithosphere, 2021: 5810181

[57] HuS HZhang GZhangMJiangX LChenY F. Deformation characteristics tests and damage mechanics analysis of Beishan granite after thermal treatment. Rock and Soil Mechanics, 2016, 37: 3427–3436 (in Chinese)

[58] Zhang Y M, Gao Z R, Wang X Y, Liu Q. Image representations of numerical simulations for training neural networks. Computer Modeling in Engineering & Sciences, 2023, 134(2): 821–833

[59] Zhang Y M, Gao Z R, Wang X Y, Liu Q. Predicting the pore-pressure and temperature of fire-loaded concrete by a hybrid neural network. International Journal of Computational Methods, 2021, 19(8): 2142011

[60] Shi C, Yang W K, Yang J X, Chen X. Calibration of micro-scaled mechanical parameters of granite based on a bonded-particle model with 2D particle flow code. Granular Matter, 2019, 21(2): 38

[61] HanZZhangL ZhouJYuan GWangP. Uniaxial compression test and numerical studies of grain size effect on mechanical properties of granite. Journal of Engineering Geology, 2019, 27: 497−504 (in Chinese)

[62] Zhang X P, Wong L N Y. Choosing a proper loading rate for bonded-particle model of intact rock. International Journal of Fracture, 2014, 189(2): 163–179

[63] BonettoFLebowitz J LBelletL R. Fourier’s Law: A Challenge to Theorists, in Mathematical Physics 2000. London: Imperial College, 2000

[64] Plevova E, Vaculikova L, Kozusnikova A, Ritz M, Martynkova G S. Thermal expansion behaviour of granites. Journal of Thermal Analysis and Calorimetry, 2016, 123(2): 1555–1561

[65] Zhu Z N, Yang S Q, Wang R, Tian H, Jiang G S, Dou B. Effects of high temperature on the linear thermal expansion coefficient of Nanan granite. Acta Geodaetica et Geophysica, 2022, 57(2): 231–243

[66] Xiang P, Xu H C, Ji H G, Li Q, Wang H. Thermal property of granite in deep strata and its effect on thermal zone of surrounding rock. Shock and Vibration, 2022, 2022: 1–9

[67] Yurtseven H, Desticioglu M. Critical behaviour of the heat capacity near the α−β phase transition in quartz. High-Temperature Materials and Processes, 2013, 32(2): 189–194

[68] Saxena S K. Earth mineralogical model: Gibbs free energy minimization computation in the system MgO-FeO-SiO₂. Geochimica et Cosmochimica Acta, 1996, 60(13): 2379–2395

[69] Miao S Q, Li H P, Chen G. Temperature dependence of thermal diffusivity, specific heat capacity, and thermal conductivity for several types of rocks. Journal of Thermal Analysis and Calorimetry, 2014, 115(2): 1057–1063

[70] Pereira A H A, Miyaji D Y, Cabrelon M D, Medeiros J, Rodrigues J A. A study about the contribution of the α−β phase transition of quartz to thermal cycle damage of a refractory used in fluidized catalytic cracking units. Cerâmica, 2014, 60(355): 449–456

[71] Sun Q, Lu C, Cao L W, Li W C, Geng J S, Zhang W Q. Thermal properties of sandstone after treatment at high temperature. International Journal of Rock Mechanics and Mining Sciences, 2016, 85: 60–66

[72] Zhao X G, Zhao Z, Guo Z, Cai M, Li X, Li P F, Chen L, Wang J. Influence of thermal treatment on the thermal conductivity of Beishan granite. Rock Mechanics and Rock Engineering, 2018, 51(7): 2055–2074

[73] Zhao X G, Xu H R, Zhao Z, Guo Z, Cai M, Wang J. Thermal conductivity of thermally damaged Beishan granite under uniaxial compression. International Journal of Rock Mechanics and Mining Sciences, 2019, 115: 121–136

[74] Wang Y H, Leung S C. A particulate-scale investigation of cemented sand behavior. Canadian Geotechnical Journal, 2008, 45(1): 29–44

[75] Wang F, Konietzky H, Fruhwirt T, Li Y W, Dai Y J. Impact of cooling on fracturing process of granite after high-speed heating. International Journal of Rock Mechanics and Mining Sciences, 2020, 125: 104155

RIGHTS & PERMISSIONS

Higher Education Press