A numerical and experimental study on the tensile behavior of plasma shocked granite under dynamic loading

This paper presents a numerical and experimental study on the mechanical behavior of plasma shocked rock. The dynamic tensile behavior of plasma shock treated Balmoral Red granite was studied under dynamic loading using the Brazilian disc test and the Split Hopkinson Pressure Bar device. Different heat shocks were produced on the Brazilian disc samples by moving the plasma gun over the sample at different speeds. Microscopy clearly showed that as the duration of the thermal shock increases, the number of the surface cracks and their complexity increases (quantified here as the fractal dimension of the crack patterns) and the area of the damaged surface grows larger as well. At the highest thermal shock duration of 0.80 seconds the tensile strength of the Brazilian disc sample drops by approximately 20%. In the numerical simulations of the dynamic Brazilian disc test, this decrease in tensile strength was reproduced by modeling the plasma shock induced damage using the embedded discontinuity finite element method. The damage caused by the plasma shock was modeled by two methods, namely by pre-embedded discontinuity populations with zero strength and by assuming that the rock strength is lowered and conform to the Weibull distribution. This paper presents a quantitative assessment of the effects of the heat shock, the surface microstructure and mechanical behavior of the studied rock, and a promising numerical model to account for the pre-existing crack distributions in a rock material.


Introduction 1
Dynamic fracture behavior of rock materials is one of the most important concepts in rock engineering applications such as percussive drilling, rock burst, and rock blasts. Unfortunately, however, creating the exact environment of these applications in the laboratory is difficult due to the large scale of samples and equipment required to carry out the tests. The rock itself as a material further complicates the testing since its behavior depends on many factors, such as chemical composition [25], structure [26], texture [37], and grain size [34] alongside testing condition parameters such as temperature [7,11,35,39,41], confining pressure [14], specimen size, etc. [6,9,13,18,17,32]. Studying the 1 dynamic behavior of rocks becomes even more challenging in environments such as deep drilling for geothermal energy, where the confining pressure and temperature change as the depth changes. Considering the amount of data needed for a complete description of the rock behavior, a model is required, which can take all these variables into account. A substantial amount of effort has already been put into this matter. For example, Saksala [27] simulated the percussive drilling with a triplebutton indenter using a finite elements based approach, where the rock fracture was modeled with a continuum viscoelastic-damage model. Saksala et al. [29] continued their work and created a 3D model where the cracks are embedded as strong displacement discontinuities. From the works similar to this paper, Cai and Kaiser [5] studied the effects of pre-existing cracks on the tensile behavior of rocks.
Not only creating the exact condition for dynamic testing of rock is a difficult task to perform, numerical modeling of these conditions is challenging as well. The major problem is modeling the dynamic crack propagation [28]. In addition, Zhang and Zhao [43] stated that a strong loading rate sensitivity can be observed in the dynamic loading of rock materials. Many approaches have been used over the past years to simulate the rock behavior under dynamic loading condition. Most of these approaches were based on the finite element method (FEM) and discrete or distinct elements method (DEM). In FEM-based models the cracks are usually represented as inelastic localized deformation zones, see refs. [36,38]. Zhu et al. [45] used a FEM-based numerical code of rock failure process analysis (RFPA) to simulate the behavior of Brazilian disc samples of rocks under combined static and dynamic loading conditions. The RFPA code was developed first by Tang [33]. In this code the rock fracture is modeled by a simple rate-independent isotropic damage model, and the Weibull distribution is used to describe the heterogeneity at a mesoscopic scale. Continuum models have the advantage of higher computational efficiency, but the results obtained from these models do not correspond perfectly to the experimental results. The reason for this can be found from the nature of these kind of models, which are rate-independent. On the other hand, DEM works quite well in the dynamic rock fragmentation and fracture because of the discontinuum assumption. The examples of using this method can be found, for example, in refs [16,19,44]. Combining FEM and DEM provides a unique opportunity to have the features of continuum and discontinuum techniques in one code. This approach was tested by Mahabadi et al. [21,22], whose work shows promising results for fragmentation analysis under dynamic loading conditions.
In this work, we apply the modeling technique presented in refs. [27,29] to the numerical simulation of the Brazilian disc test on plasma shocked granite. This method is based on a combination of isotropic rate-dependent damage model and an embedded discontinuity model. The thermal shock induced cracks are accounted for as randomly oriented displacement discontinuities, which are pre-embedded in the ordinary finite elements. In this method, the exposure time of the plasma shock is represented simply by assuming that a longer exposure time corresponds to a larger percentage of finite elements having a pre-embedded discontinuity.
In our previous papers [23,24], we presented a systematic procedure for the evaluation of the effects of microstructure and testing conditions on the dynamic tensile behavior of granitic rocks, in which the properties of rock were modified using heat shocks. The heat shocks were applied using an oxygen-acetylene flame torch. Considering the fact that using a flame torch is not the best solution from the technical point of view, a plasma shock was used to deteriorate the rock microstructure and properties in this work.

Experimental procedure
The tested material in this study was coarse grain Balmoral Red granite with the mean porosity value of 0.44% and the quasi-static compression strength of about 180 MPa. Table 1 summarizes the chemical composition for this granite. The microstructure of the rock is homogenous and shows no texture. Evidently, the mechanical properties are isotropic and the distribution of minerals is homogenous as well. There are four methods to determine the tensile strength of rock materials: dynamic direct tension test, spalling method, semi-circular bending test (SCB), and Brazilian disc method [40]. Dynamic direct tension test can be the most accurate test since the stress state is 1D, but this method has its downsides. The experimental setup of this method is complex. The requirement for aligning the system is strict to avoid bending of the sample during the test. Preparation of the sample is difficult as well because of the complex shape of the specimen and poor machinability of the granite material. Moreover, a proper high strength glue is needed for joining the sample and the pressure bars [40].
The SCB method is a simpler approach compared to the dynamic direct tension test, however, it should be noted that the result obtained with this method is the flexural strength which is different from tensile strength [43]. The spalling method is also used to measure the tensile strength of the rock materials, but the compressive stress wave may affect the sample properties before the tensile wave reaches the sample. Attenuation of the wave in this method is another drawback. To overcome these difficulties, the Brazilian disc (BD) method seems to be a good solution and yet, this method has its own downsides. The sample should break along the loading direction and close to the center of the sample for a valid BD test [15,31]. However, the instrumentation and carrying out the BD test is easier compared to the above mentioned methods. For these reasons, in this work the BD method was employed to measure the tensile strength of the studied material.
The Brazilian Disc [3] samples with a diameter of 40.5mm were cut from a rock plate with a thickness of 21mm. Figure 1 shows the natural surface of the samples and the modified surface for Digital Image Correlation (DIC) analysis. Before and after applying the heat shocks, optical microscopy was used to analyze the surface structure of the rock. In addition, liquid penetrant was applied on the surface of the samples to characterize the surface cracks before applying the heat shock. The samples were tested at the same strain rate but with a different thermal shock. Three different thermal shocks were applied on the surface of the samples by moving the tip of the plasma gun at different speeds over the surface of the samples. The plasma gun with a power of 50 kW at the distance of 6.5 cm from the surface of the samples was passed with the speeds of 50, 75, and 100 mm/s, providing thermal shocks with durations of 0.80, 0.55, and 0.40 seconds. The temperature of the surface of the samples was measured using an infrared thermometer immediately after applying the plasma shock. The recorded temperatures were 100 °C, 88 °C, and 76 °C for 0.80, 0.55, and 0.40 second plasma shock, respectively. After the thermal shocks, the samples were cooled down in air to room temperature. Unfortunately, it was not possible to obtain the heating and cooling rate of the specimen during the thermal shock process. After applying the heat shock, liquid penetrant was reapplied on the surface of the samples to observe the changes in the surface cracks caused by the thermal shocks. The samples were put under a LEICA CLS 150 XE microscope for microscopy analysis. The lights of the microscope were replaced by an ultraviolet (UV) light source. Later on these images were used to obtain the fractal dimension of the surface crack patterns by comparing the images taken before and after the heat shock. To obtain the fractal dimension, the box counting method was used. This method is commonly used when the fractal dimension of an object cannot be calculated by numerical or analytical methods. The accuracy of this method is highly dependent on image resolution, as the pixels are the grids used in the box counting method. For detailed description of this method, the reader is referred to [10].
High strain rate tests were conducted using a compression Split Hopkinson Pressure Bar (SHPB) device. As suggested by International Society of Rock Mechanics (ISRM), in a static BD test the bars should be modified by using special steel loading jaws to provide an arc contact angle of 10˚ at the failure of the disc sample [4]. Moreover, two adhesive papers should be used to wrap the sample on its periphery. However, this design causes some problems during dynamic testing, such as interferences in the wave propagation, which can result in errors in the measured data. Dai et al. [8] showed that if the SHPB system is properly aligned, the sample can be tested without the special jaws and adhesive papers.
The pressure bars were made of AISI 4340 steel with a diameter of 22 mm and the length of 1200 mm. The striker bar of the same material with the length of 120 mm and the impact speed of 10 m/s was used in the tests. Adjustable stanchions with supporting bearings were used to adjust the movement and alignment of the bars. An air gun was used to impact the striker on the incident bar. Three optical IR sender-receiver pairs were used to measure the impact velocity of the striker. These sensors were also used as the triggering signal for the oscilloscope. Two active strain gages were attached to the center of incident and transmitted bars. Kyowa CDV 700A series signal conditioner was used to amplify the signals of the strain gages, which were recorded on a 12-bit 10MSample Yokogawa digital oscilloscope. Detailed information about the instrument can be found for example in ref. [1]. A numerical dispersion correction based on the work of Gorham [12] was employed to correct the changes in the signals caused by the dispersion of the waves as they travel through the bars. A disc of soft deformable copper with the thickness of 0.5 mm was used to increase the rise time of the incident pulse and to improve the dynamic equilibrium. In this case, the striker bar impacts the pulse shaper before the incident bar, resulting in the generation of a nearly non-dispersive ramp pulse reaching the incident bar. This way the dynamic force balance on the specimen is facilitated [8]. Each test was repeated five times with the impact speed of 10 m/s. Figure 2 demonstrates the schematic view of the high strain rate testing setup.
Digital image correlation was used to analyze the rock fracture during the high strain rate test. As the natural rock surface does not provide enough contrast for the correlation algorithm, the surfaces of the samples were coated by a thin layer of white paint, and black speckles were applied on the surface by a permanent marker. This modified surface is shown in Figure 1b. The images were recorded at 25 kfps with the size of 1024*812 pixels. Two Decocool lights were used for illuminating the surface of the samples. The obtained images were analyzed by StrainMaster (DaVis) 3D-DIC software using the subset of 17 pixels and the step size of 7 pixels. The images taken by the microscope are imported to Matlab as RGB images. As the dominant color of the cracks seen in the microscope images is red and the background tends towards the blue color, the threshold was set for a specific amount of red color and a minimum amount of blue color to differentiate the affected area from the background. If a pixel meets the criteria for the thresholding, the pixel is considered as an affected point and is presented by digit "1". Otherwise the pixel is set as "0". This process produces a binary matrix, which in effect forms a map of the affected area. The obtained fractal dimensions are highly dependent on the threshold value, the resolution of the image, and the noise of the background. Considering these factors that affect the fractal dimension when the matrix is formed, the pattern will be checked by naked eye to decide if it represents the affected zone compared with the original image or not. Figure 3 shows how the thermal shock applied by the plasma gun changes the surface of the sample. It is evident from Figure 3a that before applying the heat shock there are only few surface cracks and surface defects visible. After applying the heat shock, the surface of the sample becomes damaged and some parts of the sample surface are removed.  Figure 3 Microscope images a) before and b) after a 0.80 s plasma shock, and the damaged surface and identified cracks c) before and d) after the plasma shock. The arrow indicates the direction of the movement of the plasma gun.

Fractal dimension analysis
The fractal dimension was obtained for 15 samples before and after applying the heat shock, i.e., five samples for each heat shock of 0.40, 0.55, and 0.80 seconds. The average fractal dimension for the 15 non-heat shocked samples is 1.04 ± 0.08. The average values for the heat shocked samples, in turn, are 1.39 ± 0.17, 1.60 ± 0.08, and 1.74 ± 0.05 for the 0.40, 0.55, and 0.80 second heat shocks, respectively. It should be noted that it is not possible to exactly correlate the mechanical properties of the rock material with the fractal dimension because the fractal dimension is obtained from the surface of the samples only. However, there seems to be a correlation between the relative change in the fractal dimension and the mechanical properties of the rock. Table 2 summarizes the calculated fractal dimension values for each test. It is evident that as the duration of the thermal shock increases, the normalized fractal dimension increases as well. This is simply due to the fact that with lower speed the plasma gun stays for a longer time on the surface of the samples, thus more energy is applied on the surface of the samples and higher temperature is achieved. Since granitic rocks are composed of different minerals with different heat capacities and thermal expansion coefficients, increasing the temperature in a short priod of time leads to the chipping of small particles from the surface of the samples and increases the amount of the damaged area, which leads to a higher value for the fractal dimension. As mentioned before, feldspar and quartz play the most important role in this phenomenon and they are the main cause for the thermal crack development [2].

Split Hopkinson Pressure Bar Results
The results of the dynamic BD tests are shown in Figure 4. Figure 4a shows the stress as a function of time for the specimens without any heat shock, while Figures 4b, 4c, and 4d show the results for the specimens' heat shocked for 0.80, 0.55, and 0.40 seconds, respectively. The average tensile strength for the five non-heat shocked samples is 29 ± 3 MPa. The corresponding value for 0.40s thermal shocked samples is 29 ± 3 as well, indicating that the thermal shock had almost no effect on the tensile strength of the tested material. However, as the duration of the thermal shock increases, the value of the tensile strength drops to 25 ± 4 MPa and 23 ± 3 MPa for the 0.55 and 0.80 second thermal shocks, respectively. The decreases of strength by 20% by applying the 0.8 heat shock shows the efficiency of using a plasma gun as the heat source. It has been shown [23] that applying the heat shock using flame torch leads to the decreases of strength by 58% in duration of 60 seconds. Therefore, the 20% decreases in strength of the rock only in duration of 0.8 second shows the effectiveness of this method. By increasing the duration of the heat shock, the strong peak in the stress vs time curves tends to get lower. This observation with the decrease in the tensile strength of the rock can be related to the fact that by increasing the duration of the thermal shock, the damaged volume under the surface increases.  the plasma gun is powerful (50 kW) and the plasma shock itself can provide a great amount of energy in a very short time, there is still a minimum amount of time needed for this energy to affect a reasonable thickness of the rock material. A thermal shock with a duration of 0.55 seconds decreases the strength of the rock from 29 ± 3 MPa to 25 ± 4 MPa (14% decrease). The decreasing trend of the tensile strength continues as the fractal dimension increases and for the point of 1.72 (0.80 s heat shock), the tensile strength reaches the average value of 23 ± 3 MPa (20% decrease).

Material model for the rock
The material model for the rock fracture is based on a combination of isotropic rate-dependent damage model and an embedded discontinuity model. A detailed treatment, validation and mesh sensitivity studies of the model can be found in [28,29]. Therefore, instead of pointless repetition of validation, parameter and mesh sensitivity simulations, only new simulations related to present application are presented here. However, the model is briefly described for the convenience of the reader. Figure 6. Schematic for the three-phase rock material model. The behavior of rock materials under mechanical loading can be divided into three separate phases, as illustrated schematically in Figure 6. The first is the linear elastic phase followed, upon reaching the elastic limit, by a nonlinear pre-peak hardening involving micro-cracking. As the ratedependent peak stress is reached, a macro-crack will develop triggering the third phase accompanied with non-linear softening. A distinctive feature of the present model is that the strain-rate dependency is included in the model component governing the pre-peak non-linear hardening part of the material response. By this choice, the post-peak softening process leading to the formation of a macrocrack is better captured. Physical justification for this approach is that the microcracking process takes a finite time to develop from micro-scale to macro-scale. Moreover, the classical viscoplasticity formulations with a single yield governing the whole material behavior with a single yield surface predict the final failure surface in a too blurred and messy manner, see [30].
An isotropic rate-dependent compliance damage model is chosen for the description of the nonlinear pre-peak hardening process (the first part is governed by linear elasticity). Here, the term "damage" is to be understood in the damage mechanics sense while in its other occurrences with a reference to the damage induced by the mechanical or thermal loading, it refers to the intuitive comprehension of the term, i.e. a material failure of any kind. For present purposes, such a model is characterised by σ (⊗ denotes the tensor or dyadic product) is the positive part of the stress tensor defined by the principal stresses σi, and directions ni, σt is the elastic limit stress, D is the compliance tensor being initially obtained from the initial compliance tensor De = E0 −1 (E0 being the undamaged stiffness tensor), while κ κ  , are the internal variable and its rate, respectively, λ  is the viscodamage increment, and q is the hardening variable. Variable

[ ])
, 0 ( ∞ ∈ µ µ in (1) is related to the classical scalar damage variable as ) 1 /( µ µ ω + = which varies from 0 to 1. Moreover, s h , are the hardening and viscosity moduli, respectively, and g is a parameter controlling the initial slope of the hardening curve. The inverse of the compliance tensor, Ed, is the damaged stiffness modulus. Finally, operator: in (1) denotes the standard double contraction between tensors. The embedded discontinuity part of the model describes the post peak softening process leading to the final failure of the material. It is based on the decomposition of the displacement field into a regular part and the displacement jump due to a crack. When a body in R 3 is discretized with 4-node tetrahedral elements (see Figure 7), the displacement and strain fields can be written as where the displacement jump is denoted by d α , while i N and e i u are the standard interpolation functions and nodal displacements (with summation on repeated indices), respectively. Moreover, HΓ is the Heaviside function, and δΓ is the Dirac delta function. Function ) (x ϕ is defined for each element with an embedded discontinuity in order to obtain the nodal values in that finite element. The purpose of this function is to restrict the effect of the displacement jump within the corresponding finite element so that the essential boundary conditions remain unaffected. For more details of the finite element implementation of Eq. (2), see [29].
A bi-surface discontinuity model accounting for mode I and II fracture types developed in [28] reads: where n and 2 1 , m m are the unit normal and two tangent vectors for the crack plane (see Figure   7), σ is the bulk stress, and σt and σs are the elastic limits in tension and shear, respectively.
Furthermore, κ κ  , are the internal variable and its rate related to the softening law for discontinuity. The softening slope parameter g is defined with the mode I fracture energy GIc as A fixed crack concept is adopted here, i.e., n remains the same after introducing a discontinuity in an element. A crack is embedded perpendicular to the first principal direction into an element, when criterion dlim deqv ε ε ≥ is met. Here, equivalent damage strain (damage strain), while + d ε is the positive part of it (defined analogously as the positive part of the stress above), and εdlim is the limit damage strain.
In order to have a continuous response upon transition from the continuum model to the discontinuity model, the hardening variable q and the parameter controlling the post peak slope, g above, are set as ) , is the dynamic tensile strength, and GIcdyn is the dynamic mode I fracture energy. This energy is needed to account for the correct dissipation of energy under high strain rates since. The final component of the model is an isotropic damage model. It is introduced in order to alleviate the crack spreading and locking problems related to the embedded discontinuity kinematics implemented with the low order finite elements. In the present formulation, it is introduced into the post-peak softening process when almost all, say 80-90 %, of the load bearing ability of the element with a crack is used. Thereby, the final expression for the stress of an element with a crack reads where ω is the isotropic damage parameter. For more details on this part of the model, see [28]. It should be emphasized that this component of the model is included only in order to deal with the poor performance of the low order finite element. As it becomes active at the very late stage of the postpeak softening process, it does not affect the overall material behavior description of the model. In this model, embedded discontinuities (cracks) with any desired orientation n can be introduced into the elements of the mesh before the analysis. The pre-existing crack distributions can thus be easily modeled. However, it should be kept in mind that the model accounts only for failure types that meet the damage strain criterion presented above. Therefore, the model does not account for the failure in purely compressive stress states.

Split Hopkinson pressure bar simulation model
The SHPB test setup with the BD sample of rock is modeled as follows. The incident and transmitted bars are modeled with standard 2-node bar elements. The incident pulse is modeled as an external stress pulse, σi(t), applied to the end of the incident bar. The contacts between the BD sample and the incident and transmitter bars are modeled in a standard manner by imposing contact constraints between the bar end nodes and the edge nodes of the discretized BD sample. The contact constraints are imposed with the forward increment Lagrange multiplier method, and the explicit modified Euler time integrator is employed in solving the response of the system in time [27]. The local element level solution, i.e., integration of Equations (1) and (3), is presented in [28].

Numerical simulations
The model described above is applied here in the simulation of the dynamic BD test on plasma shocked Balmoral red granite. The plasma shock itself and its effect on the rock microstructure are not modeled in this study. Instead, the plasma shock induced damage is accounted for by two different methods. In the first method this damage is assumed to be representable by the pre-embedded discontinuities with zero strength, and in the second method the damage is modeled as reduced and Weibull distributed mechanical properties. The material properties and model parameters used in the simulations here are given in Table 3.
The material property values for E, ν, ρ and GIc in Table 3 are valid for granite in general. Unfortunately, direct tensile test data on Balmoral red granite was not available for the authors. Therefore, the other model parameters, such as εdlim, g , and σt, were calibrated by trial and error via single element model simulations to match the response of Balmoral granite under static and dynamic direct tension tests. The quasi-static indirect (i.e. measured in the Brazilian disc test) tensile strength of Balmoral granite is ≈ 8.2 MPa, which is matched with the values given in Table 3. As to the other parameter values in Table 3, they were chosen as follows. Continuum viscosity s , dynamic fracture energy GIcdyn , and shear surface parameter β were calibrated by trial and error simulations so that the correct dynamic tensile strength, 29 ± 3 MPa, and the post-peak behavior of an intact BD sample made of Balmoral red granite in SHPB experiment was predicted. Viscosity for discontinuity was chosen so that its value does not affect the post peak response significantly but still secures a stable and robust stress integration, for details see [28] and [29]. Elastic shear limit, σs, was simply set as two times the elastic tensile limit. Finally, shear surface parameter n value ½ provides the correct unit (MPa) for shear stress (with unitless β ). The four-node tetrahedral mesh and an example of the crack normal distributions used in the simulations are shown in Figures 8a and c. The initial crack population was generated as follows: First, a zone where the plasma gun induced cracks are located was defined with a width of 15 mm. This area was divided into nine layers, as shown in Figure 8b. Then, the elements in these layers were searched. The number of cracks in each layer is dictated by two crack density percentages p1 and p2. In the first layer, which is the layer directly affected by the plasma gun, the number of elements having an initial crack is Ncr1 = p1Ne1 (Ne1 is the number (amount) of the elements in this layer). The number of cracks in the next layers up to layer 9 is simply NcrN = p2Ncr(N-1) (2 ≤ N ≤ 9). The elements in each layer having a crack are randomly chosen in order to obtain a spatially random distribution of cracks. Finally, oriented embedded discontinuity was assigned for each element (in the mesh) designated as containing a crack. Figure 8c shows an example of the crack population with p1 = 0.8 and p2 = 0.75, resulting in 2471 cracks in total with layer 1 having 703 cracks and layer 9 only 72 cracks. The strength of a discontinuity representing an initial crack is set close to zero.
The orientation of the pre-embedded cracks is biased so that only the directions in the first quadrant are allowed. Namely, the crack normal vector for each pre-embedded discontinuity is generated by nd = [r1, r2, r3] T , where the components ri are random numbers from a uniform distribution between 0 and 1. This orientation may not be the case in reality but serves for these preliminary modeling purposes. It is noted that an unbiased crack orientation scheme was also tested, but it did not provide the desired weakening effect within the present modeling approach that accounts presently only for tensile type of fracture.
The first set of simulations was carried out by varying the percentage p1 while keeping p2 fixed to 0.9. The results with GIcdyn = 5GIc for the final failure modes and tensile stresses (as a function of time) are shown in Figure 9. To facilitate the comparison, the tensile stress vs. time curves are calculated similarly as in the experimental procedure (i.e. by σx(t) = P(t) /(πRH) where R and H are the radius and thickness of the disc, and P = (P1 + P2)/2 is the average of the contact forces between the incident and transmitted bars and the disc). According to the results shown in Figure 9, all the failure modes display the major axial splitting crack but the details differ. The most notable difference is that the secondary radial cracks are less pronounced as the number of cracks increases. As for the tensile stresses shown in Figure 9f, a clear tendency of decreasing tensile strength can be observed upon increasing the crack number. With parameters p1 = 0.8 and p2 = 0.9, the tensile strength is 23 MPa, while without cracks it is 31.5 MPa. In this respect, the present approach seems to be able to capture the rock weakening effect of the microcrack populations. However, more experimental research is needed in order to decide whether this modeling approach correspond to the plasma shock induced damage and cracks.
One more simulation was carried out with crack density percentages p1 = 0.8 and p2 = 0.75. This scheme results in 2463 cracks in total with layer 1 having 703 cracks and layer 9 only 72 cracks. Thereby, this scheme presumably represents more realistically the plasma shock induced crack systems when only one side of the sample is treated with the plasma gun. The results with these parameter values are shown in Figure 10b. Another simple approach to account for the plasma shock induced damage is to assume that the mechanical properties are deteriorated randomly. Here this approach is tested with the assumption that the tensile elastic limit and the Young's modulus conform to the Weibull distribution. The 3parameter cumulative Weibull probability distribution function reads ( ) where mw is the shape parameter, x0 is the scale parameter corresponding to the mean value of x, and xu is the location parameter. The parameter values used in this paper are mw = 3, xu = 0, for both σt and E and x0 = 0.37E0 for E and x0 = 0.2σt0 for σt (where E0 and σt0 are values given in Table 3).
These values for the mean value parameter x0 are taken from Ref. [42]. Therein, experimental results for the mechanical properties of granite are reported after thermal treatment up to 600 °C. At his temperature, the Young's modulus and the tensile strength were reduced to 37 % and 20 %, respectively, from the values at room temperature. As for the elastic limit in shear, it is simply taken as two times the elastic limit in tension.
This distribution is used as follows. After generating a uniformly distributed random data between 0 and 1 (Pr(x) in (5)), corresponding to the number of elements in the plasma shock affected zone of the mesh, and solving for x in (5) for each element, a spatially statistical distribution of a material property is obtained in the damaged zone. The zone of the plasma shock affected elements with this method is the same as with the pre-embedded discontinuities approach. Figure 8d shows the case for the Young's modulus with p1 = 1, p2 = 0.75. Moreover, the limit of the damage strain, εdlim, is also lowered to 20 % of the value given in Table 3. The simulation results with the Weibull approach are shown in Figure 10c-e for three different numerical specimens. According to the results shown in Figure 10, the failure modes display similar overall features with deviations only in details. Most notably, the amount of damage with the initial crack population and Weibull approaches is larger at the bottom (z = 0) surface of the disc. Moreover, the deviations between the results for the different numerical specimen predicted with the Weibull method are negligible. This finding, which means that the brittle materials become more and more deterministic under increasing loading rate, is in line with the theoretical considerations and Monte Carlo simulations carried out with a probabilistic damage model based on Weibull distribution.
As for the tensile stresses shown in Figure 10f, the peak stress of the intact rock is 31.2 MPa, while the maximum stress with the initial crack and Weibull approaches are 26.3 MPa and 26.4 MPa, respectively. This means that the 15.7 % (initial cracks) and 15.4 % (Weibull distribution) weakening in the tensile strength with these methods are very close to the experimental value (14 %) with the longest tested exposure time. In this respect, both of the tested approaches are capable of capturing the rock weakening effect of plasma shock induced damage.
In passing, it is noted that mesh sensitivity study is not included here since it is already done with respect to the present application (albeit without the pre-embedded cracks or the weakened strengths) in [29] for 2D and [28] for 3D case. Therein, the present numerical is approach is shown to predict mesh objective results in uniaxial tension test as well as Brazilian disc test both in 2D and in 3D cases.

Discussion
The results obtained from the SHPB tests show a maximum tensile strength of 29 ± 3 MPa, which matches considerably well with the results of the first approach in the numerical simulation (31.75 MPa). After a 0.40 second heat shock, as mentioned before, the tensile strength does not essentially change and correlates to the situation in the simulation where p1 = 0.2 and p2 = 0.9. A 0.55 second heat shock decreases the strength of the rock to 25 ± 4 MPa. This value is a good match with the simulation condition where p1 = 0.6 and p2 = 0.9. When p1 = 0.8 and p2 = 0.9, the maximum tensile strength obtained from the numerical simulations (23 MPa) matches well with the SHPB result when the rock samples were heat shocked for 0.80 seconds. Comparison of the results of the SHPB tests with the second approach with an assumption that the mechanical properties are deteriorated leads to a different situation. The peak stress of the intact rock in the simulations is 31.2 MPa, while the maximum stress with the initial crack and Weibull approaches are 26.3 MPa and 26.4 MPa, respectively. This means that 15.7 % (initial cracks) and 15.4 % (Weibull distribution) weakening in the tensile strength with these methods are very close to the experimental value (14 %) with the longest used exposure time.
Even though the quantitative analysis provides an easier comparison to understand the behavior of the rock in different testing conditions, it does not provide a full description of the mechanical behavior. For this, one has to pay attention to many details such as the strain distribution during the loading, fragmentation, etc. DIC has proven to be a useful tool for the comparison of the simulation and experimental results. Figure 10 shows the high speed images of a sample during dynamic loading. The displacement in the direction perpendicular to the loading (Y-direction) between adjacent subsets, or ∆ = 0 , are also shown in the original images. Figure 10a shows the displacement over the surface of the sample before the test starts. No displacement can be observed before the image of Figure 10c with the time stamp of 40 µs. This observation indicates that the fracture initiates inside the samples and then propagates to the surface only after the sample has already reached its maximum stress. After this point, the crack propagates with high speed, followed by crushing of the sample at contact points with the loading bars. Unfortunately, the speckle pattern applied on the surface of the samples for the DIC analysis prevents us from analyzing characterizing the actual initiation and propagation (intergranular or transgranular) of the cracks.
Finally, it is instructive to compare the numerically predicted failure mode development to the experimental one shown in Figure 11. Figure 12 shows the simulated development of the norm of displacement jump at different times in the case of the initial crack method.
When comparing the simulated results in Figure 12 to the experimental ones in Figure 11, it can be observed that there is a satisfactory agreement. Even though both sets of results are quantitative, the comparison is not trivial since the DIC images display displacements on the surface of the sample, while the simulation results are presented in terms of the displacement jump norm over the whole volume of the sample. However, qualitatively the experimental results in Figures 11 d and

Concluding Remarks
The surface of the rock material was modified by applying heat shocks on the surface of Brazilian disc samples using a plasma gun. Later on, Surface damage caused by the heat shocks was quantified using optical images to obtain the fractal dimension of the surface crack patterns. The fractal dimension results show the increase in the amount of surface damage as well as the increase in the complexity of the crack network. It is clear that the heat shock affects a larger volume than just the surface and that there are cracks also under the surface, which are not visible with visible or ultraviolet light. However, this method can provide reasonable estimates for the effects of the heat shock on the tensile strength of the studied rocks.
Afterwards, the samples were tested using the Split Hopkinson Pressure Bar device. The tensile strength of the tested material does not show any change for the 0.40 second heat shock, even though the fractal dimension results indicate surface damage also for this heat shock duration. This means that the affected volume was very small and the duration of heat shock was not long enough to transfer the energy deeper into the sample. However, when the duration of the heat shock increases to 0.5 seconds, the tensile strength of the material drops by 14% from its original value. Increasing the duration of the heat shock up to 0.80 second leads to a 20% decrease in the tensile strength of Balmoral Red. Two different approaches were made to numerically simulate the dynamic behavior of the rock material. The first approach is based on the number of initial cracks. The results of the simulations based on this approach match with results obtained from the SHPB tests. The second approach is based on the assumption that the mechanical properties of the rock are degraded and conform to the Weibull distribution. The Weibull parameters can be adjusted with the data on thermally affected rock material properties. As for the tensile stresses, the peak stress of the intact rock is 31.2 MPa, while the maximum stress with the initial crack and Weibull approaches are 26.3 MPa and 26.8 MPa, respectively. This means a 15.7 % (initial cracks) and 14.1 % (Weibull distribution) weakening in the tensile strength, which are very close to the experimental value (14 %) with the longest tested exposure time. In this respect, both of the tested approaches are capable of reproducing the experimental weakening effect of the plasma shock induced damage.
Finally, digital image correlation was used to observe the crack initiation and propagation on the surface of the samples during the high rate tests. Time stamps from the tests were recorded on the images, and the data from the SHPB tests was synchronized with the obtained images. The obtained results from DIC correlate well with the numerical simulations showing the displacement values. However, it should be noted that the DIC results show the displacement on the surface of the samples, while the simulation results include the whole specimen volume.