Numerical modelling of dynamic spalling test on rock with an emphasis on the influence of pre-existing cracks

This article deals with numerical modeling of rock fracture under dynamic tensile loading and the related prediction of dynamic tensile strength. A special emphasis is laid on the influence of pre-existing natural microcrack populations as well as structural (artificial) cracks. For this end, a previously developed 3D continuum viscodamage-embedded discontinuity model is employed in the explicit dynamic finite element simulations of the spalling test. This model is capable of modelling the effect of natural microcracks populations always present in rocks as well as to capture the strain rate hardening effect of quasi-brittle materials. In the numerical simulations of spalling test on Bohus granite, it is shown that the model can predict the pullpack velocity of the free end of the intact rock sample and the effect of structural cracks with a good accuracy. According to the simulations, the effect of microcrack populations, modeled here as pre-embedded discontinuity populations, is weaker than the corresponding effect under quasi-static loading.


Introduction
In many rock engineering applications, such as blasting and percussive drilling, loads are dynamic and transient, i.e. associated with high amplitude and short duration stress pulse. In these applications, information on the mechanical properties of rock material is important in assessing the stability of rock structures. One of the key material properties of quasi-brittle materials, such as rock and concrete, with respect to stability and fracture analyses is the dynamic tensile strength. There are many experimental methods to measure the dynamic tensile strength of rock directly and indirectly, see [1] for an extensive review.
One of the direct tests to measure the dynamic tensile strength of rocks and concrete is the spalling test using the Split Hopkinson Pressure Bar (SHPB) [2]. In this test, a compressive stress wave generated by the SHPB device travels through the incident bar and the cylindrical rock specimen attached to it. The stress wave reflects as a tensile wave at the free end of the rock specimen and, as the tensile strength is only a fraction of the compressive strength, break the specimen by transverse splitting mechanism, i.e. spalling, into two or more pieces. The dynamic tensile strength of the rock specimen can be then calculated, e.g. from the measured velocity of the spalls [3].
In the computer analyses of rock structures, the major challenge is related to the numerical description of rock fracture processes. Many numerical models, based on the finite element method and the discrete element method, have been developed during the last decades for the simulation of rock fracture. For a review article on numerical methods in rock mechanics in general, see [4]. A constitutive model for rock under dynamic loading conditions is an important asset, not only in the numerical analysis of above mentioned rock engineering applications but also in the development and analysis of reliable experimental techniques. The above described spalling test, for example, relies heavily on the assumption that the incident compressive stress wave does not damage the rock material. The validity of this assumption can be investigated by numerical modelling, as done in [5]. Therein, the effect of the rock heterogeneity on the tensile strength was studied in 2D finite element simulations of dynamic spalling using a damage model (RFPA code). According to these simulations, the compressive wave damaged the rock material (its numerical representation to be exact) when numerical heterogeneity was increased and the compressive wave amplitude was high enough.
Another aspect that can be studied by numerical simulation is the influence of preexisting cracks on the dynamic tensile strength. This issue was studied by saadati et al. [6] experimentally and numerically with a probabilistic damage model based on the Weibull distribution implemented in FEM. They obtained 3D simulation results for pullback velocity and dynamic tensile strength with fairly good agreement with the spalling experiments on Bohus granite. However, they did not present the simulated failure modes. Forquin and Erzar [7] applied the same modelling technique in numerical and experimental study of dynamic fragmentation of concrete. They presented the fracture modes in terms of damage variable distributions in the finite element mesh. Finally, the experimental and numerical study on the dynamic failure of sandstone by Das and Shen [8] is mentioned. They implemented the rate dependent Taylor-Chen-Kuszmaul model in Abaqus finite element code and simulated the spalling test on Piles Creek sandstone. The failure modes were presented in terms of damage contours.
All the above mentioned numerical studies fall to category of continuum methods based on finite elements. This class of numerical methods suffer from poor representation of discontinuities, e.g. cracks. Therefore, as FEM has indisputable advantages in genuine continuum problems as well as unsurpassed maturity as a numerical method, possibilities to extend its applicability to the problems involving discontinuities, such as cracks, by various enriching techniques have been searched during the last two and half decades. Two main categories of these methods are the element-based enrichment [9] (Embedded FEM) or the nodal-based enrichment [10] (Extended FEM). While having quite similar accuracy properties in crack propagation simulations, the former method is computationally superior due to the local nature of the enrichment which allows for elimination of the extra degrees of freedom (related to the displacement discontinuity) by static condensation [11].
In the present paper, the embedded discontinuity FEM is chosen for numerical modelling of failure in the spallation test. For this purpose, the combined viscodamageembedded discontinuity model developed by Saksala et al. [12][13] is employed. It will be shown that not only the experimental pull-back velocity and hence the dynamic tensile strength but also the failure modes in the spalling test, in contrast to the above mentioned studies, can be predicted with this model. A special emphasis of the present study is on the influence of the pre-existing structural cracks and natural microcrack populations al- ways present in rocks. This aspect was also studied by Saadati et al. [6] but, as mentioned above, they did not present the predicted failure modes. Moreover, the structural cracks were modelled as zones of elements in the mesh with (almost) zero tensile strength. In the present work, the structural crack geometry is realistically presented and the inherent microcracks are represented as pre-embedded discontinuities with negligible strength.

Constitutive description of rock
The theory of the rock material model is presented here. In this model, the rate-dependent material behavior under dynamic loading is divided into three phases of response illustrated schematically in figure 1. When the elasticity limit is reached at the end of first linear elastic phase, microcracking starts in the second phase and the response remains isotropic nonlinear hardening until the rate-dependent peak stress is reached. Upon reaching the peak stress, a macrocrack is introduced triggering the third non-linear softening phase. A distinctive feature of this model is the rate-dependency, which is incorporated in the pre-peak part of the constitutive description. A physical justification to this choice is that micro-cracking processes are rate-dependent: it takes time, say incubation time [14], for a fracture process to develop from micro-scale to macro-scale. Thereby, higher strain rates lead to higher peak stresses.
In this section, the viscodamage model to describe the pre-peak rate-dependent hardening response of rock is first briefly described. Then, the bi-surface embedded discontinuity model to control the post-peak softening process is sketched. As the model is already presented in the aforementioned references, only the necessary (for comprehension) components are given here.
Viscodamage model for rate-dependent pre-peak hardening description In the classical method to incorporate rate-dependency into a continuum constitutive model, a constant viscosity term is added to static strength (see [15] for example). If the pre-peak and post-peak response is governed by a single yield function, this approach leads to unrealistically wide and messy failure bands [16]. For this reason, Saksala et al. [12] developed a rate-dependent extension of the three-phase damage-embedded discontinuity model by Brancherie and Ibrahimbegovic [17]. In this extension, the rate-dependency is included in the pre-peak nonlinear hardening process (as illustrated in figure 1) described by an isotropic viscodamage model. By restricting the rate-dependency in the part of the model governing the pre-peak response, the post-peak softening process resulting in a crack can be more accurately captured in dynamic applications [12] [13].
The viscodamage model describing the pre-peak hardening phase is formulated analogously to the viscoplastic consistency approach by Wang et al. [18]. Moreover, the compliance damage format is employed here. The main ingredients of such a model arē In equation (1),φ is the viscodamage loading function, σ + is the positive part of the stress tensor, σ, σ t is the elastic limit for tensile stress, E 0 is the Young's modulus of intact rock, D is the compliance (flexibility) tensor, the initial value of which is obtained as the inverse of the intact elasticity tensor by D e = E −1 0 ,κ,κ are the internal variable and its rate, respectively, andq is the hardening variable controlling the evolution of peak stress. In equation (2), variableλ is the viscodamage increment, and variableμ ∈ [0, ∞) represents the additional material flexibility due to damage. An exponential hardening law with a linear viscosity term is spefified for variableq in (3). Therein,h,s are the hardening modulus and constant viscosity moduli, respectively, andḡ is a parameter controlling the initial slope of the hardening curve. Finally, Equation (4) states the Kuhn-Tucker loading-unloading conditions. The stress and other variable updates for this model are given in [19].

Finite element with a discontinuity plane
Let a body Ω ∈ R 3 be discretized with 4-node, linear tetrahedral elements. This element is well suited for explicit analyses involving transient loads and stress wave propagation and facilitates the FE implementation of the embedded discontinuity theory considerably because the displacement jump remains constant (element-wise). Assume now that localization of deformation, through coalescence of the growing micro-cracks, has taken place so that the body is split into two disjoint parts. Thus, some of the finite elements are also split into two parts by a displacement discontinuity Γ d defined by the normal n and the tangent vectors m 1 , m 2 , as illustrated in figure 2a.
With the linear tetrahedral element, the discontinuity surface becomes a plane inside the element. Assuming infinitesimal deformation kinematics, the displacement and the strain fields for this element can be written as In equation (5), α d denotes the displacement jump, N i , u e i are the standard interpolation function and nodal displacement vector for node i, respectively. In equation (6), H Γ d , δ Γ d are the Heaviside function and Dirac's delta function at Γ d , respectively, while function ϕ restricts the effect of the displacement jump within the corresponding finite element so that the essential boundary conditions remain unaffected [20]. It should be noticed that in arriving at (6), the gradient of the displacement jump, ∇α d , is zero for linear tetrahedron. Moreover, the term containing the Dirac's delta function, (6) is valid only when x ∈ Γ d . Elsewhere, i.e. when x ∈ Ω ± e this term is zero so that it can be neglected at the global level of solving the discretized equations of motion.
Function ϕ is chosen so as to make its gradient as parallel as possible to the discontinuity normal. This choice, which makes the element tangent stiffness matrix wellconditioned, is obtained by criterion [21] ∇ϕ = arg max The three possible modes to be examined in this criterion are illustrated in figure 2b. The finite element formulation of the balance of linear momentum leading a simple implementation with the linear tetrahedron can obtained by the EAS (Enhanced Assumed Strains) concept where the enhanced modes are constructed in the strain space, orthogonal to the stress field [22]. This procedure leads to the following expression for the weak traction balance equation where V e is the volume of the element, A e is the area of the discontinuity surface and σ, σ Γ d are the stress in the element (excluding the discontinuity) and at the discontinuity, respectively. With the linear tetrahedron, the integrands in (8) are constants inside the element yielding the strong traction equation where E d = D −1 is the degraded elasticity tensor. The similarity of Equation (9) to the stress expression in plasticity theory based on the decomposition of strain into the elastic and plastic parts suggests that the methods of computational plasticity can be exploited in solving the traction and the displacement jump.

Bi-surface model for solving the displacement jump
A plasticity inspired bi-surface model separating the tensile and shear modes for the solution of the displacement jump and and the traction update is briefly described here. The loading functions, softening laws and the evolution equations of this model are summarized as where t Γ d is the traction vector, σ s , σ t are the elastic stress limits in tension and shear, respectively, κ,κ are the internal variable and its rate related to the softening law for discontinuity, andλ t ,λ s are the crack opening (mode I) and sliding (mode II) increments. The softening slope parameter g is defined with the mode I fracture energy G Ic as g = σ t /G Ic and constant viscosity modulus s is assumed. Finally, n and β are parameters that control the shape of the shear loading function. The displacement jump and other variable updates for this model are given in [19]. Viscosity is introduced into the embedded model for stabilizing purposes only, not in order to incorporate the strain rate dependency (which is done by the continuum viscodamage model). Including viscosity allows for very steep slopes of the softening response, see [12] [13].
Transition from pre-peak hardening to post-peak softening A fixed crack concept is adopted which means that n remains the same after introducing a discontinuity in an element. A discontinuity is embedded perpendicular to the first principal direction into an element, when the following criterion is met: where + d is the positive part of the damage strain d and dlim is a given limit of the equivalent damage strain deqv . In order to secure a continuous response upon transition from the viscodamage model to the discontinuity model, the hardening variable q and the parameter controlling the post peak slope, g, are set as q =q, g = σ tdyn /G Icdyn where σ tdyn is the dynamic tensile strength, and G Icdyn is the dynamic mode I fracture energy.
Thus, the parameters by which the model response is calibrated to match the quasistatic experiments are dlim , σ t ,ḡ and G Ic . In dynamic tensile experiments, the terms to be calibrated are G Icdyn ands.
The final component of the model is the isotropic damage model. It is employed in order to alleviate the locking and spreading problems related to the embedded discontinuity kinematics implemented with the low order finite elements [22]. The method is to introduce the isotropic damaging at the end phase of the post-peak softening process when, say, 80 − 90 percent of the strength is exhausted. Thereby, the damage variable operates on the stress, computed after solving for the displacement jump in Equation (9), by the standard nominal-effictive stress relation (σ = (1−ω)σ). However, as the crack spreading and locking occur primarily in applications with bending effects, this component of the model is not elaborated here in more detail.

Numerical simulations
Since the model is extensively validated in different applications in references [12] [13], only the spallation test simulations are presented in this section. Moreover, the model performance is demonstrated therein at the material point level also. Therefore, instead of pointless repetition of already published results, only new simulations are presented.

Spallation test simulation model
The spallation test, explained in Introduction, based on the SHPB device is modelled by the principle illustrated in figure 3. The compressive stress wave induced by an impacting projectile is simulated as an external stress pulse, σ i (t) obtained from the experiments. The incident bar are modeled with two-node standard bar elements, and the rock sample is meshed with the 4-node tetrahedrons. Finally, the contacts between the bar end and the rock sample are modeled by imposing kinematic contact constraints between the bar end node and the rock mesh surface nodes. These constraints are of form u bar,z −u n,z = b n where u bar,z and u n,z are the axial degrees of freedom of the bar node and a rock contact node n, respectively, and b n is the distance between the bar end and rock boundary node. The contact constraints are imposed with the forward increment Lagrange multiplier method (see [15] for details).
The FE method discretized equations of motion are solved explicitly in time. The modified Euler method is chosen for this end. Accordingly, the response (velocity and displacement) of the system is predicted as follows u t+∆t = u t + ∆tu t+∆t (21) where M is the lumped mass matrix, and u,u,ü are the displacement, velocity and acceleration vector respectively.

Material properties and model parameters
The material properties and model parameters used in the simulations are given in table 1. The material properties are valid for Bohus granite which was the rock tested in [6]. The values of the model parameters not related to Bohus granite were found by single

Spallation test simulations: intact rock samples
The spalling experiments on intact samples of Bohus granite reported in [6] are simulated here. The cylindrical rock specimen has dimensions D = 45.7 mm and L = 140 mm. The incident bar and projectile were made of aluminum and had the same diameter as the rock sample. The projectile length was L = 50 mm and a typical impact velocity was 7.3 m/s that is used here in the simulations. The finite element mesh is shown in figure 4 The simulation results along with the experimental results are shown in figure 5. The simulated rear face (free end) velocity profile agree quite well with the experimental one, as can be observed in figure 5b. The dynamic tensile strength of the rock can be calculated based on the pullback velocity, ∆v pb , as [3] σ tdyn = 1 2 ρC 0 ∆v pb (22) where ρ is the material density,  As for the failure mode, the experimental results for this particular case were, unfortunately, not available. However, the studies by Forquin and Erzar [7] and Das and Shen [8] provide suitable cases for comparison. In the former study, the authors carried out spalling tests on high strength concrete with samples of the same size as here. As the rear face velocity in a typical test was about 12 m/s, the reported double-spall failure mode compares reasonably well with the present failure mode, as can be observed in figure 5a and c. There are two main spalls in the experimental failure mode while the predicted one has an extra crack located at about 40 mm from the free end. This might be the cause for the deviation between the predicted and the experimental rear face velocity profiles. The latter study reports experimental results for Piles Creek sandstone at the impact velocity of 11 m/s. The experimental failure mode displayed the same double-spall formation with the first crack located at 76 mm from the free end and a dynamic tensile strength of 22 MPa.

Spallation test simulations: influence of cracks
The influence of two types of cracks is studied numerically in this section. First, the effect of inherent microcrack populations, modelled as spatially randomly located, oriented preembedded discontinuities with negligible strength, is tested. Then, the effect of structural macrocracks similar to those studied in [6] is investigated. The structural cracks are illustrated in figure 6. In case of inherent microcrack populations, two simulations with microcrack density of 20 % and 40 % are carried out. The orientation of all of the preembedded discontinuities is such that their normals are set parallel to the loading axis. The simulation results with the inherent microcrack populations are shown in figure 7. The microcrack populations influence the details the failure modes as well as the rear face velocity profiles, as attested in figure 7. Calculating the dynamic tensile strength with the Novikov theory (22) gives σ tdyn = 18.8 MPa (∆v pb = 3.5 m/s) in case of 20 % of elements having a microcrack and σ tdyn = 15.8 MPa (∆v pb = 2.9 m/s) when 40 % of elements have a microcrack. These strengths are 85 % and 71 % of the dynamic strength predicted for the intact specimen. It is of interest to compare these results to quasi-static tensile test simulations reported in [13]. Therein, the effect of microcrack populations was studied with exactly the same approach as here and the influence in terms of the drop in tensile strength was found to be almost linearly dependent on the microcrack density (the percentage of elements having an embedded discontinuity). Therefore, it can be concluded that the effect of microcracks is considerably less pronounced under dynamic than quasi-static loading conditions. This numerical finding is in line with the theoretical considerations by Denoual and Hild [23]. Therein, the authors found by numerical Monte Carlo simulations using a probabilistic damage model based on Weibull distribution that brittle materials become more deterministic as the loading rate increases. As there is a transition from single-to-multiple fracture modes upon increasing loading rate, this result means that under dynamic loading conditions the inherent flaws in the brittle materials do not affect the fracture processes as strongly as under quasi-static loading.
Finally, the simulations with the structural cracks illustrated in figure 6 are carried out. It should be noted that the location of the cracks was not exactly specified in [6]. Therefore, deviations between the simulation results and the experiments are expected. This is indeed the case, as can be observed in figure 8. The predicted rear face velocity shows a clear pullback velocity and some fluctuations thereafter while the experimental one has a more smooth behavior. However, the maximum velocity and the general trend are, excluding the oscillation in the predicted curve, still in quite a good agreement with the experiment. The dynamic tensile strength calculated with the equation (22) gives σ tdyn = 18.3 MPa (∆v pb = 3.4 m/s) which is 83 % of the predicted intact value above. Thus, again, the effect of (structural) cracks is quite small. It should be noted that the area ratio of the load bearing cross section at the structural crack to that of the sample is A crack /A intact = 0.54. Hence, neglecting the notch effect, the estimated strength reduction based on the load bearing area is 54 %. The predicted failure mode displays two clear macro-cracks located at the second and the fourth structural cracks (referencing to the free end). Unfortunately, the experimental failure mode was not available.

Conclusions
A previously developed continuum viscodamage-embedded discontinuity model for rock was employed in the simulations of spalling test for determining the dynamic tensile strength of brittle materials. A special emphasis was laid on the effect of microcrack populations and structural cracks on the dynamic tensile strength of the Bohus granite. The simulations for the intact samples were in a good agreement with the experiments and the present modelling approach based on embedded discontinuity finite elements better captures the failure mode than the genuine continuum models.
With the present modelling technique, the effect of microcrack populations always present in rocks can be easily modelled as pre-embedded oriented discontinuities. In the simulations, the effect of microcrack populations on the dynamic tensile strength was milder than expected based on the earlier quasi-static simulations. This feature is in agreement with the theoretical considerations that brittle materials become deterministic upon increasing loading rate. This phenomenon is an experimentally observed fact as well. The effect of the studied structural cracks was similar. Therefore, the present modelling technique has predictive capabilities in the numerical modelling of rate-sensitive quasi-brittle materials.