Numerical modeling of thermo‐mechanical failure processes in granitic rock with polygonal finite elements

This paper considers numerical modeling of intensive heating induced thermo‐mechanical failure processes in granitic rock. For this end, a numerical method based on polygonal finite elements and a damage‐plasticity model is developed. A staggered scheme is employed to solve the global thermo‐mechanical problem. The rock failure is described by a Rankine‐Mohr‐Coulomb plasticity model with separate scalar damage variables for tension and compression. Consistent tangent operator is derived for this model. Special attention is given to the temperature dependence of the thermo‐mechanical material properties of heterogeneous rock. In the numerical examples, the method is first verified with an analytical solution of thermal stresses in a hollow cylinder, and then qualitatively validated with the problems of thermal cracking of concentric cylinders and uniaxial compressive test on rock under elevated temperatures. Finally, the method is applied in novel simulations inspired by the degradation of sauna stones under slow heating‐rapid cooling and the comminution by rapid heating‐cooling cycles.

been devoted to thermo-mechanical failure processes of rocks. [11][12][13][14][15][16][17][18][19][20][21] While these studies have their merits, they do not present a method versatile enough, capable of solving the coupled thermo-mechanical problem under both mechanical and thermal loadings of short and long duration, while accounting for the rock mesostructure and the heterogeneity thereof. This is the topic of present study.
A continuum approach based on the finite element method (FEM) and a plasticity-damage model is chosen as the numerical method. It should be noted that the FEM has been enriched to better describe discontinuities. The enrichment methods include the extended FEM 22,23 and the embedded discontinuity FEM. 13 These methods are, however, substantially more challenging from the computational and implementation points of view, especially in problems involving multiple cracks. There are also other approaches more suitable to model rock fracture processes, such as the discrete element method. However, these approaches are inherently explicit in their time discretization, which makes their usage practically impossible, due to the conditional stability of the explicit time marching, in thermo-mechanical problems requiring small elements to adequately model the rock mesostructure and long heating times. Two examples of such approaches are the thermo-mechanical study using peridynamics by Wang et al. 21 and the hybrid finite element-discrete element method by Joulin et al. 12 In both works, a time step of the order of 1 ns is dictated by the method.
In the present modeling approach, a temporally unconditional implicit staggered scheme to solve the global thermomechanical problem under intensive thermal loading is developed. The rock material is described as a heterogeneous damaging (visco)plastic material with temperature dependent thermo-mechanical material properties. The main novelty is the application of polygonal finite elements in modeling thermo-mechanical rock fracture processes for the first time. Polygonal finite elements have been used in general fracture analyses by, for example, Huynh et al. 24,25 and by Khoei et al. 26 Saksala and Jabareen 27 showed that the polygonal finite elements perform better than standard finite elements in modeling failure processes of heterogeneous rock. The modeling approach is verified, qualitatively validated, and then applied in novel simulations of sauna stone degradation as well as thermal comminution of rocks by rapid heating-forced cooling cycles.

THEORY
This section presents the theory of the computational routines for solving the thermo-mechanic problem related to the simulations of rock failure due to thermal loading and the mechanical uniaxial compression tests. First, the damage-(visco)plastic material model for rock is presented. Second, a staggered implicit scheme for solving the thermo-mechanical problem is outlined. Third, an explicit time marching method for solving the mechanical tests is outlined. Finally, the theory of the polygonal finite elements and the method to describe the rock heterogeneity are described. The theory is based on the small deformation framework, justified by the brittle nature of rock fracture under normal laboratory conditions, enabling the additive split of the total strain = e + vp + (1) into elastic, viscoplastic and thermal parts.

Damage-(visco)plasticity model for rock
The purpose of the present study is not to introduce a new constitutive model for rock but to predict numerically the damage in rock due to thermal loading and in the basic uniaxial compression test. Therefore, the simplest possible failure models have been chosen, that is, the Mohr-Coulomb (MC) and Rankine models. The MC criterion-based viscoplasticity model in the 2D case is specified by: where , , and are the components of the stress tensor, , c0 is the initial compressive strength, and are the internal friction and dilation angles respectively, and MC is the viscoplastic potential accounting for non-associated flow. Furthermore, MC anḋM C are the constant viscosity modulus and the rate of the plastic increment, respectively. In the perfectly plastic case MC = 0. Finally, Equation (5) is the consistency conditions meaning that the present formulation of viscoplasticity is the Wang's viscoplastic consistency format. 28 The Rankine model, used as the tensile cut-off, is written similarly by: where t0 is the intact tensile strength while the meaning of the rest of the symbols is analogous to those in the MC model. The purpose for writing these models in the xy-stress space, instead of the principal stress space where their expressions are extremely simple, is to avoid the transformations formulae between the principal and the xy-coordinate systems when deriving the tangent stiffness matrix below. It should also be noted that the viscoplastic consistency format allows to use the robust stress return mapping algorithms of computational plasticity. 13,28 The damage part of the model employs separate scalar damage variables in tension, t , and compression, c . As the damage is driven by (visco)plastic strain, the stress states leading to plastic flow are the ones that induce the damage as well. Hence, no separate damage loading functions are needed in this formulation. The specific model components are: = t0 ℎ ∕ Ic , = c0 ℎ ∕ IIc (10) where parameters A t and A c control the final value of the damage variables. The parameters β t and β c , which control the initial slope and the amount of damage dissipation, are defined by the fracture energies G Ic and G IIc , and h e is a characteristic length of a finite element. The equivalent viscoplastic strain in tension, vp eqvt , is defined, in the rate form, as the trace of the viscoplastic strain rate tensor,̇v p , using the Macauley brackets so that tensile damage evolution occurs only if the volumetric viscoplastic principal strains are positive. The equivalent viscoplastic strain in compression, vp eqvc , being defined with the deviatoric part,̇v p , oḟv p , is similar to that of the J 2 -plasticity for metals. Moreover, Equation (12) is the Koiter's rule for bi-surface plasticity witḣR,̇M C being the viscoplastic multipliers in tension and in compression, respectively. Moreover, the colon in (11) denotes the double contraction operator for tensors, that is, ∶ = . The final model components are the nominal-effective stress relation, which specifies how damage variables operate on the stress, and the stiffness recovery scheme applied upon stress reversal. Here, the relations used by Lee and Fenves 29 and Lubliner et al. 30 are adopted by: wherēis the nominal stress, that is, the one returned on the failure surface when the trial stress violates the criteria. Moreover, t and c are stiffness recovery functions depending on the principal stresses,̄, of the nominal stress. Finally, parameters t and c control the degree of recovery. They are set here as t = 0 and c = 1, which means that t = 1 and c = (̄). This choice means that unilateral effect is applied only upon load reversal from tension to compression. This is certainly a realistic scheme for rock since a fully developed shear band, modeled here by the compressive damage variable, cannot bear any tensile loading. The damage and plasticity parts of the model are combined in the effective stress space, 31 which means that the plasticity and damage computations can be performed separately so that the stress return mapping is first performed in the effective stress space. Then, the damage variables are updated and, finally, the nominal stress is calculated by (13). Moreover, this formulation poses no extra restrictions on the model parameters. 31

Consistent tangent stiffness tensor for perfectly plastic-damage isothermal case
Next, the consistent tangent stiffness operator is derived for the mechanical part in the perfectly plastic case, that is, the viscosity is ignored, to be used with the staggered implicit scheme. This means that during the solution of the mechanical part, the temperature is kept constant. Therefore, the differentiation is carried out with respect to mechanical variables only. Under these conditions, the starting point of the derivation is the nominal-effective stress relation, which is perturbed slightly (variated): where the particular derivatives can be readily obtained from (9) and (11). Term c in (16) is more complicated involving the derivatives of the principal stresses. Starting from (14) and (15), while denoting (̄) = ∑ 3 =1 ⟨̄⟩/ ∑ 3 = 1 |̄| = <> ∕ || , one gets: where 12 is a vector containing the principal stresses. It should be noted that due to the choice of the yield criteria, that is, the Rankine and MC criteria, the gradient of the first principal stress must be calculated in any case so that the increase in computational labor is not an issue.
Next, the variation of the plastic strain tensor needs to be expressed in terms of the variation of total strain. The starting point for this is the Koiter's rule (12) for bi-surface plasticity, which is variated as follows: Now, the variations of the nominal stress and the plastic increments arē where E is the elasticity tensor, and the variation of thermal strain is zero due to the staggered solution scheme. Moreover, the consistency conditions in (5) and (8) have been exploited to solve for the plastic increments (variations in this context) MC and R . With these results in hand, and after some tensor algebra, the variation of plastic strain can be written in terms of the variation of total strain: where II is the fourth order unit tensor. Now, substituting (22) and (20) into (16),̄from (20) into (18) and that in turn into (16), and doing again some tensor algebra gives the final form of the elasto-plastic-damage consistent tangent tensor with epd = + (1 − c ) t̄⊗ UC ∶ During unloading, the tangent stiffness is . Moreover, when only one yield criterion is violated, the evolving parts related to the other criterion are set to zero, for example, Δ MC = 0, c = 0, and d = 0 when only the Rankine criterion is violated. Moreover, the scheme accounting for the stiffness recovery can be neglected by setting c = 1 and UC = 0. In view of Equations (23) and (24), the final form of the tangent stiffness is quite complex even in this isothermal case.

Staggered implicit solution scheme for the governing thermo-mechanical equations
The weak form of the equations governing the thermo-mechanical problem consists of the equations of motion and heat balance, which are written as̈= where ρ and c are the density and the specific heat capacity of the material, is the rate of change of temperature,̈is the acceleration vector, b is the volume force, q is the heat flux vector related to temperature gradient ∇θ and the conductivity k by the Fourier's law, and mech expresses the mechanical heat production through dissipation and strain rate. This last term is set to zero, that is, mech ≡ 0, due to the sheer dominance of the external heat influx. More precisely, the temperature rise due to thermo-elastic and thermo-plastic effects are of order 0.1 • C and 2 • C, respectively, due to adiabatic heat generation during rock fracture under uniaxial compression, 15 while that due to external heating is hundreds of degrees.
The finite element discretized form of the Equations (27) and (28) where u is the nodal displacement vector, θ is the nodal temperature vector, and ext is the external force vector. Moreover, int is the internal force vector, is mass matrix (lumped by row sum technique), is capacitance matrix, is the conductivity matrix, h and h are the contributions from boundary convection (such as natural and forced cooling), and is the vector of thermal loading. These are defined by where A is the standard finite element assembly operator, and B e is the kinematic matrix (mapping the nodal displacement into element strains). Furthermore, e and e are the temperature and displacement interpolation matrices, k is the conductivity, q n is the normal component of the external heat flux, Q int in the internal heat generation, and e is the gradient of e . Finally, h and are the convection coefficient and the ambient temperature. The temperature dependence of the material properties is defined later. Moreover, Equation (32) indicates that the specific heat and conductance depend on damage (to be defined later).
In the present application involving thermal loading, the inertia effects are negligible. For this reason, Equation (30) is solved as a quasi-static problem (̈≡ 0), while the heat Equation (29) is solved with an implicit time marching. Applying the backward Euler schemė= ( +Δ − )∕Δ to Equation (29), written at time + Δ , the Newton-Raphson scheme to solve for +Δ reads 32 where n is the iteration counter, and the nonlinear parts of the tangent stiffness matrix arẽ The corresponding scheme for solving the displacement increment with the quasi-static version of Equation (30) is written as Where −1 epd is the tangent stiffness matrix in Equation (26). Now, the staggered implicit scheme is as follows. First, the Equation (35) is solved for the new nodal temperature while freezing (keeping constant) the displacement field. Then, Equation (39) is solved for the new displacement field while freezing the temperature field. It should be mentioned that in the present case where the internal heat generation due to structural and dissipation effects is zero at the material point level, the adiabatic and isothermal split approaches are identical. Despite this aspect of the present approach, there is a two-way influence between the mechanical and thermal parts as the heat Equation (29) depends on material damage, which in turn depends on displacement.

Explicit dynamics approach to simulate the mechanical tests
The mechanical uniaxial tests on heat treated and intact numerical rock samples are carried out solving the equation of motion (30) with explicit time marching. The explicit modified Euler method 33 is chosen for this end. Accordingly, the system response is calculated as where ,̇and̈are the global displacement, velocity and acceleration vectors respectively at time t, and the rest of the symbols are as above. This resort to explicit time integration is due to the convergence problems of the implicit method in Equations (39) and (40) with the tangent stiffness operator (26). The aim is to simulate the uniaxial compression test until the failure mode is fully developed and, simultaneously, the load bearing capacity of the rock sample is fully exhausted. This a very challenging task for an implicit method due to the presence of the unilateral conditions, extremely steep softening response and non-associated flow rule, which renders the tangent stiffness matrix unsymmetric.
F I G U R E 1 Illustration of the triangular areas used in the definition of Wachspress shape function (A), and the triangulation of the reference regular polygon with three integration points in each triangle, and the isoparametric mapping to a physical element (B)

Polygonal finite elements
Polygonal finite elements offer some benefits in comparison to the usual triangular and quadrilateral elements. These include, in many cases, greater flexibility in meshing arbitrary geometries, better accuracy in the numerical solution, better description of certain materials, and less locking-prone behavior under volume-preserving deformation. 34 Saksala and Jabareen 27 compared the performance of the polygonal elements based on Wachspress interpolation functions to traditional triangular and quadrilateral elements in softening problems and applied them in numerical modeling of heterogeneous rocks with good results. For this reason, the polygonal elements are chosen here as well.
The finite element formulation based on the Wachspress interpolation functions implemented in Matlab by Talischi et al. 35 is adopted. It is based on the standard isoparametric mapping from a reference element to the physical element, as illustrated in Figure 1. The mathematical expression for a barycentric Wachspress shape function at node i of a reference n-gon reads where A(a, b, c) denotes the signed area of triangle a, b, c ( Figure 1A). The numerical integration scheme is based on a sub-division of the reference polygon into triangles and applying a three-point quadrature for each triangle (resulting 3n integration points for each n-gon), as illustrated in Figure 1B. Polygonal finite element meshes are generated by the PolyMesher code developed by Talischi et al. 36 This code generates 2D Voronoi diagrams (tessellations) consisting of centroidal (or alternatively non-centroidal) Voronoi cells.

NUMERICAL SIMULATIONS
The thermal treatment simulations and the mechanical test simulations are carried out in 2D plane strain conditions. First, however, the material properties of the rock forming minerals and their temperature dependence is specified. Then, the verification and validation simulations are carried out. Finally, the thermo-mechanical simulations concerning uniaxial compression test and heating-cooling cycles are performed.

Rock mineral properties and their temperature dependence
Numerical granitic rock is assumed to consist of three minerals: Quartz (33%), Feldspars (59%) and Biotite (one of the black Micas) (8%). The heterogeneity is described by random clusters of finite elements. These clusters, representing the three rock forming minerals, are assigned with the mineral material properties. The mechanical and thermal properties of the minerals are given in Table 1. These values are taken from Mahabadi, 37 Vázquez et al. 10 and Zhao et al. 38 . The uniaxial intact compressive strength for all minerals is c0 = 2 cos ∕(1 − sin ) 0 = 137 MPa. As the temperature expectedly rises hundreds of degrees during heating, it should be considered in modeling. Due to the lack of data on the constituent minerals, the temperature dependence of uniaxial tensile and compressive strengths is modeled after the experiments for granite as an aggregate of minerals. The data collected for tensile and compressive strengths of many granites by Wang These fits, valid for ∈ [293, 823] K, say thus that when the temperature is 550 • C, the uniaxial compressive and tensile strengths are 70% and 50% of their respective values at the room temperature. As to the Young's modulus, the data collected by Toifl et al. 7 is employed here. Accordingly, Young's modulus of Plagioclase Feldspar and Muscovite Mica does not depend on temperature in this range. Therefore, as Biotite is one of the black Micas, Biotite and Feldspar are taken independent of temperature in their Young's modulus. However, Quartz show nonlinear dependence fitted by: where 0 and 0 are the initial temperature (293 K) and the corresponding Young's modulus of Quartz, respectively. Moreover, and are the semi-minor and semi-major axes of the ellipse (47), which is plot against the data in Figure 2A. The thermal expansion coefficient is the single most important material property with respect to thermal loading. The constituent minerals behave very differently in this respect so that special attention must be payed to correctly fit the data in literature. Again, Quartz is the deviant behaving in a strongly nonlinear manner while Feldspar and Biotite show practically linear temperature dependence within the present range of interest. 7 The Quartz data by Polyakova 39 is approximated by the 6 th order polynomial written along with the linear fits for Feldspar and Biotite as:  Figure 2B. The reason for the nonlinear behavior of Quartz is the α-β-transition at 573 • C. This phase change is ignored in the present study. The specific heat capacity and the thermal conductance depend also on temperature, especially those of Quartz. The thermal conductance of Feldspar and Biotite is set temperature independent based on the experiments collected by Toifl et al. 7 The data for Quartz therein is well matched with a second order polynomial by This fit is plotted in Figure 2C. Finally, the temperature dependence of the specific heat capacity is approximated by linear fits for each mineral as follows: where the fit for Quartz data (52) is plotted in Figure 2D, and the values at 293 K are those in Table 1 while f bmax = 1140 J∕kgK is valid for both Feldspar and Biotite. It should be noted that when the temperature falls outside the range, the respective parameters have the values obtained at the end points of the range. When a material deteriorates, through microcrack or void initiation, the thermal properties approach those of air. This feature is taken here into account by relations where ( ) and ( ) are the functions defining the temperature dependence of each mineral (when applicable) and c = (̄) as discussed above. Thereby, the thermal conductance and the heat capacity deteriorate in pace with the material damage. Moreover, upon load reversal resulting in void/crack closure, the thermal properties are recovered.

Thermal stresses in a hollow cylinder: verification problem
The first numerical example concerns the verification of the method against the analytical solution of thermal stresses in a hollow cylinder due to prescribed temperatures (Dirichlet boundary condition) at the inner and outer rings. The problem geometry is shown in Figure 3A.  parts as well as temperature dependence of the material properties, to solve the thermo-mechanical problem for stresses and temperature. Due to the relatively low values of the thermal conduction coefficient here, the simulation time is set to T sim = 1E4 s to reach the steady state solution (55)-(57). A time step of 1000 s is used. The simulation results and the polygonal mesh are shown in Figure 3.
As the problem has radial symmetry and the analytical solutions are presented in radial coordinates, the stresses are plotted by extracting them from the nodes closest to x-axis. According to the results in Figure 3, the analytical solutions are very well predicted with the present method. There are few deviations in the radial stress component in Figure 3E. These deviations are due to the fact that the stress components are not solved at nodes but at Gauss points, from which they are interpolated and averaged to the nodes. In any case, it can be concluded that the present approach is verified against the linear elastic case, that is, the implemented numerical method correctly and accurately solves the initial/boundary value problem it was designed to solve. It should also be noted that the staggered scheme is clearly (temporally) unconditionally stable in this uncoupled, linear elastic case as it was possible to use a time step of 1000 s.

Thermal cracking of concentric cylinders: qualitative validation problem
Next, the model is validated qualitatively against experimental crack pattern in concentric cylinders due to heating of the inner cylinder. The experimental background is in the cracking of reinforced concrete structures under thermal stress investigated by Abdalla. 41 This problem has also been used for qualitative validation of thermo-mechanical models based on peridynamics by Wang et al. 21 and hybrid finite-discrete elements code by Joulin et al. 12 Heating of the steel reinforcement inside the concrete cylinder leads to radial cracking of the concrete, as schematically illustrated in Figure 4B. However, the rock material as described in Section 3.1 is applied here instead of concrete. Moreover, the simulation is carried out by heating the steel phase of the model using the volumetric internal heat generation Q int in Equation (34). The reinforcement properties are set as follows: E = 200 GPa, ν = 0.3, ρ = 7800 kg/m 3 , α = 3E-5 K −1 , k = 5 W/mK, c = 700 J/kgK. It should be noted that the thermal expansion coefficient and the thermal properties may not be those of steel but modified ones due to the fact that the outer cylinder is granite rock, not concrete, which has lower thermal expansion coefficient than granite. These modifications are justified by the similarity concrete and rock materials in all respects relevant to this problem. Figure 4. show the simulation results with a heating duration of 100 s and heat power Q int = 2E7 W/m 3 . The initial temperature is 293 K. The temperature, the first principal stress and tensile damage distributions are shown in Figure 4C-E at t = 13.5 s. The temperature rise is about 35 • C in the middle of the reinforcement. This temperature is enough to initiate tensile damage evolution (note the range of the color bar), which displays radial crack-like patterns emanating from the reinforcementrock interface. At the end of heating (t = 100 s), the maximum temperature is about 125 • C, which is slightly higher than the maximum temperature, 100 • C in the experimental study by Abdalla. 41 The final crack pattern, represented by the damage distributions in Figure 4G and 4H, did not change notably from what it was at 13.5 s-only the values grew larger. Figure 4H shows that some compressive damage has also developed. However, all the compressive damage occurred in elements that has tensile damage as well, that is, both the Rankine and MC criterion were violated at those Gauss points. It should also be noted that the effect of material damage is attested in the temperature distribution ( Figure 4F and Figure 5B) as a lateral weak discontinuity over the main crack-like damage formation reaching the outer ring of the rock cylinder. This is a realistic feature, which most of the previous numerical studies lack. Finally, the effect of rock mesostructure is tested. Figure 5 shows the simulation results with two additional rock mesostructures. The damage patterns are essentially the same as that with the mesostructure in Figure 4A, albeit with differing details. In any case, the general trend in this problem is that the temperature rise initiates several tensile cracks at the reinforcement-rock interface while only one of them reaches the outer edge releasing the lateral stress state, thus preventing the rest of the cracks to propagate further. As this seems to happen in the experiment as well, it can be concluded that the present approach correctly predicts the failure mode in this kind of problem, that is, the numerical method is qualitatively validated against an experiment.

Uniaxial compression test under elevated temperatures: qualitative validation problem
Uniaxial compression test at elevated temperatures measures the temperature dependence of rock compressive strength. In this experiment, the rock sample is slowly heated to the desired uniform temperature after which it is cooled down to room temperature before being subjected to the mechanical test. 11 It should be emphasized here that heating is performed so slowly that no thermal stresses would be generated inside an ideal homogeneous solid. However, rocks are heterogeneous, granite particularly so, and the temperature dependence of the constituent minerals even accentuates it, as clearly manifested in the present case of Quartz being the deviant. Therefore, thermal cracking, modeled here as damage growth, is inflicted on the specimen, which manifests as degraded strength and stiffness in the consequent mechanical test. It should also be reminded that for brittle materials, the compressive strength is not a material property per se, but an emerging property contributed by the specimen size, boundary conditions and, most importantly, by the meso-and microstructure (inherent microcracks, grain size and shape as well as heterogeneity) of the specimen.
Here, the uniaxial compression test is performed on two numerical rock samples (2000 polygons) shown in Figure 6 using the explicit scheme in Section 2.4 applying a constant velocity, 0.05 m/s, boundary condition at the upper edge. The numerical samples display different polygon shapes: the NumRock1 has more centroidal Voronoi cells (polygons), while NumRock2 is a Voronoi tessellation generated for almost random seeds (a single iteration is performed on random seeds by the Lloyd's algorithm, 35 thus better representing a natural rock mesostructure. Figure 6 shows the simulation results for the final failure modes and the corresponding average stress-strain curves.
The failure mode realized with NumRock1 displays a typical experimental "shearing along single plane"-failure mode classified by Basu et al. 42 Notably, the tensile damage variable reach also values close to 1 in the failure plane elements. The corresponding stress-strain response exhibits a non-linear pre-peak part even though the model is linear elastic up to failure. This is a realistic feature and represents microcracking due to stiffness heterogeneity at elements where the compressive strength is reached before the specimen loses its load bearing capacity as an aggregate of grains. The compressive strength, 117 MPa is lower than the nominal strength 137 MPa of each mineral.
The second numerical specimen, NumRock2, shows a shear failure mode with two branches initiating at the upper edge where the specimen has a geometric imperfection ( Figure 6D) with a single node in the mesh deviating 0.036 mm downwards from the upper edge. The resulting uniaxial strength, 97 MPa, is 17% lower than that of NumRock1. This feature demonstrates thus the effect of geometric imperfections of the specimen on the measured strength. 43 Indeed, the strength of NumRock2 without the geometry defect is 118 MPa, that is, practically the same as that of NumRock1 albeit with a differing failure mode. As the geometric defect has extremely small dimension, 36 μm, and the mesostructure variation itself caused only a 0.85% difference in strength, this result suggests that the experimental deviation of the compressive strength of different rocks stems mostly from the specimen imperfections, in addition to imperfections in the boundary conditions. It should however be noted that the present approach to model the mesostructure is a very simple one and cannot capture the grain boundary behavior, which surely influences the rock response under uniaxial compression. In any case, further elaboration of this topic is beyond the scope of the present paper.
Next, the samples are subjected to thermal treatment, that is, heating from 20 • C to 300 • C, with a duration of 6000 s and Q int = 1E5 W/m 3 . The results are shown in Figure 7.
The strong heterogeneity of the numerical rock samples has caused substantial thermal damage, as attested in Figure 7. Naturally, most of the damage is of tensile type (note the range of the compressive damage). The geometric defect in NumRock2 has no bearing in this kind of loading type. Finally, the cooled down numerical samples are subjected to uniaxial compression test. The results are shown in Figure 8 (NumRock2 with the geometry imperfection).
Uniaxial compression induces similar failure modes, that is, shear banding with two branches, on both numerical rock samples. However, the initiation spot at the upper edge is different, being the location of the geometric defect for Num-Rock2. The stress-strain responses exhibit more pronounced nonlinear pre-peak parts for the heat-treated cases with the curves deviating from the intact specimen responses at around 50 MPa of the average stress. The stiffness recovery scheme in Equations (14) and (15) was clearly at play here. The resulting compressive strengths ( Figure 8C) 88 MPa for NumRock1 and 83 MPa for NumRock2 mean 25% and 14% reductions compared to the intact strengths. These results are within the experimental scatter, which is disturbingly wide at this temperature. The granite tested by Yin et al. 9 show 35-29% reduction while the data collected by Gautam et al. 5 has mostly 10% reduction at 300 • C. Some granites even become stronger, that is, show negative reduction in UCS, upon heating up to 300 • C. 5,11 In any case, further elaboration of this topic is beyond the scope here. Finally, NumRock1 is heated up to 550 • C in 3h to demonstrate that certain granites may not be able to bear slow heating without disintegration.
The strong temperature dependence of Quartz thermal expansion coefficient is reflected in the damage distributions ( Figure 9A) inflicted on the numerical rock. The tensile damage variable reach values beyond 0.9 at many elements (Gauss points) and remains below 0.1 at many others, which can be interpreted as a disintegration of the specimen into debris of various grain size. This result is not without an experimental witness, as exemplified in Figure 9B, which shows a complete disintegration of a Finnish granite sauna stone after slow heating of about 4 h. 4 The sauna stone of size 50 × 50 × 50 mm disintegrated before reaching 600 • C. This behavior, which is certainly not typical for sauna stones, can be explained by the anomalous behavior of Quartz thermal expansion upon reaching the Curie point (573 • C) combined with a possible poor condition of the granite. In any case, this experiment lends credence to the present model in its present calibration of material and model parameter values.

3.5
Thermal cracking of rock sample due to heating-rapid cooling cycle: application problem Final numerical example concerns an application of the verified and validated method in modeling rock fracture due to heating-forced cooling cycle. In the first example, the numerical rock sample, NumRock1, is heated to a target temperature of 325 • C in 3600 s and then subjected to forced cooling at the left vertical and top edges. Forced cooling by water is modeled as a convection boundary condition with the heat transfer coefficient h w = 14000 W/m 2 K. 20 This example simulates the thermal loading on sauna stones (the top ones in the sauna stove) due to heating and throwing of water at 20 • C. The simulation results with the cooling duration of 5 s are shown in Figure 10.
Rapid cooling by water at the left and top edges generates a tensile stress state, which in turn induces more, mainly tensile, damage close to these edges (compare Figure 10B to Figure 7A). The temperature dropped very fast to the room temperature, as seen in Figure 10C. During computations, it was thus necessary to adjust the time step in order to capture the steep temperature drop and to achieve convergence. It should be reminded that when the problem is nonlinear due to material failure, the overall staggered scheme is not unconditionally stable in time since the Newton-Raphson scheme in the mechanical part does not converge with an arbitrary large residual (imbalance vector).
This example helps to understand why rapid cooling, or quenching, of sauna stones by water inflicts damage on them. However, this specific case exaggerates the phenomenon as the sample disintegrated when heated up to 550 • C. Good quality sauna stones, usually made of Olivine Diabase, Peridotite or even Ceramics, should endure 1-2 years while retaining their original shape without turning into gravel. Nevertheless, all stone types undergo surface damage in use realized as chipping, which is observed in a weekly manner when cleaning the sauna floor.
The second case considers the same rock sample (NumRock1), but this time heating is applied faster, as a surface influx, while the cooling is performed similarly. The application in mind is comminution by thermal shock provided with a rapid heating-forced cooling cycle. The heating and cooling are applied at both vertical edges of the sample. The heating time is 30 s with the external heat flux q n = 1E5 W/m 2 is applied. The results are shown in Figure 11. Rapid surface heating generates a steep lateral thermal gradient ( Figure 11A), which induces substantial thermal stresses exceeding the tensile strength of the minerals. This results in lateral crack-like tensile damage formations that do not reach the edges of the specimen since the stress state therein is compressive. However, the rapid cooling by water reverses the stress state facilitating further propagation to reach the edges, as attested in Figure 11F. This kind of heat treatment thus suggests a potential non-mechanical comminution method. Further elaboration on this topic is, again, outside the scope of the present method development stage of research.

CONCLUDING REMARKS
This paper developed a continuum based numerical scheme to model rock fracture under thermo-mechanical loadings. The rock material model was based on MC and Rankine criteria with separate scalar damage variables to account for the asymmetry of rock in tension and compression. This relatively simple model correctly captured the experimentally observed rock failure modes in the uniaxial compression and the thermally loaded hollow cylinder. Moreover, the predicted thermal weakening effect at 300 • C on the compressive strength of granite was within the experimental bounds, which, however, are exceptionally wide in this case. Furthermore, the simulation of compressive test on a specimen with an extremely small geometric imperfection, which lowered the specimen strength 18% from that of geometrically perfect specimen, suggests that experimental variation in the compressive strength of rocks is mostly due to these imperfections. Due to the sheer dominance of the external heat influx, it was possible to neglect the heat generation due to structural and dissipation effects, which rendered the underlying thermo-mechanical problem uncoupled. Thereby, the staggered approach, based on the backward Euler scheme, to solve the global equations was unconditionally stable in time for the linear elastic case, which is an advantage over the particle methods, which are based on explicit time integration and are thus temporally conditionally stable. These methods are, however, superior to the present method in fracture description. In any case, the present approach is capable of solving thermally induced rock fracture problems of wide range of heating times.
It should also be mentioned that the effect of material deterioration on the thermal properties was accommodated in the model, which added to the reality of the method through weakly discontinuous temperature field over the crack or damaged elements in the present case.
The rock heterogeneity was described as random clusters of polygonal finite elements representing the granite constituent minerals. The temperature dependence of the mineral thermo-mechanical properties was carefully modeled, which resulted in even more pronounced heterogeneity of the material upon temperature rise due the deviant behavior of Quartz in comparison to the Feldspar and Biotite. The α-β-transition of Quartz at 573 • C was ignored, which restricted the applicability of the method to temperatures below 573 • C. Extension to include the α-β-transition is thus a topic of further development of the method.
The method was finally applied in simulations of sauna stones deterioration under slow heating-rapid forced cooling. The model calibrated for granite predicted notable damage after heating the numerical rock up to 300 • C and then cooling it rapidly by water. Heating the sample slowly to 550 • C lead to disintegration of the sample into gravel, which has also an experimental witness. In the second case, the rectangular sample was heated at both vertical edges with a surface flux to 300 • C in 30 s and then cooled fast by water at the same edges. This thermal shock lead to fragmentation of the sample by multiple lateral cracks, suggesting thus a comminution method for granite by rapid heating-cooling cycles.

A C K N O W L E D G E M E N T S
This research was funded by the Academy of Finland under Grant number 298345.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.