Effect of Thermal Cracking on the Tensile Strength of Granite: Novel Insights into Numerical Simulation and Fractal Dimension

Abstract

This study investigates the effect of thermal cracking on the tensile strength of granite through a combination of experimental testing and numerical simulations. The primary objective is to understand how thermal stress, induced by heat treatment at various temperatures (25 °C to 600 °C), influences crack initiation, propagation, and tensile strength changes. The granite specimens were subjected to Brazilian splitting tests after heat treatment, and the load–displacement curves and tensile strength variations with heat treatment temperature were analyzed. A grain-based model (GBM) was developed to simulate the complex cracking behavior, incorporating the mineral compositions and thermal expansion properties of the granite. The fractal dimension of the cracks was quantified using the box-counting method, and the relationship between fractal dimension and tensile strength was discussed. The results show that the GBM can effectively simulate the microcracking behavior and tensile fracture properties of heat-treated granite, accounting for mineral composition and thermal expansion. Thermal cracks are mainly intergranular tensile cracks, which increase in number with higher temperatures, while under mechanical loading failure is primarily due to intragranular tensile cracks. Higher heat treatment temperatures lead to denser crack networks with greater fractal complexity, reducing tensile strength and creating more tortuous crack propagation paths.

1. Introduction

As modern engineering projects increasingly involve high-temperature environments, understanding the mechanical properties of granite, a common rock in underground engineering, under the influence of high temperatures has become a critical research focus. Thermal cracking in rocks has been widely studied due to its implications for the stability and safety of rock engineering structures, particularly in structures exposed to high temperatures, such as in geothermal energy extraction, nuclear waste repositories, and tunneling projects [1,2,3,4,5]. However, the thermal cracking behavior of granite remains challenging to predict accurately due to its polycrystalline structure, varying grain sizes, and mineral composition [6,7,8]. This complexity is further amplified when granite is exposed to extreme temperatures, leading to thermal expansion, phase transitions in minerals, and eventual mechanical degradation [9,10].

One of the key challenges is to understand how microstructural changes in granite under the effect of high temperatures, particularly thermal cracking, affect its macroscopic mechanical properties. Previous studies have primarily focused on understanding how heating affects the compressive and tensile strength of various rocks [11,12]. Experimental studies have shown that higher temperatures tend to weaken rock, resulting in microcrack formation and increased porosity, which subsequently lowers mechanical strength [13,14,15]. While thermal cracking is often cited as a major factor in tensile strength reduction, the specific interaction between different crack types and their collective influence on mechanical failure is underexplored. Moreover, the precise mechanisms through which thermal cracking influences tensile strength, particularly the role of inter- and intragranular cracks, are still poorly understood.

Numerical modeling has become crucial for simulating microcracking behavior in heat-treated granite. In particular, the grain-based model (GBM) has enabled more detailed simulations of how the complex grain structure of granite responds to thermal stress. Zhao [16] explained the temperature-dependent mechanisms affecting the mechanical properties of granite through a particle-based model, revealing that strength reduction is primarily due to rising thermal stresses and the formation of tensile microcracks. Hu et al. [17] investigated the impact of grain size heterogeneity and temperature on the mechanical and microcracking behavior of Lac du Bonnet granite using a GBM, finding that as temperature and grain size heterogeneity increase, mechanical properties decline, microcracking intensifies, and the controlling factor of failure shifts from thermal cracks to grain size heterogeneity. Guo et al. [18] used a thermo-mechanical coupled numerical model based on computed tomography and a grain-based model (CT-GBM) to analyze the high-temperature microcracking behavior of granite, revealing that thermal damage primarily occurs through tensile failure, with grain boundary cracks, especially between feldspar grains, becoming dominant at higher temperatures. Pan et al. [19,20] investigated the mode-I and mode-II fracture behavior of heat-treated granite using notched semicircular bend (SCB) specimens and cracked straight-through Brazilian disc (CSTBD) specimens. They studied the evolution of thermal cracks in the heating and cooling phases by the GBM numerical method and found that thermal cracks can change the direction and path of the propagation of load-induced cracks. However, the fractal nature of crack propagation in heat-treated granite has not been investigated, leaving a gap in quantifying the complexity of crack patterns and how they relate to mechanical properties such as tensile strength.

This study aims to provide a more detailed understanding of the effect of thermal cracking on the tensile strength of granite through a combination of experimental and numerical approaches, using both qualitative and quantitative analysis. First, granite specimens were heat-treated at various temperatures (25 °C to 600 °C) and subjected to Brazilian splitting tests to evaluate their tensile strength. Subsequently, a numerical model considering mineral composition and thermal effects was developed using the GBM method. The model was validated based on the experimental results, which simulated the initiation and propagation of thermally induced and mechanically induced cracks, and characterized the evolution of cracks in granite specimens subjected to varying heat treatment temperatures. Finally, the fractal dimension of the cracks was quantified using the box-counting method, and the relationship between fractal dimension and tensile strength was discussed.

2. Experimental Materials and Methods

The granite blocks used in this study were collected from the batholith of Macheng City, Hubei Province, China. The mineralogical characteristics of the granite specimen under the polarizing microscope are shown in Figure 1. The main mineral components are plagioclase, K-feldspar, quartz, and biotite, and detailed mineralogical information can be obtained from previous studies [11,19].

Granite blocks were machined into disc specimens with a diameter of 50 mm and a thickness of 25 mm, and then heat-treated at different temperatures using a muffle furnace. A temperature range of 25 °C to 600 °C was chosen, which is very similar to the thermal environment of granite at shallow to moderate depths and would reflect the effects of mineral phase transformation or decomposition above about 573 °C. This is especially relevant for underground engineering and geothermal energy applications, where rock is often exposed to similar temperature conditions. The heating rate was set to 5 °C/min to heat the specimen sufficiently and to avoid as much as possible the thermal shock produced by a rapid heating rate [21,22,23]. A constant temperature was maintained for 2 h after reaching the target temperatures, which were 25, 150, 300, 450, and 600 °C, respectively. The heated specimens were immersed in circulating cold water for 1 h for rapid cooling. After that, Brazilian splitting tests were carried out on the heat-treated granite specimens using a rock mechanics testing system, model TFD 20D, which is produced by Changchun Keyi Test Instrument Co., Ltd. (Changchun, China) to obtain the evolution characteristics of tensile strength with heat treatment temperature. The test was conducted using displacement-controlled loading with a set loading rate of 0.1 mm/min. The heat treatment process and mechanical testing system are shown in Figure 2.

The tensile strength of the Brazilian disc specimens under the load P was calculated using the following equation:

$${\sigma }_{\mathrm{t}}={\displaystyle \frac{P_{\max }}{\pi Rt}}$$

(1)

where ${\sigma }_{\mathrm{t}}$ is the tensile strength, MPa; $P_{\max }$ is the peak load, kN; R and t are the radius and thickness of the Brazilian disc specimen, respectively, mm.

3. Numerical Modeling Methodology

3.1. Grain-Based Model

Granite is a polycrystalline material, with varying grain sizes and types, leading to complex mechanical behavior. In this work, the simulations were conducted using Particle Flow Code in two dimensions (PFC 2D) version 6.0, a commercial software developed by Itasca Consulting Group, Inc. (Minneapolis, MN, USA) designed to simulate the mechanical behavior of materials, particularly granular, rock, and rock-like materials [24]. It models material behavior by representing the medium as an assembly of circular or polygonal particles that interact through contact forces and bonds, making it highly suitable for studying fracture processes, stress distribution, and particle motion. The GBM within the PFC 2D was employed to simulate natural rock behavior, which allows the rock to be modeled as an assembly of discrete grains, each with unique mechanical properties, such as stiffness and strength [25]. This mirrors the physical properties of real granite, enabling a realistic simulation of how cracks develop and propagate. Furthermore, since bonds between grains can break and form cracks, the GBM can reproduce complex fracture networks and failure mechanisms, such as tensile, shear, and mixed-mode fractures [26,27,28]. The specific steps to build a GBM are as follows:

(1) Based on the mineral composition of the granite specimens, an initial particle packing with a radius in the range of 1.0–1.8 mm was generated within a rectangular frame, which determined the size of the various mineral crystals (Figure 3a).

(2) Based on the position and radius of the initial particles, Voro++ was used to generate rigid blocks to form polygonal grain structures (Figure 3b). Voro++ version 0.4.6 is an open-source software library for computing the Voronoi diagram, a widely used tessellation with applications in mesh generation and numerical modeling.

(3) The initial particle packing was removed and then regenerated with a polygonal grain structure with a smaller radius (0.15–0.25 mm) to generate disc numerical specimens (Figure 3c). Before assigning micro-parameters, floating particles were removed from the numerical specimens and a confining pressure was applied to ensure that the particles settled into stable contact with one another.

(4) The intergranular contacts were grouped according to the location and grouping of the grains. The parallel bond contact model (PBC) and smooth joint contact model (SJC) were assigned to intra- and intergranular contacts, respectively. The modeled numerical specimens contained a total of 12,280 particles, 22,609 PBCs, and 7213 SJCs (Figure 3d).

After micro-parameters were assigned to the model, the confining pressure applied to the specimen was unloaded, and thermal simulation was achieved by assigning different thermal expansion coefficients to the mineral particles. A steady-state heating method was used to ensure uniform heating of the specimen. The steps of thermal simulation were as follows:

(1) The wall unit around the numerical specimen was removed.

(2) The numerical specimens were heated in cycles with a gradient of 20 °C until the target temperature was reached (i.e., 150, 300, 450, and 600 °C). The time step for each cycle was based on the model reaching static equilibrium.

(3) After reaching the target temperature, the temperature of the numerical specimen was rapidly reduced to room temperature, and the model was brought to static equilibrium again by iterative calculations.

3.2. Model Calibration

Based on the results of the Brazilian splitting test, particle properties, bond strengths, and contact models were adjusted to better match experimental data and fracture behavior. The model calibration was performed mainly for specimens at room temperature, and the numerical results for heat-treated specimens were realized by thermal simulation, which eliminates the need for recalibration at various temperatures.

The most commonly used method for calibrating micro-parameters in PFC is the trial-and-error method, which involves iteratively adjusting model parameters until the simulated rock behavior matches the target mechanical properties observed in real specimens. This process begins with selecting initial parameter values, running simulations, and comparing the output to desired properties. If the simulated results deviate from the target values, parameters are adjusted, and the simulation is repeated. This cycle continues until an acceptable match is achieved, ensuring that the micro-parameters accurately represent the mechanical behavior of rocks. Therefore, the load–displacement curves and tensile strengths of the room-temperature and heat-treated specimens were combined, and the micro-parameters of the numerical specimens were iteratively adjusted using the trial-and-error method until satisfactory micro-parameters were obtained.

Table 1. Micro-parameter calibration results of the numerical model.

Item	Micro-Parameter (Unit)	Value
		Plagioclase	Quartz	K-Feldspar	Biotite
Basic parameters	Mineral content (%)	70.7	17.0	10.0	2.3
	Minimum radius of particle (mm)	0.15	0.15	0.15	0.15
	Particle radius ratio	1.66	1.66	1.66	1.66
	Particle density (kg·m−3)	2600	2650	2600	2850
	Particle friction coefficient	1.2	1.2	1.2	1.2
Mineral grains (PBC model)	Contact normal to shear stiffness ratio	1.7	1.0	1.6	1.1
	Particle–particle contact modulus (GPa)	1.0	2.0	0.7	0.5
	PBC normal to shear stiffness ratio	1.7	1.0	1.6	1.1
	PBC modulus (GPa)	1.0	2.0	0.7	0.5
	PBC radius multiplier λp	0.60	0.60	0.60	0.60
	PBC tensile strength (MPa)	33.6	36.4	59.5	28.0
	PBC cohesion (MPa)	48	52	45	40
	PBC friction angle (°)	30	30	30	30
Grain boundaries (SJC model)	SJC radius multiplier λs	0.60
	SJC normal stiffness factor	0.60
	SJC shear stiffness factor	0.80
	SJC tensile strength (MPa)	10
	SJC cohesion (MPa)	30
	SJC friction angle (°)	30
	SJC friction coefficient	1.2
Thermal parameters	Specific heat (J·Kg−1·K−1)	1015
	Thermal conductivity (W·m−1·K−1)	3.5
	Thermal expansion coefficients (10−6·K−1)	14.1	24.3	8.7	3.0

presents the calibration results for the micro-parameters.

3.3. Model Validation

The model was validated using the load–displacement curves and tensile strength obtained from the Brazilian splitting tests, and a comparison between the numerical and experimental results is shown in Figure 4. It is worth noting that in PFC 2D the default thickness of the numerical specimen is 1 m, while the thickness of the real specimen is 25 mm, so the numerical result of the load is about 40 times the experimental result [29]. The match between numerical and experimental results is relatively good for all specimens except the 150 °C heat-treated specimen. The reasons for this mismatch are multiple, including simplifying assumptions about the constitutive relations and micromechanical parameters on which the numerical simulations are based, as well as the complexity of mineral grain growth, phase transformations, and thermal cracking under high-temperature conditions [19,20]. Overall, the developed numerical model is effective in simulating Brazilian splitting tests on heat-treated granite specimens.

4. Results and Discussions

4.1. Thermal Cracking Behavior

Table 2. Evolution of thermal crack distribution.

T (°C)	Heating-Induced Cracks	Cooling-Induced Cracks
150
300
450
600

Red, green, and black represent intergranular tensile cracks, intragranular tensile cracks, and intergranular shear cracks, respectively.

shows the evolution of thermal crack distribution during the heating and cooling phases, where red, green, and black colors represent intergranular tensile cracking, intragranular tensile cracking, and intergranular shear cracking induced by heat treatment, respectively. The variation in the number of different types of thermal cracking with heat treatment temperature is shown in Figure 5. No thermal cracking was produced when the heat treatment temperature was 150 °C. Existing studies have shown that thermal expansion at this temperature results in a closer bonding between mineral particles, which is one of the main reasons for the enhanced strength of the specimens [30,31,32]. When the heat treatment temperatures were 300 °C and 450 °C, the mutually independent cracks formed during the heating phase propagated further and coalesced with the newly sprouted cracks in the cooling phase, resulting in about five times as many cracks in the cooling phase as in the heating phase. When the heat treatment temperature was 600 °C, the phase transformation of quartz crystals led to a higher number of thermal cracks in the heating phase than in the cooling phase. In addition, thermal cracking mainly consisted of inter- and intragranular tensile cracks, with intergranular tensile cracks being dominant. As the temperature increased, the number of inter- and intragranular tensile cracks increased, with a small number of intergranular shear cracks appearing at 600 °C.

4.2. Mechanically Induced Cracking Behavior

The crack propagation paths of the specimens after the Brazilian splitting test are shown in

Table 3. Comparison of crack propagation paths between experimental and numerical results.

T (°C)	Test Image	Numerical Image	Binary Image
25
150
300
450
600

In the numerical images, red represents thermal cracks, while yellow, blue, and green represent intragranular tensile cracks, intergranular tensile cracks, and intergranular shear cracks induced by loads. In the binary images, black represents thermal cracks and mechanically induced cracks.

. In the numerical images, red represents thermal cracks, while yellow, blue, and green represent intragranular tensile cracks, intergranular tensile cracks, and intergranular shear cracks induced by loads. The crack propagation path of the numerical specimens matched the experimental results, that is, splitting occurred along the loading centerline. However, there is some deviation between the numerical and experimental results of crack shapes. The reason for this is that granite is a naturally inhomogeneous material with strong anisotropic characteristics and exhibits more complex fracture behavior under heat treatment. Under Brazilian splitting load, the failure of the specimen is dominated by tensile fracture, with the number of intragranular tensile cracks dominating, followed by intergranular tensile cracks, the presence of a small number of intergranular shear cracks, and the absence of intragranular shear cracks. In addition, thermal cracking affects the crack propagation path of mechanically induced cracks, as evidenced by the fact that the crack propagation path becomes more tortuous with increasing heat treatment temperature.

Figure 6 shows the number and percentage of mechanically induced cracks at different heat treatment temperatures. Within the temperature range of 25–300 °C, the number and percentage of mechanically induced cracks of various types remained relatively constant. When the heat treatment temperature was 450–600 °C, the mechanically induced cracks tended to propagate more along the path of widely distributed thermal cracks, and the strength and interparticle bonding properties of the heat-treated specimens decreased significantly. The propagation and coalescence of those cracks created mechanically induced cracks on a larger scale, leading to a reduction in the number of cracks.

Figure 7 presents the rose diagrams for the orientations of different types of mechanically induced cracks, where the radius indicates the crack count, and the direction is defined as the angle measured clockwise from the horizontal in 10-degree intervals. Under Brazilian splitting load, the specimens split in the vertical direction, so that most of the mechanically induced cracks were concentrated in the 70–100° range of orientation.

In particular, the orientation of intergranular shear cracks is predominantly concentrated between 60° and 120°, whereas the orientations of inter- and intragranular tensile cracks are more dispersed.

4.3. Fractal Dimension Characteristics of Cracks

4.3.1. Calculation Method for Fractal Dimension of Cracks

The fractal dimension quantifies the complexity of a fractal object by measuring how the detail in a pattern changes with the scale at which it is measured. The box-counting method is widely used for calculating the fractal dimension, especially for irregular shapes or patterns [33,34]. The basic idea of the box-counting method is to cover the object with boxes of varying sizes and analyze how the number of boxes needed to cover the object changes as the size of the boxes decreases. As shown in Figure 8, the image of the numerical specimen containing cracks is first converted to a binary image using MATLAB version R2021b. Then, a grid of boxes of size δ is overlaid on the binary image and the number of boxes N(δ) containing a part of the object is counted, and the process is repeated with different box sizes. If the object has fractal characteristics, the relationship between N(δ) and δ can be given by the following equation:

$$N(\delta )\sim {\delta }^{-D}$$

(2)

where D is the fractal dimension of the object.

Therefore, the D obtained by the box-counting method can be calculated as follows:

$$D=-\underset{\delta \to 0}{\lim }{\displaystyle \frac{\ln N(\delta )}{\ln \delta }}$$

(3)

The ln N(δ)–ln δ relationship curves for the cracks of heat-treated specimens after loading are shown in Figure 9. The ln N(δ)–ln δ curves of each group of specimens are straight lines, and the correlation coefficient R² is close to 1. The test data show a good linear correlation. This indicates that the crack distribution of the heat-treated specimens has obvious fractal characteristics, so the fractal dimension can be used to quantitatively characterize the cracks [35,36].

4.3.2. Relationship Between Fractal Dimension of Cracks and Tensile Strength

Table 4. Calculation results for the numerical model based on the GBM.

T (°C)	Accumulated Number of Cracks			Fractal Dimension			Tensile Strength (MPa)
	Heating	Cooling	Loading	Heating	Cooling	Loading
25	0	0	806	0	0	1.32	6.73
150	0	0	890	0	0	1.33	7.50
300	80	453	1297	0.42	0.84	1.42	6.70
450	382	2260	2977	0.80	1.26	1.54	5.80
600	3288	5291	5697	1.33	1.49	1.65	2.46

lists the cumulative number of cracks and the fractal dimension for each group of specimens. When the fractal dimension is less than 1, the cracks have a point-like sporadic distribution. When the fractal dimension is greater than 1, the cracks are spatially distributed between one and two dimensions, with complex point-like and line-like discontinuous distributions. In the heating and cooling phases, the fractal dimension gradually increases from 0 to more than 1 with the increase of the heat treatment temperature, accompanied by the propagation and coalescence of thermal cracking. In the loading phase, the fractal dimension gradually increases with the increase of the heat treatment temperature, and the morphology of the mechanically induced cracks becomes more complex, indicating that the specimen is more damaged or broken.

The relationship between fractal dimension and tensile strength is shown in Figure 10. It can be seen that the fractal dimension of the mechanically induced cracks of the heat-treated specimens at 150 °C is close to that of the room-temperature specimens, and the tensile strength increases slightly. As the heat treatment temperature increases, the tensile strength gradually decreases, and the fractal dimension of the cracks gradually increases. It is interpreted that the thermal cracks increase the degree of damage of the specimens, leading to the fact that the mechanically induced cracks are more inclined to propagate along the path of the thermal cracks. Higher heat treatment temperatures result in the formation of more intensive thermal cracking, leading to more fragmentation of the specimen after failure and larger fractal dimension of the cracks. Therefore, the tensile strength can be characterized by the fractal dimension characteristics of the cracks.

5. Conclusions

Brazilian splitting tests were conducted on heat-treated granite specimens at different temperatures, followed by modeling and validation of a GBM containing four different mineral particles. The number and distribution characteristics of microcracks during the heating, cooling, and loading phases were analyzed, and the relationship between the fractal dimension of the cracks and the tensile strength was discussed. The following conclusions were drawn.

(1) The microcracking behavior and tensile fracture properties of heat-treated granite specimens can be effectively simulated by the GBM, which accounts for different mineral compositions and thermal expansion characteristics.

(2) Thermal cracking primarily consists of inter- and intragranular tensile cracks, with intergranular tensile cracks being dominant. As the heat treatment temperature increases, the number of both inter- and intragranular tensile cracks increases.

(3) Under Brazilian splitting load, the failure of the specimens was dominated by intragranular tensile cracking, followed by intergranular tensile cracking, with a small amount of intergranular shear cracking.

(4) The distribution of both thermal cracking and mechanically induced cracking in the heat-treated granite specimens exhibited clear fractal characteristics. As the heat treatment temperature increased, the fractal dimension of the cracks gradually increased, while the tensile strength gradually decreased.

(5) Higher heat treatment temperatures resulted in denser thermal cracks, causing the crack propagation paths of the mechanically induced cracks to become more tortuous, and the fractal dimension of the total cracks to increase.