Numerical Modelling of Microwave Heating Assisted Rock Fracture

This paper presents a numerical study on the effects of microwave irradiation on the mechanical properties of hard rock. More specifically, the weakening effect of microwave heating induced damage on the uniaxial compressive and tensile strength of granite-like rock is numerically evaluated. Rock fracture is modelled by means of a damage-viscoplasticity model with separate damage variables for tensile and compressive failure types. We develop a global solution strategy where the electromagnetic problem is solved first separately in COMSOL multiphysics software, and then provided into a staggered implicit solution method for the thermo-mechanical problem. The thermal and mechanical parts of the problem are considered as uncoupled due to the dominance of the microwave-induced heat source. The model performance is tested in 2D finite element simulations of heterogeneous numerical rock specimens subjected first to heating in a microwave oven and then to uniaxial compression and tension tests. According to the results, the compressive and tensile strength of rock can be significantly reduced by microwave heating pretreatment.

Cumulative equivalent plastic strain in compression and tension G Ic , G IIc Mode I and II fracture energy Elasticity tensor EM Electro-magnetic problem

Introduction
In the mineral processing industry, the comminution of hard rock such as granite is a process that demands a considerable amount of time and energy (Ek 1986) and accounts for a large percentage of the total energy consumption in the whole mineral recovery process. Moreover, traditional mechanical fracturing of rocks employs equipment that is subjected to wear, which increases the costs of the operations. Therefore, new energy-efficient comminution methods, where rocks and ores are weakened by a non-mechanical pre-treatment, have been developed in the last decades. These methods include magnetic, electrical and thermal agents (Somani et al. 2017;Klein et al. 2018).
In particular, microwave heating applications of rock have been developed since the 1960s (Ford and Pei 1967;Maurer 1968) within the mineral processing industry, with the scope of reducing energy demand during comminution of ores (Teimoori and Hassani 2021). The main advantage of microwave heating, with respect to conventional heating, is that it is a selective process. In fact, different minerals react differently against microwave irradiation, some of them being either absorbent with regards to radiation (e.g. pyrite, hematite, magnetite) or transparent (e.g. calcite, quartz) (Chen et al. 1984;Tinga 1989;Walkiewicz et al. 1991). The energy absorbed by rock through the exposure to microwave radiation is converted to heat. The amount of heat produced by microwaves in a particular rock depends in general on water content, size (Gao et al. 2011), mineral composition (Nicco et al. 2020), and schistosity (He et al. 2019) of the rock. Moreover, microwave radiation conditions such as power level, exposure time and distance from rock have an influence on temperature intensity and distribution (Wang and Forssberg 2005).
Microwave irradiation-induced heat promotes volume expansion and deformation of mineral grains, which become constrained by the adjacent transparent grains. Moreover, due to the different thermal properties of composing minerals, mineral particles will expand differently. This disparity in expansion between different minerals causes stress concentrations. Therefore, subsequent thermal stresses lead to local inter-and transgranular fractures and possibly to macrofailure (Jones et al. 2005). As a result of microwave pretreatment, rock or ore samples become more amenable to mechanical breakage (Fitzgibbon and Veasey 1990;Marland et al. 2000), allowing for energy saving (Hara et al. 2019).
In laboratory testing, microwave heating is performed most commonly through microwave ovens, either with a single mode cavity (only one mode of energy is excited within it) (Teimoori et al. 2019) or a multi-mode one (Hartlieb et al. 2018;Hassani et al. 2011;Lu et al. 2019). An example of a multimode applicator is a domestic microwave oven. Its basic components are a magnetron, a transmission line or waveguide to convey the power, and a closed metal box (cavity). The reflection of microwaves generates a variety of working modes in the resonant microwave chamber. Multimode cavities are the most popular ovens, essentially because they have the advantage of being mechanically simple, cheap, and able to process specimens having different sizes and geometries. In these ovens, non-uniform heating is frequently experienced due to interference of electromagnetic waves propagating inside the cavity which leads to hotspots formation. Their presence constitutes an important drawback in many microwave applications (e.g. in food industry), but it may become an advantage in microwave pretreatment of rocks for breaking purposes, and this phenomenon may be furthermore enhanced by the thermal runaway effect (Jerby et al. 2019).
Although many experimental studies have been carried out, Hartlieb et al. (2016) have pointed out that few authors deal with the numerical investigation of stresses, cracks and damage generation. In numerical studies, three methods for solving the equations that govern the electromagnetic problem have been widely used: finite difference time domain (FDTD), finite integration technique (FIT), and finite element (FE) methods. Jones et al. (2005) simulated thermal stresses generated in a rock particle embedded in a gangue (microwave transparent) mineral, and exposed to the high electric field strength microwave energy. Meisels et al. (2015) studied the propagation and absorption of a microwave beam and subsequent temperature rise and stress in a two-component rock. Toifl et al. (2016) assessed microwave-induced stresses in inhomogeneous hard rocks at a microstructure level through a 3D simulation procedure. Hassani et al. (2016) modelled, using COMSOL Multiphysics software, microwave heating of a block of rock in a cavity, finding temperature profiles that were in close agreement with their experimental results. Wang and Djordjevic (2014) calculated thermal stresses in a two-dimensional circular plate, containing two-phase minerals, exposed to short-pulse microwave energy. Li et al. (2019) used COMSOL Multiphysics software to simulate a thin section of pegmatite exposed to microwave irradiation at different heating intervals. They investigated the stress-strain responses of different minerals of the pegmatite and their implications for rock strength degradation. Finally, some recent papers have shown promising numerical results on ores strength reduction by means of microwave heating for energy saving purposes during comminution (Whittles et al. 2003;Jones et al. 2007) (see also Korman et al. (2015); Komar Kawatra (2016) for the influence of uniaxial compressive strength on energy consumption during comminution).
Despite their genuine merits and advancement in numerical modelling of microwave irradiation effects on rocks and ores, these studies, and most of the relevant literature, to the best of our knowledge, lack the modelling of the whole process of microwave-induced damage on a rock sample and its weakening effect on rock strength. The present study attempts to fulfil this gap of knowledge. More in detail, the aim of this paper is to evaluate numerically, by 2D simulations, the weakening effect of microwave irradiation on the uniaxial compressive and tensile strength of polycrystalline rock specimens inserted in a multimode oven. This task demands solving a three-stage electromagnetic-thermomechanical problem, where a rock specimen is subjected first to microwave heating pretreatment and then to conventional mechanical tests. To this end, the electromagnetic problem is solved first in commercial FE software COMSOL Multiphysics.
Then, a staggering implicit approach is developed to solve the thermo-mechanical problem, where the microwave-induced volumetric power density serves as a thermal loading input source. Rock fracture is modelled by means of a damageviscoplasticity model. Rock heterogeneity, an essential factor in rock fracture, is represented by means of Voronoi tessellations. In the numerical simulations, the effect of temperature dependence of material parameters and variation of the position of horn antenna are tested, along with the influence of different random rock mesostructures is taken into account.

Rock Constitutive Model
The theory of the rock constitutive model is outlined here. The model, following Saksala (2010), consists of a viscoplastic part and damage part. The viscoplastic part describes the stress states leading to rock failure and determines the inelastic deformation. Accommodation of rock strain-rate sensitivity is rendered through viscosity. The damaged part accounts for stiffness and strength degradation by separate isotropic damage variables in compression and tension due to the asymmetry of rock behavior. Finally, we emphasize that the continuum approach based on viscoplasticity and damage mechanics is chosen for the sake of computational efficiency and implementational simplicity reasons.

Bi-Surface Viscoplastic Consistency Model
The viscoplastic part of the constitutive model is formulated with the viscoplastic consistency approach by Wang (1997); Wang et al. (1997). The stress states leading to inelastic strains and damage are indicated by a bi-surface criterion consisting of Drucker-Prager (DP) criterion for compressive (shear) failure and the modified Rankine (MR) criterion as a tensile cut off. The model components for perfectly viscoplastic case (softening is accounted for by the damage part) are (1) where f DP and f MR are the yield surfaces, g DP is the plastic potential for DP surface (for MR case g MR = f MR ), I 1 is the first invariant of the stress tensor , J 2 is the second invariant of the deviatoric part of stress tensor , i is the ith principal stress, ⟨⋅⟩ are the Macaulay brackets, i.e. the positive part operator, DP and k DP are the DP parameters, s DP and s MR are the viscosity moduli (here s DP = s MR = s , nonzero only in the explicit global solution method), c and t are the dynamic cohesion and tensile strength depending on the viscoplastic increments ̇D P and ̇M R , respectively. The DP parameters DP = 2 sin ∕(3 − sin ) and k DP = 6 cos ∕(3 − sin ) are expressed in terms of the friction angle . Finally, the parameters DP and k DP g are similar to the ones in Eq. (1) 1 , but here they depend on the dilatancy angle instead of the friction angle.

Damage Formulation and Combining the Damage and Viscoplastic Parts
The highly asymmetric response of rocks in tension and compression is accounted for by separate scalar damage variables. The damage in both cases is chosen to be driven by the viscoplastic strain. Consequently, the damage part of the model controls both the strength and the stiffness degradation. In both cases, the standard phenomenological isotropic (scalar) damage model is employed with an exponential damage function. Thus, the damage part of the model is defined by equations where c and t are the damage variables in compression (relating to DP criterion) and tension (relating to MR criterion), respectively. Their maximum values are controlled by the parameters A c and A t . Parameters c and t are defined based on the fracture energies G IIc and G Ic , uniaxial compressive strength c and tensile strength t , and by a characteristic length of a finite element h e , which here is the square root of the element area. This approach ensures an objectivity of the energy dissipation irrespective of the mesh refinement (Feenstra and de Borst 1995). The cumulative equivalent viscoplastic strain in compression and tension, vp eqvc and vp eqvt , respectively, are computed within the algorithmic implementation incrementally, and thus defined by the rate of visco-plastic strain tensor ̇ vp (defined in Eq. (2) (2) with the Koiter's rule of bi-surface plasticity and surfaces both active) and its principal values ̇v p i . Furthermore, the colon in Eq. (2) is the double contraction operator for tensors, i.e. ∶ = A ij A ij . Macauley brackets are applied in the expression of ̇v p i so that only positive principal strains drive the tensile damage.
The damage and viscoplastic parts of the model are combined in the effective principal stress space. Within this approach, the viscoplastic and damage processes are considered independent. Thus, the viscoplastic calculations and stress integration are performed first, independently of damage, in the effective stress space. As a consequence, the robust methods of computational plasticity can be used in the stress integration (Grassl and Jirásek 2006). Finally, the nominaleffective stress relationship is based on the models proposed by Lubliner et al. (1989) and Lee and Fenves (1998) where is the thermal strain tensor defined as where is the thermal expansion coefficient, is the temperature, and is the second order identity tensor.

Stress Integration of Rock Constitutive Model
To solve the model in Eq. 1, the cutting plane return mapping is used. When both yield criteria are violated ( f MR > 0 and f DP > 0 ), the consistency conditions (Eq. (1) 5,6 ) require, at the end of each time step, the following equations to be verified: f MR (̄ t+Δt ,̇t +Δt ) = 0 and f DP (̄ t+Δt ,̇t +Δt ) = 0 . The yield surfaces are formulated to depend on the increments Δ t+Δt i , and the conditions to be fulfilled become Now, by expanding these equations by the first term of the vector valued Taylor series, a linear system is obtained for solving the MR and DP viscoplastic increments: where the components of the gradient are obtained by applying the chain rule for the derivatives. As mentioned before, the viscosity moduli s MR , s DP are nonzero only in the explicit global solution method. After solving the viscoplastic increments by Eq. (6), the stress, viscoplastic strain and internal variables are updated in a standard manner . For simplicity, the stress integration is carried out in the principal stress space. This is feasible as the return mapping preserves principal directions for isotropic material. Computations at the local integration point level are finally summarized in Box 1. The computations thus start with the total strain tot from the global solution and end up in the new stress .

Simulation Method for Microwave Heating Induced Failure
The simulation method for the problem of microwave heatinginduced cracking in rock is outlined in this section. First, the electromagnetic-thermomechanical problem is given in the strong form. The electromagnetic and thermomechanical parts are solved independently with a one-way coupling relationship, since, at this preliminary stage of method development, the dependence on temperature of the dielectric properties is neglected, as similarly as in Toifl et al. (2016); Wang et al. (2017). The radiation-induced heat is calculated after solving the electric field, and it is used as an input load in the thermomechanical problem. It appears that the latter reduces into an uncoupled problem due to the highly dominating role of the radiation-induced internal heat over the heat generation due to dissipation effects. Second, the nite element formulation of the strong form is briefly presented. Third, fully explicit (dynamics) and implicit (quasi-static) solution methods based on a staggered approach are developed for solving the nite element discretized version of the thermo-mechanical problem.

Strong form of the Uncoupled Electromagnetic-Thermomechanical Problem
Modelling of rock breakage due to microwave heating requires solving the underlying electromagnetic-thermomechanical problem. The electromagnetic, thermal and mechanical parts of the problem have immensely different time scales (Ottosen and Ristinmaa 2005). By assuming that the electromagnetic cycle times (in the transient formulation of a time-harmonic EM problem) are short, compared to thermal timescales, the problem is broken down into steps. The first step is computing the EM heat losses. For sinusoidal excitations, the EM problem can be solved in the frequency domain to get the cycle-averaged losses. In the second step, these losses become a constant heat source in the time-dependent heat transfer problem. The electromagnetic part of the problem here is modelled as a wave (electromagnetic radiation with frequencies in the range of 2.4 and 30 GHz for microwaves) propagating through air. The rock sample in the oven is considered as a lossy dielectric object that absorbs part of the wave and scatters the rest in different directions. This is solved in COMSOL Multiphysics software, where the Maxwell equation for the electric field strength EM is solved as a time-harmonic problem ( EM (x, y, z, t) = EM (x, y, z)e j t ). Mathematically, this equation is (Monk 2003 where r is the relative permeability, k 0 the wavenumber, is the angular frequency, 0 is the permittivity of free space, c is the electrical conductivity, j 2 = −1 , and r is the relative permittivity. The real part of the complex relative permittivity, r = � r − j �� r , is the dielectric constant ′ r , and ′′ r is the imaginary part, which accounts for dielectric losses. For air, = 0 , air r = 1 − j … 0 , r = 1 . Here, material properties are considered to be independent of temperature due to the lack of reliable and up to date dielectric properties versus temperature curves. For a 2D analysis, the electric field inside the oven chamber is EM (x, y, t) = E EM (x, y)e j t z . Perfect electric conductor (PEC) boundary condition × EM = is set on the cavity walls where is the normal vector on the walls. Dimensions of the waveguide connected to the RF source (magnetron) are such that the waveguide operates with its dominant TE 10 mode (Pozar 2011). An electromagnetic wave loses part of its energy when travelling through a lossy dielectric medium. This electromagnetic energy is converted into thermal energy within the medium. Thus, P v is the time-averaged volumetric power density absorbed by a dielectric in an electro-magnetic field, calculated in this case as (Haus and Melcher 1989;Jackson 1998) The power density term realizes the one-way coupling of the EM field with the temperature field in the balance of energy equation.
After solving the electromagnetic part of the problem, the strong form of the global balance equations for the thermomechanical problem can be written as where is the material density, is the body force (per unit mass) vector, and c t ( ) is the specific heat. The constitutive equation for the heat flow, which relates the heat flux to the temperature, is Fourier's law = −k( )∇ , where k( ) is the conductivity. Finally, Q int is the thermo-mechanical coupling term that accounts for the mechanical heating due to thermoelasticity and thermo-plasticity in the material. Now, this coupled thermo-mechanical problem can be considerably simplified due to the nature of the microwave heating-induced cracking of rock. Specifically, the order of magnitude of the temperature rise due to microwave heating in this study is hundreds of degrees centigrade. Therefore, the thermo-elastic effects in the bulk continuum material can be ignored as their order of magnitude is about 0.1 • C (see Ottosen and Ristinmaa (2005)). Moreover, it is noted that the temperature rise in a fully developed shear band during a typical laboratory compression test on rock is only a few degrees centigrade (Saksala et al. 2019). Hence, Q int ≡ 0 hereafter. Therefore, the solution of the problem defined by Eq. 9 becomes considerably simpler as the material point level solution of temperature rise can be neglected.

Finite Element Discretized form for the Uncoupled Thermo-Mechanical Problem
The finite element formulation for the heat balance equation becomes, as in Ottosen and Ristinmaa (Ottosen and Ristinmaa 2005), where is the nodal temperature vector, , and are the capacity matrix, the conductivity matrix and the external force, respectively, defined by the following expressions: where c t ( ) and k( ) are the temperature dependent specific heat and conductivity, which depend on temperature, is the temperature interpolation vector, P v is the volumetric power density in the solid induced by microwave heating, and e is the gradient of the temperature interpolation vector. Finally, A is the standard finite element assembly operator.
The mechanical problem is governed by the finite element discretized equation of motion where is the lumped mass matrix, ext is the external force vector, e u is the kinematic matrix, and is the stress vector defined in Eq. (3).

Solution Methods for the Global Thermo-Mechanical Problem
An isothermal split is then applied to the uncoupled thermo-mechanical problem so that the thermal and mechanical parts are solved independently. An implicit-implicit quasi-static scheme is thus developed using the staggered approach to solve the global thermomechanical part of the problem. Here, the mechanical part is solved as a quasi-static problem, which is justified by the negligible inertia effects during microwave heating. Staggered schemes are treated more comprehensively in Ottosen and Ristinmaa (2005); Ngo et al. (2014); Felippa and Park (1980); Martins et al. (2017). For modelling uniaxial compression and tension tests that are to be performed after solving the microwave heating problem, an explicit time marching scheme is adopted.

Implicit-Implicit Staggered Scheme for the Thermo-Mechanical Problem
The thermal part of the problem is solved in the same way as in Pressacco and Saksala (2020), so only the essential steps of the calculations will be reported here. The backward Euler scheme, ̇n +1 = ( n+1 − n )∕Δt , is used in the FE discretized heat equation (Eq. (10)), evaluated at time step n + 1 . The linearization of the equation in n+1 leads to an algebraic equation system (Ottosen and Ristinmaa 2005) where i−1 ,tan is the tangent stiffness matrix for the thermal problem, Δ i is the temperature increment for the current iteration i, i−1 is the residual vector. The superscript i − 1 refers to the previously converged iteration and the subscript n refers to the previous time step. After reaching convergence for the thermal problem, the temperature is kept fixed, and the solution of the mechanical part of the problem can start.
An implicit Euler scheme is used also for solving the mechanical part. To derive the material tangent stiffness matrix tan , Eq. (3) for the stress is written in its rate form, denoting = (1 − t )(1 − c ) , as since ̇ = 0 , and ̇ = due to the isothermal split (the temperature is kept fixed during the solution of the mechanical part).
Material tangent stiffness matrix is derived using the consistency condition in a standard manner. Since the stress return operations here are performed in the principal coordinates system, it is necessary to perform a transformation of coordinates in order to obtain the material tangent stiffness matrix in the xyz coordinate system. After some algebraic manipulations, Eq. (14) can be written as where is the transformation matrix and ′ tan is the material tangent stiffness tensor (principal stress space): where G −1 ij are the components of the inverse of in Eq. (6). The global equation for solving the displacement increment Δ i thus becomes: where i−1 u,tan is the assembled global tangent stiffness matrix and i−1 u = − i−1 int is the residual vector, with i−1 int being the internal force vector. Finally, Fig. 1 shows the flow of the staggered implicit solution procedure.

Explicit Scheme for the Uniaxial Compression and Tension Tests
The explicit modified Euler time integration scheme is chosen for this problem. It solves the equation of motion (Eq. (12)) for the acceleration as then the response is predicted by Hahn (1991) As all explicit time integrators, this scheme is conditionally stable with respect to time step Δt . The time step is restricted by the Courant limit and can be easily calculated.
The explicit time integration is used here due to the convergence problems of the implicit scheme illustrated in Eq. (17), with the tangent stiffness operator defined in Eq. (16). 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 lost. This is a challenging task for the implicit method due to the extremely steep softening response and non-associated flow rule, which renders the tangent stiffness matrix unsymmetric. Moreover, the effective stress space formulation separates the plasticity and damage computations and may thus be detrimental for the stability and convergence since the strong softening is due to the uncoupled damage part of the model: the coupling comes only upon deriving the tangent stiffness operator. Furthermore, the tangent stiffness in Eq.
(16) has one term involving the Macauley brackets. This term may switch on and o (zero and nonzero) during uniaxial compression test simulation since the material is heterogeneous, which means that there are local tensile stresses exceeding the tensile strength, as shown later in the numerical examples. Finally, it should be noted that steep post-peak softening response requires decreasing the time step significantly (from the pre-peak stage) to capture the response so that the computational economy due to the unconditionality of the implicit method is, in comparison to the explicit method, lost to a considerable extent.

Numerical Simulation: Results and Discussion
This section displays the numerical simulations of microwave heating and subsequent mechanical tests. The effects of variation of certain model parameters, along with the mesostructure will be tested. This being a preliminary model development, all simulations are performed in 2D case under plane strain conditions. However, it should be admitted that the 2D results are not fully realistic and more quantitative predictions require full 3D analysis. Moreover, due to the qualitative aspect of this analysis, proper validation of the numerical simulations will be postponed in future studies of the method.

Rock Mesostructure Description
Material heterogeneity (i.e. mineral grains with their different size, shape distribution and material properties, and potential preexisting microcracks and voids) largely influences the mechanical response of rock during both thermal and mechanical loading. Therefore, material heterogeneity is an essential feature that needs to be considered when studying fracture behavior of rocks. Since many rocks (like granite for example) display polygonal grain texture, the heterogeneity is included here explicitly by representing mineral mesostructure by a Voronoi tessellation of convex polygons. Then, the mineral constituents of rock are described by random clusters of finite elements in the mesh. In this study, the Matlab code PolyMesher by Talischi et al. (2012) is used for the creation of the grain texture. This code generates centroidal 2D Voronoi tessellations from either random or a user-specified sets of seed points. The numerical rock samples generated with this method (with dimensions 25 × 50mm , Fig. 2), are composed of three minerals (i.e. 30% quartz, 55% feldspar and 15% biotite), to reproduce a granite-like rock (Haldar and Tišljar 2014). A certain percentage of the grains is randomly selected and assigned with a specific mineral with their respective mechanical and thermal properties.
The minerals and their percentages in the specimen are selected in order to represent not any specific real rock, but a generic granite-like rock. Mineral properties values (for intact rock at room temperature) are given in Table 1. Dielectric properties can be found in Zheng et al. (2020), and they are kept constant during the solution of the electromagnetic problem for simplicity. Mechanical properties for the different minerals are taken mainly from Mahabadi (2012), whereas the values for density, thermal conductivity and specific heat are taken from Schön (2011), Clauser and Huenges (2013), and Waples and Waples (2004), respectively. Finally, the viscosity (used exclusively in compression and tension test simulation with the explicit solution algorithm) is set to 0.005 MPa ⋅ s/m , which is large enough to stabilize the computations and small enough not to cause significant strain rate effects at low-velocity loading.
Temperature dependency of material parameters is taken into account, since temperature rise due to microwave irradiation is hundreds of degrees. It is remarked that possible phase changes, such as --transition of Quartz at = 573 • C , are ignored. In the present study, thermal conductivity, specific heat, Young's modulus, tensile and shear strengths are assumed to depend linearly on temperature in the temperature range from 20 to 500 • C-an assumption justified for granitic concrete (Tufail et al. 2017) and granites (Homand-Etienne and Houpert 1989;Vosteen 2003). Mathematically, this linear dependence is defined as However, the temperature dependence of the thermal expansion coefficient of granite (Heuze 1983) is nonlinear and can be approximated by the second order polynomial with where 0 is the value of the coefficient of thermal expansion of minerals at initial temperature ( 20 • C).

Uniaxial Compression Tests on Intact Numerical Rock Samples
Compression tests are carried out on intact numerical rock samples first to find their failure modes and compressive strengths. These will be the reference cases to which the microwave pretreated samples will be compared. During mechanical tests, the sample is simply supported at the bottom, and a constant velocity of 0.1 m/s is applied at the top side, resulting in a strain rate of 2 s -1 . This loading rate, at which the DIF (dynamic increase factor) is still very modest (less than 1.5) (Zhang and Zhao 2014), is chosen since it is within the range of strain rates in drop weight test machines representing thus the low-velocity impact comminution. The explicit time integration in Eqs. 18-20 is used here. The results of these simulations are shown in Figs. 3 and 4 . Figure 3 shows the tensile and compressive damage variables ( t and c , see Eq. (2)) for Mesostructure 1 at failure. The predicted final failure mode is the experimentally observed shearing along two planes failure mode with some long axial tensile cracks emanating from the main shear failure plane (Basu et al. 2013). The stress-strain response displays a volumetric strain in accordance with the typical experimental behavior (Mardalizad et al. 2018) which is first compact and then dilatant (Jaeger et al. 2007). Based on this simulation, the present damage-viscoplasticity predicts the salient features, i.e. the correct failure mode and the correct stress-strain response, of rock under uniaxial compression.
Shear failure modes, with differing details, can be observed also for the other two intact numerical samples (Mesostructure 2 and 3 in Fig. 2). The average axial stressstrain curve shows variation in compressive strengths values, with a maximum difference between Mesostructure 1 ( 127 MPa ) and Mesostructure 3 ( 145 MPa ). Therefore, the present mesoscopic modelling approach can reproduce, to a certain extend, the experiment deviation in compressive strength of rock (Tang and Hudson 2010). However, it is remarked that under dynamic loading the response of brittle materials becomes more deterministic with higher load rates, i.e. increasing loading rates give smaller deviations from the mean value (Denoual and Hild 2000).

Microwave Pretreatment of Numerical Rock Samples
In this section, the effect of microwave heating pretreatment is tested on intact numerical samples with the boundary conditions defined in Sect. 4.1. The effect of the temperature dependence of material properties, horn position, different power inputs, and the mesostructures (Fig. 2) is tested. Temperature-dependent material properties will be used in the microwave heating pretreatment simulations if not otherwise stated. The oven power is set to 60 kW , and the frequency is 2.45 GHz.

Microwave Pretreatment of Numerical Rock Samples: Reference Case
For microwave irradiation pretreatment, the numerical rock sample is inserted in an oven represented by a conductive box connected to a rectangular microwave horn antenna located at the left side of the cavity. The oven setting is shown in Fig. 5. The sample is continuously supported at the bottom during heating. The 2.45 GHz frequency corresponds to the one commonly found in kitchen ovens, whereas the power level is much higher (most commercial ovens have a power between 0.8 and 1.8 kW ). This choice comes from the need to shorten the total computational time. Moreover, Fig. 3 Uniaxial compression test simulation results for intact rock samples (Mesostructure 1). Failure modes are represented as damage patterns and variables (a, b). Axial stress-strain curve (c) the horn structure, which is not rigorously optimized as an antenna, is introduced to direct EM radiation towards the sample. Actually, in view of the relatively low dielectric loss factor of rock-forming minerals (compared to those of food), the kitchen oven power is not enough for comminution. Cavity size is 600 × 600 mm . All the specimen edges are considered insulated, and the effect of the surrounding air on the specimen is neglected. The initial ambient temperature in the sample is 20 • C. The first simulation of microwave heating pretreatment is conducted with Mesostructure 1 (Fig. 2) using the implicit solution scheme of Sect. 3.3.1. This simulation serves as a reference case to which other simulations are compared. The sample is heated for 60 s in total. To solve the electromagnetic problem, rock mesostructure, after being created in Matlab, is exported to commercial FE solver COMSOL Multiphysics software through LiveLink for Matlab module. Then, the polygonal Voronoi cells are meshed with the regular CST elements, as displayed in a detail in Fig. 2c. Finally, after solving the electromagnetic problem in COM-SOL, the nodal values of the produced heat are exported from COMSOL to Matlab, where the solution of the thermomechanical problem is achieved. Figures 6, 7 and 8 show the simulation results. At the end of the pretreatment, after 60 seconds of heating, the maximum temperature reached in the sample is 419 °C, and the average nodal temperature is 230 °C (Fig. 7d). The curve in Fig. 7c shows the microwave-induced heat in the sample as the power density integrated over the sample volume. The predicted tensile damage patterns in Fig. 8c and g show three crack-like zones inside the specimen. At the higher power level, the value of the tensile parameter is close to 1 therein so that these zones can be interpreted macrocracks caused by microwave heating. The compressive damage is negligible here. Moreover, cracks qualitatively similar to those predicted in Fig. 8 were observed in the experiments by Hao et al. (2020). Therein, the authors heated a cuboid sample of magnetic ore in a microwave muff furnace with a lower power of 2 kW but a longer irradiation time of 5 min.

Influence of Temperature Dependence of Material Properties
In this section, the effect of temperature dependence of some material properties (thermal conductivity, specific heat, Young's modulus, tensile strength and cohesion) is tested by repeating the simulation in the previous section (4.3.1) with temperature-independent material properties. The simulation results along with the reference case are shown in Fig. 9 The total internal heat is the same in both cases since the dependence of dielectric properties on temperature is neglected in this study. The final average nodal temperature in the case of temperature-independent properties is 245 • C , around 20 • C higher than the temperature-dependent properties case. However, in the reference case the tensile damage t is concentrated in three inner horizontal cracks (Fig. 9b, c), whereas in the case of temperature-independent properties the damage t is more diffused in the whole specimen and it does not reach the same value as in the reference case (Fig. 9b). Moreover, in the case of temperatureindependent properties compressive damage is not initiated at all.

Influence of Horn Position
In this section, the effect of horn position is tested by repeating the reference simulation of Sect. 4.3.1 with the horn antenna now set on top of the oven chamber. The power in this case is set to 234 kW , almost four times the one in the reference case. This power level is chosen to obtain results comparable to the ones of the reference case in terms of produced internal heat and average nodal temperature in the sample in a feasible amount of time. The simulation results along with the reference case are shown in Figures 10, 11, and 12.
As mentioned before, multimode microwave heating systems commonly present a number of spatially distributed hot and cold spots inside the cavity. In fact, microwave radiation entering a multimode cavity encounters many reflections to form a complex standing wave pattern, controlled by the dimensions of the cavity and the load. Moreover, the presence of different resonant modes leads to a inherently non-uniform electric field distribution and consequently to non-uniform EM heat generation. Non-uniformity of the thermal field is typically observed in large samples of low thermally conductive materials, as heating will be completely controlled by the electrical field distribution (Teimoori and Hassani 2021). Thus, it seems reasonable that a different specimen-cavity setting produces a different thermal field.  Figure 11 shows that for the horn on top case the hotspot appears at the opposite side of the horn location. The total generated heat and the average nodal temperature curves are fairly similar (Fig. 11c, d). At the end of the pretreatment the maximum temperature reached in the lowest areas of the sample, for horn antenna on top case, is 588 • C , and the average nodal temperature is 226 • C (in the reference case the maximum temperature was 419 • C and the average nodal temperature was 230 • C).
In the second case the damage zones are located approximately where the tensile zones and the largest hotspots are (at the bottom and in the middle of the sample). A similar distribution can be observed in Xu et al. (2021). There, diabase rock samples were irradiated for 30 s with a 3kW microwave applicator.

Influence of Mesostructure
To test the effect of heterogeneity, the reference simulation case in Sect. 4.3.1 is carried out with Mesostructure 2 and 3 (Fig. 2). The simulation results for temperature distribution and induced damage patterns are shown in Fig. 13.
It can be seen that in all three cases the results are fairly similar. The average nodal temperature curves, not shown here, deviated on by few degrees between the different numerical rock samples, whereas the temperature distribution plots in Fig. 13 show some differences in the hot spot distribution. The damage patterns display the same behavior for all three specimens, with three main horizontal tensile damage areas set inside the specimen. The compressive damage was again negligible and is not shown here.

Uniaxial Compression and Tension Tests of Microwave Heating Pretreated Numerical Samples
This section presents the simulations of mechanical tests conducted on the microwave pretreated numerical samples. The mechanical test simulations are carried out on numerical samples cooled down to room temperature. Moreover, it is assumed that the thermal stresses and deformations are zero at the initial state of the mechanical tests. However, the damage induced by the microwave heating pretreatment is present. This represents the in-situ comminution case where a considerable time has passed between the microwave heat treatment and the mechanical breakage.

Uniaxial Compression of Heat-Shocked Numerical Samples: Influence of Horn Position
Uniaxial compression tests are performed on the microwave heating pretreated numerical samples as in Sects. 4.3.1 and 4.3.3. The simulation results for failure modes and stressstrain response are shown in Figs. 14 and 15. From the failure modes (Fig. 14) it can be observed that in both cases new cracks propagate from the pre-existing microwave treatment-induced cracks (damage patterns). Moreover, the failure modes display axially propagating tensile cracks and inclined shear cracks. The drop in compressive strength, referring to the intact strength of 127 MPa are 26% and 65% for the horn on the left and horn on top case, respectively. The horn on top case show much higher strength reduction due to a favourable orientation of the microwave heating induced cracks, which effectively weaken the sample in uniaxial compression (Zhang and Zhu 2020). Theoretically thus, the microwave pretreatment significantly reduces the compressive strength. This is a reasonable result as it is known that imperfections, in general, reduce the compressive strength of rock samples (Pressacco and Saksala 2020).
It should be remarked that the tensile damage due to the microwave pretreatment leads to a considerable stiffness degradation attested in the average stress-strain curves in Fig. 15. This is due to the chosen nominal-effective stress relation (3), which ignores the stress (and stiffness via Eq. (16)) recovery effects upon load reversal. Moreover, the pre-peak parts of the stress-strain curves for the pretreated samples are more nonlinear than those of the intact rock response.

Uniaxial Compression of Heat-Shocked Numerical Samples: Influence of Mesostructure
Uniaxial compression tests are conducted on the microwave pretreated numerical samples with Mesostructure 2 and 3. The simulation results for failure modes and stress-strain response are shown in Figs. 16 and 17. The maximum drop in uniaxial compressive strength (38%) is realized for Mesostructure 3 (Fig. 16c). As for the failure modes, they naturally dier both from each other, due to different mesostructures, and from the failure modes of the intact samples, due to microwave induced damage patterns. These strength reduction percentages agree, to some extent, with those obtained by Hao et al. (2020). As mentioned above, the authors heated cuboid samples of magnetic ore in a microwave muff furnace and subjected the treated samples to the dynamic impact comminution test (drop hammer test). The crack patterns (recorded by high-speed camera) due to hammer impact were somewhat similar to the crack patterns predicted here. However, the much higher loading rate therein resulted in multiple macrocrack planes. For samples cooled down at room temperature (after the microwave treatment), the weakening effects varied 20-30% . Therefore, as the crack patterns due to microwave treatments were qualitatively similar as those predicted here, these experimental results lend some validation on the present simulations.

Tension Test of Heat-Shocked Numerical Samples
In this last section, a numerical tension test is carried out on the intact sample with Mesostructure 1 and microwave pretreated sample with the same mesostructure. The simulation results for the intact and heat-shocked cases with a boundary velocity of 0.05 m/s are shown in Fig. 18. The simulation results for the uniaxial tension test at velocity 0.05m/s shown in Fig. 18 demonstrate that the model can correctly predict the experimental transversal splitting mode with a single macrocrack. The resulting tensile strength is about 8 MPa, which is close to the homogenized tensile strength of 8.45 MPa (calculated by simple rule of mixtures). The pretreated sample expectedly fails by transversal extension of the largest microwave-induced crack. The corresponding tensile strength is about 3 MPa, which is 38% of the intact value. Moreover, the stress-strain response shows a lower tangent modulus and a more ductile post-peak response. Finally, the curvy crack path in Fig. 18c is due to heterogeneous rock description. It can then be concluded here that the present type of microwave heating pretreatment would be a feasible way to weaken the rock in tensile loading.

Conclusions
A continuum approach-based method to simulate microwave heating assisted rock fracture developed and tested in this paper. The rock failure was modelled by a damage-(visco) plasticity model based on Rankine and Drucker-Prager yield criteria, while the strength and stiffness degradation was accounted for by separate scalar damage variables in tension and compression. This model, despite being fairly simple, captures the salient features, including the correct failure modes and average stress-strain responses, of heterogeneous rock under uniaxial tension and compression tests, as was shown in the simulations.
The developed staggered implicit scheme to simulate the thermo-mechanical problem due to microwave irradiation takes as an input the electric field and the related power density distribution in the rock sample solved in COM-SOL software. The temperature rise and the consequent thermally induced damage were successfully simulated in the numerical examples. The simulations yield following conclusions with respect to strength reductions and tested wave guide-sample configurations: When the horn (waveguide) is by the side of the specimen, the thermally induced cracks are lateral being thus more favourable to reduce the tensile strength than the compressive strength of the specimen. When the horn is on the top of the sample, the resulting cracks were located at the bottom of the specimen and displayed aligned orientation to the specimen axis. This irradiation-induced crack pattern was more favourable to reduce the compressive strength of the sample.
The maximum reduction in compressive strength was 26% and 65% in the horn by the side and on the top of the specimen cases, respectively, while 63% reduction in tensile strength was achieved in the horn above the sample configuration. These very optimistic figures are contributed by the long heating times ( 60 s ) combined with relatively high source power ( 60 kW ) of microwave irradiation. Moreover, the unilateral conditions were ignored in the adopted damage formulation. These should be included in further studies of the method. Finally, to fully assess the weakening effect of microwave heating pretreatment on compressive strength rock, the present method should be extended to the 3D case. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 17
Simulation results for uniaxial compression test on microwave pretreated samples (Mesostructure 1, 2 and 3). Axial stress-strain curves Fig. 18 Simulation results for tension test on microwave pretreated samples (Mesostructure 1). Failure modes represented as distributions of damage variables for intact (a) and pretreated specimen (b). Axial stress-strain curves for intact and pretreated specimen (c)