3D Numerical Prediction of Thermal Weakening of Granite under Tension

Abstract

This paper deals with numerical prediction of temperature (weakening) effects on the tensile strength of granitic rock. A 3D numerical approach based on the embedded discontinuity finite elements is developed for this purpose. The governing thermo-mechanical initial/boundary value problem is solved with an explicit (in time) staggered method while using extreme mass scaling to increase the critical time step. Rock fracture is represented by the embedded discontinuity concept implemented here with the linear (4-node) tetrahedral elements. The rock is modelled as a linear elastic (up to fracture by the Rankine criterion) heterogeneous material consisting of Quartz, Feldspar and Biotite minerals. Due to its strong and anomalous temperature dependence upon approaching the α-β transition at the Curie point (~573 °C), only Quartz in the numerical rock depends on temperature in the present approach. In the numerical testing, the sample is first volumetrically heated to a target temperature. Then, the uniaxial tension test is performed on the cooled down sample. The simulations demonstrate the validity of the proposed approach as the experimental deterioration, by thermally induced cracking, of the rock tensile strength is predicted with a good accuracy.

1. Introduction

High temperature conditions often prevail in geotechnical and mining applications [1,2,3]. For this reason, rocks, especially those containing Quartz mineral, are prone to be affected by elevated temperatures. For Quartz bearing rocks, this is due to its α-β transition at about 573 °C [4]. The temperature effects on rock have naturally drawn considerable attention in the literature, both experimentally [1,5,6,7,8,9,10,11,12,13,14] and numerically [11,12,13,14,15,16,17].

Under elevated temperatures, the mechanical properties (e.g., the Young’s modulus and strength) of rock decrease, while some of the thermal properties increase (the thermal expansion coefficient and the specific heat) and others decrease (the thermal conductance) [7,11]. The degradation of rock strength has an extreme importance from the rock mechanical applications point of view. The major mechanism of the thermal weakening of rock strength is the thermally induced cracking, which in turn originates from the rock mineral heterogeneity [13]. It is the mismatch of the elastic properties at the grain boundaries of different minerals that leads to thermal stresses and generates cracking inside and between the grains even under uniform temperature fields (i.e., no temperature gradients present). For rocks with Quartz as a one the constituent minerals, e.g., granite, which is studied in the present paper, this heterogeneity becomes increasingly pronounced due to the exceptional behavior of Quartz. Namely, the thermal expansion of Quartz is highly nonlinear up to the α-β transition temperature. In contrast, the rest of the granite forming minerals (Fledspars and Micas) exhibit linear temperature dependence in their thermal expansion [7].

Predictive modelling of rock fracture due to elevated temperatures is an important task in geosciences. It provides knowledge on the failure mechanisms impossible or too expensive to reach by experimentation or in-situ testing. However, predicting temperature effects in rock by computational models is challenging due to cracking phenomenon, which requires numerical description of displacement discontinuities. In computational mechanics, the two main methods to simulate cracking are the finite element method (FEM) and the discrete element method (DEM). The reader can refer to [18,19] for reviews on numerical rock mechanics. As a summary of these methods, it can be stated that the DEM is naturally suited for fracture and fragmentation modelling (see [20] for an example). However, a critical drawback is the computational labor required by tracking and updating the contact configurations and neighbors of the particles that, in contrast to FEM, need not to retain their neighbors during the analysis. Within the continuum approach, such as the classical FEM, fractures (cracks) can be modelled only in an averaged manner, i.e., as localized deformation, by damage and/or plasticity models. However, FEM based models are computationally efficient and relatively simple in their calibration of the material parameters. Considerable efforts have, therefore, been devoted to enhancing the fracture description capabilities of the classical FEM. The two most notable results of this work are the embedded discontinuity FEM [21] and the extended FEM [22]. The former is chosen in the present study due to its considerably simpler implementation [23].

The weakening effect due to thermal cracking on the tensile strength of granite is numerically investigated in this paper. To this end, a staggered solution method is adopted for solving the governing thermo-mechanical problem. The granite is modelled as a linear elastic and heterogenous material up to a fracture, which is described by the embedded discontinuity finite elements. A full 3D extension of a previous model by Saksala [24] is presented. The granite elasticity properties, in contrast to [24], and tensile strength are taken as heterogeneous, i.e., individual for each rock forming mineral. Moreover, only the thermal expansion coefficient of Quartz depends on temperature in the present model, which is due to the much stronger thermal expansion behavior of Quartz (as noted above) than the other granite-forming minerals. This simplified approach is numerically validated by genuine 3D simulations.

2. Numerical Methods

2.1. Failure Model for Rock

Let a 3D domain be discretized with linear tetrahedral elements with some of them intersected by displacement discontinuities (cracks). Figure 1a shows a tetrahedron sliced by two intersecting discontinuity planes. Two such crack planes are required when modelling the thermal cracking induced weakening of the tensile strength. The reason for this is that the thermally induced cracks may have unfavorable orientations (perpendicular) with respect to the loading direction in the consequent uniaxial tension test. Thereby, as cracks cannot rotate, an element with a crack normal perpendicular to the loading axis would generate spurious stresses. This is prevented by the introduction of a new crack.

As granite is a brittle material at room temperature, the small deformation assumption is justified. Under these circumstances, the displacement and strain fields for a 4-node tetrahedron with two displacement discontinuities read

$$u\left(x\right)={\sum }_{k=1}^4{N}_k\left(x\right)u_k^e+{{\sum }^{ }}_{i=1}^2{M}_{{\Gamma }_{di}}\left(x\right){\alpha }_{di} \mathrm{with} M_{{\Gamma }_{di}}\left(x\right)=H_{{\Gamma }_{di}}\left(x\right)-{\phi }_{{\Gamma }_{di}}\left(x\right),\left(i=1,2\right)$$

(1)

$$\epsilon \left(x\right)={\sum }_{k=1}^4{\left(\nabla N_k\otimes u_k^e\right)}^{sym}-{{\sum }^{ }}_{i=1}^2\left({\left(\nabla {\phi }_{{\Gamma }_{di}}\left(x\right)\otimes {\alpha }_{di}\right)}^{sym}+{\delta }_{{\Gamma }_{di}}{\left(n_i\otimes {\alpha }_{di}\right)}^{sym}\right),$$

(2)

where ${\alpha }_{di}$ is the crack opening (displacement jump) vector for crack i, and $N_k$ and $u_k^e$, being functions of placement vector x, are, respectively, the regular shape functions for the linear (4-node) tetrahedral element and nodal displacements. The crack normal is denoted by $n_i$ (for crack i). In addition, $H_{{\Gamma }_{di}}$ and ${\delta }_{{\Gamma }_{di}}$ denote the Heaviside function and its gradient (∇), the Dirac delta function, respectively. Consistent with the constant strain of the linear tetrahedron element, the crack opening vector is also assumed elementwise constant. Thereby, $\nabla {\alpha }_{di}\equiv 0$, and (2) readily follows from taking the gradient of (1). Moreover, the terms ${\delta }_{{\Gamma }_{di}}{\left(n_i\otimes {\alpha }_{di}\right)}^{sym}$ (where (⊗)^sym denotes the symmetric part of the tensor product), containing the Dirac delta function, are zero only outside the discontinuity plane and can thus be ignored when solving the discretized equations of motion at the global level.

Function $M_{{\Gamma }_{di}}$ in (1) restricts the effect of ${\alpha }_{di}$ inside the corresponding finite element. This simplifies substantially the implementation of the kinematics, as it removes the need to treat the essential boundary conditions related to crack opening at the model boundaries (see Figure 1b for the 1D case). Function ${\phi }_{{\Gamma }_{di}}$ in $M_{{\Gamma }_{di}}$ is selected from the shape function combinations of the element so that

$$\nabla {\phi }_{{\Gamma }_{di}}=\arg \left(\underset{k=1,2,3}{\max }\frac{\left|{{\sum }^{ }}_{n=1}^k\nabla N_n·n_i\right|}{\Vert {{\sum }^{ }}_{i=1}^k\nabla N_n\Vert }\right),\left(i=1,2\right)$$

(3)

This criterion thus selects ${\phi }_{{\Gamma }_{di}}$ so that $\nabla {\phi }_{{\Gamma }_{di}}$ is as parallel to $n_i$ as possible. Figure 1b illustrates the displacement decomposition (1) and the relevant functions in 1D case with a single linear truss element pulled at the right end. For this case, the regular nodal displacement is $u_{reg}=N_1{u}_1+N_2{u}_2$, and function φ is the interpolation (shape) function N₂. The effect of function $M_{{\Gamma }_d}$, i.e., the restriction of the displacement jump, ${\alpha }_d$, inside the element, is readily observed.

The FEM formulation of the above kinematics is performed within the framework of the enhanced assumed strain concept (EAS), see the details in [23,25,26]. The final equations for the spatially discretized thermo-mechanical problem are

$${\int }_{{\Omega }_e}\rho N_i{N}_j{\ddot{u}}_{j}d\Omega +{\int }_{{\Omega }_e}\sigma \cdot \nabla N_{i}d\Omega =0,i,j=1\dots N_{nodes}$$

(4)

$${\int }_{{\Omega }_e}c\rho N_i{N}_j{\dot{\mathsf{\theta }}}_{j}d\Omega +{\int }_{{\Omega }_e}k\nabla N_i\nabla N_j{\mathsf{\theta }}_{j}d\Omega -{\int }_{{\Omega }_e}{Q}_{\mathrm{int}}{N}_{j}d\Omega =0,i,j=1\dots N_{nodes}$$

(5)

$${\varphi }_{di}\left(t_{{\Gamma }_{di}}\right)=0,t_{{\Gamma }_{di}}=\sigma ·n_i,i=1,2 (\mathrm{if} \mathrm{element} \mathrm{has} \mathrm{crack}(\mathrm{s}\left)\right)$$

(6)

$$\sigma =E:\left(\widehat{\epsilon }-{{\sum }^{ }}_{i=1}^2{\left(\nabla {\phi }_{{\Gamma }_{di}}\otimes {\alpha }_{di}\right)}^{sym}-{\epsilon }_{\mathsf{\theta }}\right),$$

(7)

where $\widehat{\epsilon }={\sum }_{k=1}^4{\left(\nabla N_k\otimes u_k^e\right)}^{sym}$, and ${\epsilon }_{\theta }=\alpha \Delta \theta I$ is the thermal strain with α being the thermal expansion coefficient, $\Delta \theta$ is the temperature change, and I is the unit tensor. Moreover, ${\ddot{u}}_j$ is the acceleration vector, $N_{nodes}$ is the number of nodes, and $N_i$ is the shape function of node i. Furthermore, $\sigma$ is the stress tensor, ${\dot{\mathsf{\theta }}}_j$ is the time derivative of the temperature vector, ${\mathsf{\theta }}_j$, at node j, $\rho$ and c are the density and the specific heat capacity, respectively. Finally, k is the conductivity, and $Q_{\mathrm{int}}$ is the internal heat production.

Equation (4) is the discretized balance of linear momentum, and Equation (5) is the discretized heat balance. Equation (6)₁ is the failure criterion, where ${\varphi }_{di}$ is the loading function and $t_{{\Gamma }_{di}}$ is the traction vector. Furthermore, Equation (6)₂ is the strong form of traction balance over the crack i surface. Finally, Equation (7) is the constitutive relation for rock with $E$ being the stiffness tensor. The only boundary condition in Equations (4)–(7) is the external volumetric heat influx $Q_{\mathrm{int}}$, which accounts for the heating in an oven. All other types of heat generation are ignored as insignificant compared to the external heat influx [15].

In tension, the rock material behaves in linear elastic manner until the tensile strength is reached. Upon violation of the Rankine criterion (the first principal stress criterion), a crack (displacement discontinuity) is introduced with a normal complying to the first principal direction. A glance at Equations (4), (6) and (7) reveals their formal similarity to the corresponding equations in plasticity theory. Consequently, the stress integration of the evolution equations can be performed with the usual stress return mapping algorithms [15,23,25]. The following traction-separation model is employed to control the crack softening (opening).

$${\varphi }_{di}\left(t_{{\Gamma }_{di}},{\kappa }_{di},{\dot{\kappa }}_{di}\right)=n_i·t_{{\Gamma }_{di}}+{\left({\left(m_{1i}·t_{{\Gamma }_{di}}\right)}^2+{\left(m_{2i}·t_{{\Gamma }_{di}}\right)}^2\right)}^{1/2}-\left({\sigma }_{ti}+q_{di}\left({\kappa }_{di},{\dot{\kappa }}_{di}\right)\right),$$

(8)

$$q_{di}=h_{di}{\kappa }_{di}+s_d{\dot{\kappa }}_{di},h_{di}=-g_d{\sigma }_{t}exp\left(-g_d{\kappa }_{di}\right),g_d=\frac{{\sigma }_t}{G_{Ic}},$$

(9)

$${\dot{t}}_{{\Gamma }_{di}}=-E:\left(\nabla {\phi }_{{\Gamma }_{di}}\otimes {\dot{\alpha }}_{di}\right)·n_i,$$

(10)

$${\dot{\alpha }}_{di}={\dot{\lambda }}_{di}\frac{\partial {\varphi }_{di}}{\partial t_{{\Gamma }_{di}}},{\dot{\kappa }}_{di}=-{\dot{\lambda }}_{di}\frac{\partial {\varphi }_{di}}{\partial q_{di}},$$

(11)

$${\dot{\lambda }}_{di}\ge 0,{\varphi }_{di}\le 0,{\dot{\lambda }}_{di}{\varphi }_{di}=0, (\mathrm{i}=1,2)$$

(12)

The symbols herein are as follows: ${\kappa }_{di},{\dot{\kappa }}_{di}$ are the internal variable and the rate of it related to the softening law $q_{di}$ for the discontinuity i; ${\sigma }_{ti}$ is the tensile strength; s_d is the viscosity modulus; $h_{di}$ is the softening modulus of the exponential softening law; parameter $g_d$ controls the initial slope of the softening curve. Its calibration is via the mode I fracture energy G_Ic. Moreover, ${\dot{\lambda }}_{di}$ is the crack opening increment. The loading function ${\varphi }_{di}$ in (8) is a sum of the normal and shear terms with β being the shear retention factor. Finally, Equation (12), i.e., the Kuhn-Tucker conditions, impose the consistency akin to plasticity theory.

It should be emphasized that a crack is introduced in an element where the stress state violates the Rankine criterion, while the opening (softening) of that crack is controlled by the model in Equations (8)–(12). Moreover, even if the crack initiation is by the Rankine criterion, a crack once introduced, may open in shear mode as well. Now, another criterion for the second crack in an element already having an unfavorably oriented crack is required to carry out the tension test after the thermal treatment. The same criterion, as used in the 2D study [24], originally in [27], is used here as well. A new crack is thus generated in an element having a thermal crack upon fulfilling the condition:

$$\mathrm{If} {\sigma }_1>{\sigma }_t^{*} \& \left|n_1·n^{*}\right|<C_{\alpha } \mathrm{with}$$

(13)

$${\sigma }_t^{*}=\left(1-\left|n_1·n^{*}\right|\right){\sigma }_t^{ic}+\left|n_1·n^{*}\right|{\sigma }_{t0}$$

(14)

where $n_1$ is the normal of the first crack, ${\sigma }_1$ and $n^{*}$ are the present 1st principal stress and the corresponding principal direction. The tensile strength, ${\sigma }_t^{*}$, for the second crack is a convex combination of the strength of the initial crack, ${\sigma }_t^{ic}$, and the intact tensile strength ${\sigma }_{t0}$. The second, new, crack is introduced (only once) when the angle between the old, first, crack normal and the new principal direction is greater than $\alpha =acos\left(C_{\alpha }\right)$. $C_{\alpha }$ here is a parameter for which a value $1/\sqrt{2}$ (corresponding to $\alpha =45°$) is used in this study.

2.2. Staggered Explicit Scheme for Solving Global Thermo-Mechanical Problem

The discretized problem (Equations (4) and (5)) is solved explicitly in time with a staggered scheme [24]. Staggered solution means that, first, the new temperature vector is solved while keeping the displacement vector fixed, and second, the displacement vector is solved with fixed temperature vector. Equations (8)–(10) are solved with the standard stress return mapping. Figure 2 illustrates the flow of the solution process with the explicit Euler time integrator.

In Figure 2, the relations for solving the temperature and the displacement are the matrix forms of Equations (4) and (5). In addition, the asterisk * means that the trial traction is in the elastic region, or the trial stress is below the tensile strength. Mass scaling is used to increase the conditional time step of the explicit time stepping. This is needed since the time scale of the mechanical problem is many orders of magnitude smaller than that of the thermal problem. Finally, it is noted that the explicit time marching scheme in Figure 2 is employed in the simulations of uniaxial tension test on the intact and heated numerical rock.

3. Numerical Examples

3.1. Numerical Rock Description

Table 1. Material properties for simulations.

Property/Mineral	Quartz	Feldspars	Biotite
Young’s modulus E [GPa]	90	70	40
Poisson’s ratio ν	0.1	0.3	0.2
Tensile strength σt0 [MPa]	14	11	7
Mode I fracture energy GIc [J/m2]	40	40	28
Density ρ [kg/m3]	2650	2650	2650
Thermal expansion coeff. α0 [1/K]	$8\times {10}^{-6}$	$8\times {10}^{-6}$	$8\times {10}^{-6}$
Specific heat capacity c [J/kgK]	820	820	820
Thermal conductivity k [W/mK]	2.6	2.6	2.6
Fraction [%]	33	60	7

lists the mechanical and thermal property values for numerical granite taken from [11,28,29].

The rock heterogeneity is described by spatially random clusters of finite elements assigned to represent the constituent minerals in Table 1. The randomness in this heterogeneity description is generated as follows: An array with a length of the number of elements in the mesh containing integers 1, 2 and 3 (each coding a mineral in the mesh) is first created. The fraction (percentage) of each integer in this array corresponds to the fraction of the mineral in the rock. Next, this array is randomly permutated. A spatially random mineral mesotexture is obtained, as each entry in the array is also implicitly an element number in the mesh (1st entry is the element number 1 in the mesh etc.). Three realizations (i.e., developed with three different random permutations of the integer array just described) of numerical rock mesotextures (consisting of 206,617 elements) are shown in Figure 3. For example, the fraction of elements representing Quartz mineral in each numerical rock sample, having thus the material properties specific to Quartz in Table 1, is 0.33⋅206,617 ≈ 68,183. Finally, the shear retention factor in Equation (8) is β = 1, and the viscosity (Equation (9)) is set to s_d = 0.001 MPa⋅s/m.

With respect to temperature dependence of the rock properties, a simplified approach is adopted. Accordingly, only the thermal expansion is taken as temperature dependent. The experimental setting is as follows: First, the specimen is heated slowly to the desired temperature. Then, it is left to cool down to room temperature. The uniaxial tension test is finally carried out at the room temperature. This thermal pre-treatment generates cracks into the sample causing reduction of its tensile strength.

Heterogeneous brittle rocks have inherent microflaws. Therefore, their tensile strength is an emerging property measured for the laboratory sample, not for a material point. However, in the numerical modelling, the constitutive law is described at the material point level. Thereby, it would be circular reasoning to input the experimental tensile strength at certain temperature into the material point level failure model, and then predict that same strength—the correct outcome would ensue as a matter of course.

For this reason, only the thermal expansion coefficient of Quartz is taken to depend on temperature by

$${\alpha }_q\left(\theta \right)={\alpha }_{q0}\left(1+\frac{1.5}{823-293}\left(\theta -293\right)\right), [1/\mathrm{K}]$$

(15)

where ${\alpha }_{q0}$ (= $8\times {10}^{-6}$ K⁻¹) is the thermal expansion coefficient at room temperature. Equation (15) assumes a linear thermal expansion dependence for Quartz in between 20 and 550 °C (where (${\alpha }_{q550}$ = $2\times {10}^{-5}$ K⁻¹)), which is not realistic, as it was already noted that Quartz thermal expansion is highly nonlinear on approaching the Curie point [7]. This assumption is a feature of the present approach, which, however, will prove to predict the thermal weakening effect with a good accuracy.

3.2. Numerical Rock Heat Treatment

Uniform heating up to 300 °C and 500 °C is carried out here on the numerical rock samples. The reason for choosing these temperatures is that they often appear in the experiments [11]. Moreover, significant reduction in the tensile strength of granite only appears beyond 200 °C. Finally, after the Curie point of granite, i.e., 573 °C, the α-Quartz changes to β-Quartz involving a phase change the proper modelling of which is not trivial.

Mass scaling is used here to increase the critical time step of the explicit time stepping in Figure 2. This is enabled by the non-inertial nature of the uniform heating of rock samples. More precisely, the critical time step of the explicit time marching can be increased to 100-fold by applying a 10,000-fold density. Moreover, to secure a homogenous temperature field in the rock sample, volumetric heating is applied here with setting $Q_{\mathrm{int}}$ = $1\times {10}^{10}$ W/m³ at each node of mesh. With this setting, the desired temperature 500 °C is reached at every node of the mesh in ~20,000 time-steps, which translates to a heating time of 0.1 s. The critical time step for the numerical rock sample 1 (NumRock1) is ~$5\times {10}^{-8}$ s (with the real density). Simulation results for this sample with heating up to 500 °C is shown in Figure 4.

The uniform volumetric heating applied induces considerable cracking in the rock sample, as attested in Figure 4e, where only every 100th crack normal of the total amount of the 160,634 cracks is plotted. The crack normal orientation is quite randomly distributed, i.e., no trends or preferred orientations can be observed. According to Figure 4d, showing the temperature and cumulative crack number evolution, massive cracking starts when the temperature reaches 100 °C. After reaching about 120,000 cracks, the cumulative crack number evolution proceeds with a much milder rate.

Figure 4b,c show the resulting crack opening magnitude and the stress-like softening variable q_d in Equation (9). Values reaching −10 MPa indicate substantial loss of the traction bearing capacity of the cracks and, consequently, local weakening effect.

There is no reason to show results for the target temperature 300 °C because they are, for this same sample, somewhat similar since the simulation was stopped earlier, at ~0.06 s upon temperature reaching 300 °C. Moreover, the simulation results for the other numerical rock samples are not shown here either due to their similarity to those shown in Figure 4.

3.3. Numerical Tension Tests on Intact and Heat Treated Rock

Next, uniaxial tension tests are performed on the numerical rock samples. A constant velocity boundary condition is applied at the top surface of the numerical rock. Setting the velocity to 5 mm/s turns into a strain rate of 0.1 s⁻¹, which is low enough not to cause significant hardening effects [30]. The results are shown in Figure 5.

All the numerical samples have failed in the experimental transverse splitting mode, as attested in Figure 5a showing the magnitude of crack opening. It should be noted that there is a huge number of cracks all over the sample volume (Figure 5c), with a normal parallel (or closely so) to the loading axis, but only some of them open to create the final failure plane. Figure 5b shows the cumulative number of cracks, exceeding 100,000, as a function of the axial average strain for NumRock1. It also shows the stress-strain curves, which are strikingly similar; the tensile strength of each sample is 11.1 MPa. The reason for the insignificant deviation is likely the small element size.

Next, the heat-treated numerical samples are submitted to the uniaxial tension test. To replicate the experiments, these simulations are performed on the naturally cooled down samples. The cooling process itself, as it should not generate more cracks (being very slow), is not simulated but imposed. This means that no thermal stresses thus exist in the samples. In addition, the thermally induced cracks are assumed to be closed, i.e., their displacement jump vectors are set to zero (${\alpha }_d=0$). However, their orientations (crack normal) and residual strengths must be carried over to the initial state of the tension simulations. Consequently, the residual strength for an element with a thermal crack is set to ${\sigma }_{t0}+q_d$ where ${\sigma }_{t0}$ is the intact tensile strength of the mineral and $q_d$ is the stress-like softening variable, with the value it reached during numerical heating. The results are shown in Figure 6 and Figure 7 for the target temperatures of 300 °C and 500 °C, respectively.

Heated (to 300 °C) samples all fail with a clear single fracture plane, albeit at different locations from those without heating (compare Figure 6a to Figure 5a). Again, there is not much deviation between the stress-strain curves of the rock samples (Figure 6b) until the post-peak softening part. However, the nonlinear pre-peak part of the curves is, due to the presence of thermal cracks, considerably more pronounced than that in the non-heated case. The resulting tensile strength of the samples is ~9.6 MPa, which is 13.5% smaller than the intact (room temperature) strength. The crack normal plot in Figure 6c shows that a new crack (green color) is initiated in many elements with an unfavorably oriented thermal crack.

Heating up to 500 °C results in more severe weakening of the samples and more nonlinear pre-peak part of the stress-strain response, as attested in Figure 7b. The tensile strength of the samples is ~7.6 MPa, i.e., 31.5% percent drop in the strength. The failure modes still show a clear single transverse crack plane but there is more diffused crack opening (all over the samples) than in the previous cases.

4. Discussion

Some features of the present model and the results are elaborated here. The simulation results are first compared to the experiments collected in [11]. For this purpose, the predicted normalized tensile strengths are plotted at different temperatures in Figure 8.

The present approach predicts the weakening effect correctly. Indeed, the predictions are within the experimental deviation, and very close to the fitted average curve, for several granites. It should be noted that the experimental curve, $f_{{\sigma }_t/{\sigma }_{t0}}=0.9912{\left(1-4.10\mathsf{\theta }/2483.30\right)}^{1/4.10}$, in Figure 8, representing the average for several granites, is averaged over a wide range of test temperatures. Furthermore, only relevant deviation bars (300 °C and 500 °C) are shown in Figure 8.

Then some features of the numerical model are discussed. The present continuum model has some advantages over the discontinuum approach, e.g., the computational effectiveness, the easiness of calibration, and the physical meaning of the model parameters. Furthermore, the inferior fracture description of the finite element method is mended to a good extent by the embedded discontinuity technology.

A staggered explicit approach was employed in the simulations of the slow heating rock samples. This process is quasi-static in nature, which enables using drastic mass scaling for the mass matrix to increase the critical time step to 100-fold in the present case. When combined with intensive volume heating of the sample, realistic 3D simulations can be performed in practical CPU times with fully explicit time marching using a laptop computer.

Heterogeneity of rock leads to degradation of strength at high temperatures. In case of granite, heterogeneity becomes stronger due to the deviant behavior of Quartz at high temperatures. In this study, the tensile strength and elasticity heterogeneity of rock forming minerals was accounted for, while homogenized values were used for the thermal properties. However, probably due to small element size, the (spatially) random element cluster description of granite heterogeneity did not generate appreciable deviation in repeated numerical tests. As for the temperature dependence, it was applied only to the Quartz thermal expansion by a linear fit representing roughly the net effect between Quartz and the rest of the granite forming minerals (Feldspars and Micas). This simplified approach is advantageous in that it is (partly) based on mechanical and thermal properties easily determinable for a rock sample. The measurement of the corresponding properties for rock mineral phases calls for much more elaborate methods. Overall, this approach is up to its task, as the thermal cracking induced rock strength degradation was successfully predicted. It is worth of emphasizing that the present model predicts the weakening effect in a non-circular way, meaning that experimental data on the temperature dependence of the tensile strength of granite is not used as a model input.

Comparing this study to that by Saksala [24], it can readily be stated that the present 3D version better represents reality, unlike the 2D version [24]. For this reason, the present predictions are more accurate than those of the 2D approach, which, being based on the plane strain assumption, cannot correctly model the cylindrical geometry of the sample used in the present investigation.

However, the present approach does not consider some of the critical textural aspects of rock as a polycrystalline material. Most critically, the grain boundaries were not modelled properly, i.e., as geometric entities that can open as a crack. A possibility to include them within the present continuum approach is to use the cohesive elements between the groups of finite elements representing the grains. However, this strategy would require much more detailed data on the rock mineral texture and the mechanical properties thereof.

5. Conclusions

A 3D extension of a 2D numerical method to predict thermal cracking induced weakening effects in the tensile strength of granitic rock was developed and validated in this paper. The study yields following conclusions:

• An explicit staggered scheme is an effective method to solve the nonlinear coupled problem of thermal cracking due to uniform temperature fields. The present method, based on the embedded discontinuity finite elements to model rock fracture, is computationally fast and has physically meaningful model parameters.
• The non-inertial nature of the slow heating of a rock sample to a homogeneous temperature field enables using drastic mass scaling for the mechanical part of the coupled problem. The mass scaling increases the critical time step of the explicit time stepping for the mechanical part of the governing initial/boundary value problem.
• In terms of modelling, the major factor influencing the thermo-mechanical behavior of granitic rocks seems to be the deviant behavior of Quartz mineral. More specifically, the thermal cracking induced reduction of the tensile strength of granite can be correctly predicted by accounting for the heterogeneity of the mechanical properties (stiffness and tensile strength) and the temperature dependence of Quartz thermal expansion.
• This method, as it does not use any experimental data on the temperature dependence of the tensile strength as a model input, replicates the experimental weakening effect in a non-circular way. The practical benefit of the method is thus a fast and reliable 3D prediction of the tensile strength of granite rock at elevated temperatures with a small set of measurable parameters of a Quartz bearing rock.

Finally, some issues to be addressed in further developments of the present approach are suggested. First, this study considered only the thermal degradation of tensile strength. However, the rock mass in nature is usually under compression. Thereby, the degradation of compressive strength due to thermal cracking needs to be investigated. This is, however, not a trivial task as it requires a compressive failure criterion.