On the Strain Rate Sensitivity of Coarse-Grained Rock: A Mesoscopic Numerical Study

A numerical study on the strain rate sensitivity of coarse-grained rock fracture under dynamic loading is presented. For this purpose, the embedded discontinuity finite element method is employed as a numerical tool. Moreover, a mesoscopic description of grain boundary-grain interior structure of rock is given. Thereby, the present approach is able to account for interand intragranular failure types of rock. The numerical simulations carried out here corroborate the conception that in direct tension the dynamic increase of tensile strength of rock is a real material property. Moreover, the simulations agree with the hypothesis that in uniaxial compression the dynamic increase of compressive strength is a structural effect due to lateral inertia. Finally, the numerical simulations of the dynamic Brazilian disc test suggest that structural effects also contribute to the dynamic increase in the apparent indirect tensile strength.


Introduction
Rocks are strongly strain rate sensitive materials, which makes their behavior under dynamic loading challenging from the geotechnical engineering point of view (Forquin 2017;Li et al. 2017Li et al. , 2018Qian et al. 2009). When the strain rate increases, the strain rate sensitivity of rock is realized as an apparent increase in both tensile and compressive strengths accompanied with a transition from single crackto-multiple crack/fragmentation failure mode (Denoual and Hild 2000;Xia and Yao 2015;Zhang and Zhao 2014;Zhang et al. 1999). However, the micromechanisms behind the dynamic increase in strength of rock are still unclear and open to research (Forquin 2017).
In uniaxial tension at low strain rates, is rather easy to explain the specimen strength to be dictated by the largest favorably oriented flaw, which grows at the expense of others, leading to a single macrocrack failure mode as illustrated in Fig. 1a (Ahrens and Rubin 1998). Upon increasing loading rates, there is enough time for many flaws to grow and create zones of reduced stress around them, see Fig. 1a. Due to the limited crack propagation velocity (less than the speed of sound in the material), higher strain rates leads to higher tensile strengths and higher number of fragments. Therefore, the micro-to-macroscale development of the cracking process requires a finite time, called incubation time by Petrov and Morozov (1994). In unconfined compression, perhaps the most frequent explanation in the literature (see, e.g. Gary 2013) for the dynamic increase in rock strength is that the lateral inertia of the specimen acts as a confinement (through the Poisson effect) leading to the strain rate hardening (Fig. 1b). Idealistically, the specimen fails with a multiple axial splitting mode and the rock column-shaped fragments are forced to move laterally until they collapse or fly out (Fig. 1b). Therefore, higher loading rates create higher lateral inertia leading to higher peak stresses.
In this paper, the strain-rate effects of rock are numerically investigated. Two research hypotheses motivated by the discussion above are posed: (1) in uniaxial direct tension, the dynamic increase in strength is a genuine material (or microstructural) property, not a macrostructural effect. (2) In uniaxial dynamic compression, the dynamic increase in strength is a structural effect boosted by the specimen length/width ratio. This can be tested by increasing the specimen width to length ratio in the numerical tests. Saksala (2018) answered these questions affirmatively in a numerical study. These results are included here for completeness sake.
The dynamic testing methods for tensile strength include direct tension test, indirect tension tests based on the Brazilian disc test and the semicircular bending test, as well as the spalling test (Zhang and Zhao 2014). However, the experimental results are discrepant so that at a given strain rate, the Brazilian disc test gives higher strengths than the direct tension test. Due to the indirect nature of Brazilian disc test, it should be investigated whether the lateral inertia or some other structural effect contributes to the apparent increase in dynamic tensile strength in this test. Thus, the third hypothesis numerically tested in this paper is that (3) the dynamic increase in tensile strength in Brazilian disc test is (at least substantially contributed by) a structural effect.
In order to investigate these three research questions, a constitutive description of rock capable to account for strain rate dependent mesofailure types, i.e. intra-and intergranular failures, should be employed. For this end, the rate dependent mesomechanical model based on multiple embedded discontinuity finite elements presented by Saksala (2018) is chosen here as well. This model represents the mineral grains explicitly as Voronoi cells and grain boundaries by zones of elements with pre-embedded discontinuities. Moreover, the material failure, i.e. crack opening, is rate-dependent by viscosity. As mesoscopic approach by finite elements describing the grains and grain boundaries explicitly is taken, the numerical rock has to be coarse grained for computational feasibility reasons.
The paper is structured as follows. First, the numerical modelling theory is briefly described for the convenience of the reader. Then, the earlier results are summarized and the new simulations and analyses for the Brazilian disc test are presented. Finally, the concluding remarks close the paper.

Numerical Model
The rock fracture modelling principles as well as the grain structure description of the numerical rock are outlined in this section. The kinematics related to the embedded discontinuity method are presented only from the finite element implementation point of view. For a more comprehensive discussion on the displacement discontinuity theory, the early works on the method by Oliver (1996), Simo and Oliver (1994), Simo et al. (1993) are recommended. Moreover, further details can be found in earlier applications of the embedded discontinuity finite elements in rock fracture modelling, e.g. by Saksala (2015Saksala ( , 2016.

Rock Fracture Modelling
The embedded discontinuity method is an element based enriched finite element method suitable for modelling discontinuities, such as cracks. This method is particularly attractive as it can be formulated in the spirit of computational plasticity models (Mosler 2005;Radulovic et al. 2011). Therefore, the computational efficiency of continuum methods is retained while a superior fracture description is gained. The advantage of this method over the widely used discrete element/particle methods is the computational efficiency. However, the particle methods are naturally better suited for modelling severe fracturing and fragmentation.
The method describes crack initiation by embedding a displacement discontinuity in a finite element upon violation of, e.g. the Rankine criterion. Then, a cohesive law or a softening model that relates the traction to the displacement jump at the discontinuity models the opening of the real crack. In order to enhance the performance of the method in multiple crack description, an element with many intersecting embedded discontinuities was developed (Mosler 2005). The linear three-node triangular element, known also as the constant strain triangle (CST) (see Fig. 3), is selected here. For this element embedded with three discontinuities, the finite element discretized displacement and strain fields are (Mosler 2005;Radulovic et al. 2011) (1) . 1 Schematics of rock fracture mechanisms in quasi-static and dynamic uniaxial tension (a), and dynamic uniaxial compression (b) where k d is the opening vector (displacement jump) for discontinuity k, and N i and e i denote the regular shape functions and nodal displacements, respectively. Moreover, H k Γ d and k Γ d are the Heaviside function and the Dirac delta function. It has also been assumed that ∇ k d ≡ 0 when arriving at strain in (2). This is the usual elementwise constant displacement jump assumption. In addition, the peculiar decomposition using functions M k (1) is used because it simplifies the treatment of the essential boundary conditions by limiting the influence of k d to that element only. This is due to the fact that M k It is instructive to illustrate the displacement decomposition (1) and the related functions in 1D case. This is shown in Fig. 2 for a case of a single two-node bar element under tension. In the left, the functions involved in the decomposition are plotted. The decomposition consists of the regular nodal displacement, u reg = N 1 u 1 + N 2 u 2 , and the enhanced, discontinuous part, d M Γ d , the effect of which is restricted inside the element by function M Γ d . It should be noticed that in 1D case, the selection of function ϕ is readily identifiable as the interpolation function of node 2, i.e. N 2 .
The selection of ramp functions k appearing in M k Γ d in 2D case is as follows. For the multiple discontinuity case, k = N k , i.e. the ramp functions are identified with the interpolation functions of the opposite nodes (see Fig. 3). The crack normal thus reads i = ∇N i ∕N i . In the single discontinuity case with he normal, d , parallel to the first principal (stress) direction, is chosen according to.
This criterion gives ϕ with a gradient as parallel as possible to d .
Following Radulovic et al. (2011), the weak, variational form of the traction continuity over the crack is enforced with the enhanced assumed strains concept. After some developments the details of which are not presented here, the weak (global) form of the traction balance reads: where A e is the element area, and l i d denotes the discontinuity length with i as its normal. Due to the linearity of the CST element, the integrands in (4) are constant and, thus, it becomes the strong (local) traction balance, which projects the stress tensor, σ, on the crack normal: where E is the elastic stiffness tensor, i d is the traction vector, and ε is the strain tensor. The rock behavior is taken as linearly elastic before a discontinuity is introduced or the tensile strength of the pre-embedded discontinuities is reached.
As mentioned above, a model that relates the traction at the discontinuity to its opening to govern the softening process is needed. Here, a plasticity theory based model is employed for solving the crack opening vector and the traction update. This model consists of the following components (Saksala 2018): where i are the loading functions having contributions from tension and shear components of the traction, i is the tangent vector of a discontinuity, i ,̇iare the internal variable and its rate, and σ t and s are the tensile limit stress and the (constant) viscosity of the material. Moreover, h i is the softening modulus defined by the exponential softening law. Parameter g, defined by the mode I fracture energy G Ic , controls the amount of softening, and β is the shear retention factor. By (9), the viscoplastic multiplier ̇i is identical to ̇i . Finally, Eq. (10) are the consistency conditions, implying that viscosity is accommodated with the consistency approach by Wang et al. (1997).
The formal similarity of the model described in Eqs. (6)-(10) to plasticity models allows it to be solved with the computational plasticity techniques of stress return mapping. This means that the usual elastic predictor-viscoplastic corrector split is applied. Accordingly, the trial elastic stress for each element in the mesh is first calculated based on the new total strain coming from the  (6) is violated, a viscoplastic correction, i.e. the stress return mapping, is needed. During the stress return mapping, the trial stress, or actually a trial traction here, is iteratively returned on the corresponding loading surface or their intersection in case more than one criterion is violated. During the return mapping, the consistency recovering traction vector, crack opening displacements and internal variables are updated by Eqs. (7)-(9). For further details, see Saksala (2016) and Wang et al. (1997).

Mesoscopic Description of Rock
Depending on loading conditions, rocks display both interand intragranular failure types (Mahabadi et al. 2014). Therefore, these failure types should be accounted for in the numerical modeling aiming at realistic predictions. In this paper, the mesoscopic description by Saksala (2018) is employed. This approach describes the rock mesostructure as a cluster of grains represented by Voronoi cells meshed with the standard CST elements. The grain boundaries are modelled as narrow strips of elements with at least one node connected to an edge of a Voronoi cell representing a grain. Then, discontinuities are pre-embedded (inserted before the simulation) into the elements designated to belong to the grain boundaries, see Fig. 3. These pre-embedded discontinuities govern the intergranular fracture type. The intragranular failure is described by inserting new discontinuities (microcracks) during the analysis inside the grains (elements inside the Voronoi cells) according to the Rankine criterion, i.e. when the first principal stress exceeds the tensile strength, a crack is embedded inside an element. A single discontinuity governed by model (6)-(10) with i = 1 is allowed per a grain-interior element.

Explicit Scheme for Solving the Equations of Motion
Since strain rate effects of rock are considered, the governing equations of motion are solved with an explicit time integrator. The explicit modified Euler method (Hahn 1991) is selected for time integration. According to this scheme, the system response is calculated as where u is the global displacement vector, Δt denotes the time step, ext t is the global external force vector, int,e t is the t+Δt = t + Δṫt +Δt , internal force vector for a finite element e, and B e is the operator mapping the element nodal displacements into the element strains. Moreover, A is the assembly operator assembling (summing) the element contributions into a global, system level objects, and, finally, M is the lumped mass matrix obtained by the row summing technique. Now that the theory of the rock material modelling principles is presented, the flow of the total simulation process is illustrated graphically for a convenience of the reader in Fig. 4. The simulation starts with the situation where the rock mesostructure is generated and the discontinuities are embedded in the grain boundary elements (Fig. 3).
The stress return mapping referred in the flow chart is the one described in Sect. 2.1. Moreover, the criterion to be checked for the grain interior element, i.e. σ 1 > σ t is the Rankine, i.e. the first principal stress, criterion.

Numerical Simulations
The hypotheses presented in Sect. 1 are numerically tested in this Section. As mentioned, the results presented earlier on uniaxial tension and compression in Saksala (2018) are included for completeness sake and further discussed here. Then, the new simulations on Brazilian disc test are presented and analyzed. A Matlab code was written for performing the simulations. Table 1 tabulates the material and model parameters used in the simulations if not otherwise stated. These values may not correspond to any real rock but are chosen for demonstrative purposes only. The grain boundaries are assumed to be significantly less stiff and weaker than the grain interiors. Moreover, the viscosity modulus given in Table 1 is, on one hand, small enough not to cause significant strain-rate hardening at the strain-rates tested. On the other hand, it is large enough to have a stabilizing effect on the simulations.

Uniaxial Compression and Tension Tests
Numerical uniaxial compression and tension tests were carried out in Saksala (2018) at different loading rates on samples with two different specimen length/diameter ratio. The average grain size of the numerical rock samples shown in Fig. 5a is 5-6 mm with the average element side length of 0.65 mm. The grain boundary zones are shown in red therein. At low strain-rate, the loading is applied as a constant boundary velocity with an initial acceleration stage to diminish the elastic wave oscillations due to reflections at the boundaries of the numerical sample. However, at high strain-rates in tension, a linear initial velocity field of form v y (y) = 2v 0 y/h (with v 0 = ̇h/2) is imposed at each node in y-direction (see Fig. 5b). This setting prevents the stress wave oscillations at the top and bottom edges as well as premature failures close to the edges.
The simulation results for uniaxial tension test with Mesh1 at 1 s −1 and 10 s −1 are shown in Fig. 5c, d. The predicted failure mode, represented as localization bands of the sum of the norms of the displacement jump vectors (crack opening), is the experimental transversal splitting mode with a single crack at strain rate 1 s −1 and the corresponding tensile strength is 3.2 MPa. This value is well below the 5 MPa given for grain interiors in Table 1 but close to the assumed grain boundary strength of 3.5 MPa. Thus, the crack path zigzags along the grain boundaries. At strain rate 10 s −1 , the predicted tensile strength is the same 3.2 MPa but there are more macrocracks, which fact reflects a more ductile stress-strain response. However, at this strain rate the experimental DIF ranges from 1.5 to 2 (Zhang and Zhao 2014). Therefore, stronger rate dependency in the model is required to predict the correct tensile strength. Indeed, when the viscosity was increased from the value in Table 1 to s = 0.2 MPa s/m, the tensile strength is 5.2 MPa yielding DIF = 1.63, which is experimentally acceptable. However, the sample is fragmented, which is reflected in a somewhat overly ductile stress-strain response. In any case, these results demonstrate that the dynamic increase in strength of rock in direct tension is not a structural effect but a genuine material effect since a higher viscosity value was needed to predict the correct dynamic tensile strength.
Uniaxial compression tests are simulated at loading rates 1, 10, and 30 s −1 on the specimens shown in Fig. 5a. The results are shown in Fig. 6.
The simulation results in Fig. 6a, b predict again that when the strain rate increases, the amount of macrocracks increase, which is the experimentally observed transition from single-to-multiple cracking. The realized failure mode is the commonly observed shear banding with long, slightly aligned macrocracks spanning the specimen axially. Unlike the tension test simulations, here the maximum compressive strength increases quite substantially for both meshes with the low value of viscosity. It can also be observed that as the strain rate increases, the amount of cracks seems to accumulate to the vertical edges of the numerical specimen, which is also an experimental observation (Zhang and Zhao 2014).
The dynamic increase factors (DIF) for both meshes are calculated in Table 2. As the increase factors for Mesh2 with L/D-ratio 1 are larger than those for Mesh1, a clear shape effect is attested. At the highest strain rate simulated here,  (c) i.e. 30 s −1 , the experimental DIFs are scattered between 1.3 and 2.4 (Zhang and Zhao 2014). Therefore, the predicted values in Table 2 are within experimental bounds. Figure 7a, b shows an example of grain interior cracks in tensile and compressive loadings, respectively. In tension, the orientation trend of the cracks is horizontal, and in compression it is vertical (the major principal direction in uniaxial compression) with some deviations caused by disturbances in the stress field due to cracking or grain boundary heterogeneities.
Finally, the transversal (x) and axial (y) stress distributions in Mesh2 at the peak stress are shown in Fig. 7c. The stress component distributions are severely heterogeneous due to the grain-boundary heterogeneity and the cracking. Moreover, compressive stresses can be observed at the grain boundaries in the distribution of σ x , especially in the center area of the specimen. As the magnitude of these stresses increased upon increasing loading rate, their presence is the reason for the damage (in the sense of stronger crack opening) accumulation on the specimen vertical edges in Fig. 6. As this effect was stronger in higher strain rates and with wider specimen (smaller L/D-ratio), a specimen shape effect is clearly present.

Brazilian Disc Test Simulations
The hypothesis of structural effects (lateral inertia) contribution in the indirect tensile strength dynamic Brazilian disc (BD) test is numerically tested here. The schematic of the computational model is shown in Fig. 8.  . 7 Example of grain interior cracks in the end of tension (a), and compression (b) test simulation, and the stress components at the peak stress (~ 93 MPa) in compression test simulation at 30 s −1 with Mesh2 (c) In this modelling approach, the compressive stress pulse is accounted for as a pre-defined stress function σ i (t). The bars (with cross-sectional area A b in Fig. 8) of the split Hopkinson pressure bar device are meshed with two-node bar elements and the mesomechanical description of rock is applied to the Brazilian disc sample of rock. Finally, kinematic contact conditions are imposed between the bars and the disc. During the analysis of the impact problem, the contact forces P 1 and P 2 are solved using the method described in Saksala (2010).
The tensile stress at the center of the disc is where F is the force compressing the disc with dimensions t = 25 mm (thickness) and d = 40 mm (diameter). In the present setting where the sample is placed between the bars, the simple average of P 1 and P 2 is used for F. Both bars in the model have length and diameter of 1500 mm and 22 mm, respectively. The stress pulse simulating the compressive wave is defined as σ i (t) = A p sin(ωt) with A p being the amplitude and ω = 2π/T. In the simulations, the pulse length T = 160 µs is used throughout, but the amplitude is varied to obtain different strain rates.
The model data given in Table 1 are used in the simulations. First, a reference case simulation is carried out at strain rate 0.1 s −1 . The boundary conditions similar to those used in the uniaxial compression test are applied here as well. The simulation results are shown in Fig. 9.
The low rate BD test simulation predicts the correct axial splitting failure (see Fig. 9c). The major splitting crack has a slightly inclined trend spanning the specimen from the left side of the upper contact area to right side of the lower support area. This is due to the polygonal (22-gon) geometry of the numerical disc. This geometry results in flattening, which should actually be taken into account in calculating the indirect tensile strength with Eq. (14). However, the (14) T = 2F∕ td, flattening here is equivalent to the loading angle 2α = 16°, which yields a correction factor of 0.9744 (Wang et al. 2004).
The indirect tensile strength corresponding to the maximum force in Fig. 9d is σ T = 3.86 MPa (with the flattening correction), which is 20% higher than the predicted direct tensile strength. It should be taken into account that the mesostructural details are not identical in the numerical samples. However, the results are comparable since the grain and element sizes are of the same order (the average grain size is 5-6 mm end the average element side length is about 0.65 mm) in all of the simulations. Moreover, the experimental indirect tensile strengths are systematically, and in an increasing manner upon increasing rate, larger than the direct ones starting from strain rate 0.1 s −1 (Zhang and Zhao 2014).
Next, higher strain rates are applied using the model described in Fig. 8 with the numerical sample in Fig. 9a. Amplitudes A p = 25 MPa and 50 MPa of the compressive stress wave are used. The results are shown in Fig. 10.
At higher loading rates, the details of the predicted general axial splitting failure mode change from a single major crack mode to the multiple cracking mode, as attested in Fig. 10a, b. At the highest rate, corresponding to amplitude 50 MPa, even some radial cracks are generated. The stress rates in these simulations were 76 GPa/s and 143 GPa/s corresponding to the strain rates 1.5 s −1 and 2.8 s −1 , respectively. The dynamic equilibrium, i.e. P 1 = P 2 , was fairly well achieved in these simulations. Therefore, Eq. (14) can be used to calculate the indirect dynamic tensile strengths, which were 4.6 MPa and 6.5 MPa at the lower and higher stress rate, respectively. These give, when referenced to the low strain rate simulation above, DIFs of 1.2 and 1.7 respectively. These values are slightly lower than the typical experimental DIFs, which are ~ 1.5 and ~ 2.0 at these strain (stress) rates (Zhang and Zhao 2014). In any case, the predicted DIFs are definitely not negligible and indicate a

Concluding Remarks
A numerical study on the loading rate sensitivity of coarse-grained rock under direct and indirect tension as well as under uniaxial compression was carried out. The chosen mesomechanical approach can account for both the inter-and intragranular failure types of rock. The numerical simulations demonstrated in general that this method predicts some of the important features, such as failure modes and DIFs, of rock behavior under dynamic loading. Therefore, as the present model has predictive capabilities, some credibility could be bestowed on the numerical confirmations of the particular hypotheses tested here.
Namely, it seems that the dynamic increase in the apparent strength of rock under direct tension is a real material (or at least microstructural) effect. In contrast, the dynamic increase of rock under uniaxial compression seems to be a structural property. Finally, structural effects probably contribute to the dynamic increase in indirect tensile strength as it is measured in the Brazilian disc test. However, the structural effects are weaker in the BD test than in the uniaxial compression test so that the material effects contribute to the indirect tensile strength. This presence of the structural effect could be proposed as an explanation to the discrepancy in the experimental results, which consistently attest higher DIFs for the rock under indirect than direct dynamic tension. Strictly speaking, these results are valid for the coarse-grained numerical Fig. 9 BD test at strain rate 0.1 s −1 : the mesh with 10,662 elements (a), grain boundary/interior structure (b), final failure mode (c), and the force-time curve (d) rock specimen. Therefore, more investigations are needed to bring them to bear on other rock types.