1 Introduction

An understanding of the dynamic fracture process of rock materials is of importance for efficient rock drilling in several applications e.g. in mining and geothermal industries. It is estimated that 50% of the total cost per installed megawatt geothermal energy is associated with drilling and well construction [3], with the predominant cost factor being associated with drill bit wear [21]. Knowledge and understanding about the full-scale rock drilling processes can be gained from multiphysics simulations, where rock fracture, wear and rock cutting transport by fluid is taken into consideration. These simulations can therefore be used to optimise drilling operations in terms of rate of penetration and service life. Since rock materials are heterogeneous e.g. due to varying mechanical and geometrical properties of the constituent minerals, an accurate rock fracture model needs to be able to account for this.

The initiation, propagation and coalescence of cracks are well known to be governing factors for rock excavation processes [25, 27]. A common laboratory rock experiment is the Brazilian disc test (BDT), which is used to determine the indirect tensile strength of rock materials [12]. Furthermore, the BDT constitutes a suitable test for evaluating the initiation, propagation and coalescence of cracks in rock materials. For rock drilling, which is a dynamic process, the BDT should be conducted with dynamic loading. The split-Hopkinson pressure bar (SHPB), or Kolsky Bar [8], is an experimental configuration commonly used to evaluate the mechanical behaviour of rocks at high strain rates (10\(^{1}\)-10\(^{3}\) s\(^{-1}\)) [35, 40]. For dynamic tensile testing of rock materials, the BDT is often used due to the difficulties with direct uniaxial methods [12], and is the suggested method by the International Society for Rock Mechanics and Rock Engineering (ISRM) for obtaining the dynamic tensile strength [40]. Today, the BDT is widely used to characterise the dynamic tensile behaviour of various rock materials. Xu et al. (2021) studied the morphological and temporal aspects of fracture in igneous rocks and obtained pixel traces of the fracture path in the sample [38]. Li et al. (2021) demonstrated a novel method, using high-speed photography and digital image correlation (DIC), for obtaining the dynamic tensile stress–strain curve of granite [14]. Although being the most practical method for tensile testing of rock, the validity of the BDT has long been debated [12, 20], mainly due to improper points of crack initiation [4]. Furthermore, Xia et al. (2017) used strain gauges in a dynamic BDT of granite to identify the time of crack initiation and demonstrated that fracture initiated at 20 % of the peak stress in the sample [37]. Since the Griffith condition is only valid up until the point of fracture initiation, the excess peak stress does not represent the tensile strength in the rock material and should be accounted for. This method has also been used by Li et al (2020) where it was demonstrated how the overloading effect increases with loading rate [13]. This overloading effect has also been observed in granite and diorite by Wessling and Kajberg (2022) [35].

For numerical modelling of rock materials, various works incorporating rock heterogeneity may be found. Using the Finite Element Method (FEM), Saadati et. al. [26] combined a plasticity model of compression and a damage model for tension, where rock defects were accounted for using the Weibull distribution. This approach has also been used to study the effects of preexisting cracks in granite [19]. Also using the FEM, Saksala et. al. [29] represented different minerals in a granite rock by randomly distributing polygonal elements. Using a damage-viscoplastic model, the authors simulated a dynamic BDT and managed to capture the correct axial splitting failure. A similar approach has been used to simulate dynamic BDT in a SHPB configuration [28].

In contrast to the aforementioned continuum based methods, the bonded particle model (BPM) [23] represents a continuum through discrete, rigid spheres that are bonded together with linear elastic beams. These beams can transmit load in tension, shearing, twisting, and bending, and material deformation and fracture depend on microscopic stiffness and strength parameters. Fracture initiates at a micromechanical scale by one broken bond which then propagates through more broken bonds and coalesce with other cracks to form larger fracture surfaces. The BPM has been used to evaluate common experimental methods for rock, such as the BDT [17], uniaxial compression [11], as well as industrial processes such as rock cutting [24, 31]. Unfortunately, the BPM has been found to suffer from not being able to accurately represent rock materials with a high compressive to tensile strength ratio [22], one of the reasons being due to the spherical shape of the discrete elements used to represent the rock material [23]. However, various solutions to this problems have been identified, such as increasing the interlocking range [30], representing more complex grains through clumped particles [1, 18] or by introducing smooth-joint contacts between some particles that allow for overlapping and sliding [6].

Recently, Wessling et. al. [36] developed a heterogeneous clustering model for the BPM which takes both geometrical and mechanical heterogeneity into account. The rock grains in this model are represented as random ellipsoidal subsets of bonded particles, with mechanical properties being governed by the Weibull distribution. It should be noted that this is an holistic approach for representing a heterogeneous rock material. This means that specific mineral types are not accounted for, which has been done by other authors e.g. by using polygonal grain-based discrete approaches in 2D [10, 39].

The main objective of this work is to use the model developed in [36] to study the effects that sample heterogeneity and grain cement strengths has on dynamic rock fracture. The modelling approach is also further developed by utilising three Weibull distributions for assigning mechanical properties to the grains and boundaries. Dynamic rock fracture accounting for heterogeneity has been simulated before (see e.g. [15, 29]); however, no studies can be found regarding the repeatability of measured quantities and the initiation and propagation of cracks. To this end, the repeatability is investigated in this work. It is important to note that a constant loading rate (1400–1500 GPa \(s^-1\)) is investigated in this study. That is, strain rate dependency is not considered in this work. First, a homogeneous BPM is calibrated towards dynamic experimental data on Kuru black diorite, where the macroscopic stiffness and strength as well as fracture mode and indirect tensile strain are matched. Secondly, the grain structure of the rock is represented using ellipsoidal subsets with an average size corresponding to that of the real rock material. A parametric study is conducted on the Weibull heterogeneity index of the average grain properties and the cement strength, and the results are compared to experiments in terms of variations in predicted tensile strength, fracture path in the sample and time of crack initiation. In addition, the magnitudes and repeatability of overloading, strain rate and crack propagation time are evaluated. By ensuring that this modelling approach is able to replicate the fracture process in a dynamic rock experiment, it can be used to study larger industrial rock excavation processes. For example, the model could be applied to a rock drilling process and used to evaluate performance and service life.

2 Materials and methods

In this section, the methodology of the work is presented. This includes a description of the studied rock material and the numerical set-up of the BDT in a SHPB configuration. Furthermore, the fundamental theory of the bonded particle model (BPM) and the heterogeneous BPM as well as the calibration of these are explained.

Fig. 1
figure 1

The surface of Kuru black diorite investigated in this study [35]

Full size image
Table 1 Dynamic mechanical properties of Kuru black diorite [35]
Full size table

2.1 Rock material

The rock material studied in this work is Kuru black diorite, which is a massive rock with an average grain size of 2.4 mm. The surface of a Kuru black diorite rock is shown in Fig. 1. This rock has previously been subject to a thorough experimental characterisation with dynamic loading in a SHPB configuration [35], whose results are summarised in Table 1. The authors combined the SHPB, high-speed photography (HSP) and digital image correlation (DIC) in order to obtain the tensile strength, tensile strain and instant of crack initiation in the BDT. Note that the authors differentiate between the peak stress and tensile strength by taking the overloading effect into account.

2.2 Split-Hopkinson and the Brazilian disc test

Fig. 2
figure 2

An experimental schematic of a split-Hopkinson pressure bar

Full size image

The simulations of the dynamic Brazilian disc test are conducted by replicating the SHPB configuration used in the dynamic experimental characterisation of diorite [35]. A brief explanation of the simulated SHPB follows below. The SHPB configuration consists of a 3 m long incident bar and a 2 m long transmitted bar of maraging steel with diameters 25.4 mm, between which the sample is placed (see Fig. 2). In the experiments, a projectile is accelerated by an air gun and impacts the incident bar, which generates a compressive stress wave. This wave propagates through the incident bar and is partly transmitted through the sample and partly reflected back into the incident bar. From strain gauge measurements and by assuming one dimensional wave propagation, the forces exerted on the sample can be obtained. The full length of the bars are represented in the simulations using linear elastic FEM with 180,000 and 120,000 hexagonal solid elements for the incident and transmitted bar, respectively. By including the full length of the bars, the incident, reflected and transmitted force pulses can be obtained and used to evaluate the force equilibrium in the sample using the same method as in the experiments. Force equilibrium is necessary to fulfil according to the fundamental theory of the SHPB experiment [5]. Instead of simulating the projectile striking the incident bar, a force pulse is generated on the bar surface. This force pulse corresponds to an average of five experimental incident pulses (see Fig. 3).

The Brazilian disc test is an experimental method for obtaining the indirect tensile strength of rock materials. Here, an indirect tensile stress is induced at the centre of a disc by compressing it diametrically. The tensile stress at the centre part of the disc can then be obtained as [40]

$$\begin{aligned} \sigma _{t}=\frac{2P(t)}{\pi Dt} \end{aligned}$$
(1)

where P(t) is the load (see below). The diameter D of the simulated samples is 25 mm and the thickness t is 12.5 mm, both of which are in agreement with the suggested standards for dynamic tensile testing [40]. It should be noted that Eq. 1 is only an approximation of the tensile strength as obtained from direct methods. However, due to being more practical than direct methods, the Brazilian disc test is most commonly applied for tensile testing of rock materials. In order for the test to be accurate, the main splitting crack should initiate close to the centre part of the disc. However, due to the heterogeneity present in rock materials, the crack initiation point has been observed to deviate from the centre [4], even initiating from the loading plates [33]. Furthermore, when dynamic loading is used, the sample has been found to continue to carry significant load after crack initiation, after which the Griffith failure criteria is invalid, which results in an overestimation of the dynamic tensile strength [13, 35, 37].

Fig. 3
figure 3

Five experimental input pulses and the averaged input pulse used for the simulations

Full size image

By utilising an incident pulse in the simulations similar to the experiments (see Fig. 3), and evaluating the reflected and transmitted force pulses following the same procedure as in the experiments, direct comparisons can be performed. The force at the interface between sample and incident bar is given by

$$\begin{aligned} P_1=p_i(t)+p_r(t) \end{aligned}$$
(2)

where \(p_i(t)\) and \(p_r(t)\) are the incident and reflected force pulses evaluated in the simulations. The force at the interface between sample and transmitted bar is given by

$$\begin{aligned} P_2=p_t(t) \end{aligned}$$
(3)

where \(p_t(t)\) is the transmitted force pulse evaluated from the simulations. It is important to note that force equilibrium in the sample is crucial for achieving reliable measurements. The force equilibrium can be evaluated by comparing the incident force \(P_1\) with the transmitted force \(P_2\). An example of this equilibrium evaluation follows in Sect. 3.

The load P(t) in Eq. 1 is evaluated as the average of \(P_1\) and \(P_2\) (3-wave analysis [5]), i.e.

$$\begin{aligned} P(t)=\frac{P_1+P_2}{2}=\frac{p_i(t)+p_r(t)+p_t(t)}{2} \end{aligned}$$
(4)

The indirect tensile strain is measured using a methodology proposed by Li et. al. [14]. Here, the indirect strain is measured as an average over a rectangular region at the centre of the disc. This region is determined to be large enough to avoid local effects, such as microcracks, but small enough to capture the indirect tensile strength close to the centre of the disc. The size of this region in the simulations is the same as has been determined for Kuru black diorite previously, i.e. a width and height of 3.44 and 5.19 mm, respectively [35].

2.3 Bonded particle model

The discrete element method (DEM) represents a granular material is through a set of particles that interact with each other through contact laws [2]. The motion of each particle is obtained from explicit integration of Newton’s second law,

$$\begin{aligned} m_{i} \frac{d \varvec{v}_{i}}{d t}=\varvec{F}_i \end{aligned}$$
(5)

where \(\varvec{v}_{i}\) and \(m_{i}\) are the velocity and mass of a individual particle \(i \), respectively. The total force \(\varvec{F}_i\) acting on this particle is given by

$$\begin{aligned} \varvec{F}_i=\varvec{F}_i^{ext}+\sum ^{n_j}_{j=1}\varvec{F}_{ij} \end{aligned}$$
(6)

where \(\varvec{F}_i^{ext}\) contains the external forces acting on particle \(i \), such as gravity, and \(\varvec{F}_{ij}\) contains the contact forces with the neighbouring particles \(j \). The angular motion is governed by

$$\begin{aligned} \varvec{I}_{i} \frac{d \varvec{\omega }_{i}}{d t}=\sum _{j} \varvec{M}_{i j} \end{aligned}$$
(7)

where \(\varvec{I}_{i}\) is the moment of inertia, \(\varvec{\omega }_{i}\) is the angular velocity and \(\varvec{M}_{i j}\) is the torque applied from particle j.

In this study, the commercial simulation software LS-Dyna [16] is used, where the DEM contact forces are represented by a linear spring-dashpot model [7]. There are five parameters controlling the DEM behaviour; normal and shear spring constants, \(k_n\) and \(k_s\), normal and tangential damping, \(c_n\) and \(c_t\), and particle-particle friction \(\mu _{p-p}\). The contact between DEM and FEM, e.g. the contact between rock and SHPB bars in this study, is a penalty based contact which is controlled through particle-FEM friction \(\mu _{p-F}\) and particle-FEM damping \(c_F\). The integrations of Eqs. (5) and  (7) are conducted using explicit integration.

Fig. 4
figure 4

A homogeneous Brazilian disc filled with 840,000 discrete element particles (a), a heterogeneous Brazilian disc filled with ellipsoidal subsets of bonded particles (b) and examples of various types of ellipsoidal subsets that can be obtained (c)

Full size image

The BPM is an extension to DEM used to model a brittle material, such as rock [23]. The rock is represented through an assembly of cemented discrete particles. The particles are cemented together with linear elastic beams, or bonds, that can transmit load in tension, shearing, twisting and bending. Furthermore, these bonds break once the stress state exceeds a predescribed value. These bonds also work in parallel to the conventional DEM. Therefore, the DEM contact is recovered once the bond breaks. A bond will break if the bond tensile or shear stress, \(\overline{\sigma }^{\max }\) or \(\overline{\tau }^{\max }\), exceeds the prescribed tensile and shear bond strengths, \(\overline{\sigma }_{c }\) and \(\overline{\tau }_{c }\), i.e.

$$\begin{aligned}&\overline{\sigma }^{\max }=\frac{-\overline{F}^{\textrm{n}}}{A}+\frac{\left| \overline{M}^{\textrm{s}}\right| \overline{R}}{I}\ge \overline{\sigma }_{c } \end{aligned}$$
(8)
$$\begin{aligned}&\overline{\tau }^{\max }=\frac{\left| \overline{F}^{\textrm{s}} \right| }{A}+\frac{\left| \overline{M}^{\textrm{n}}\right| \overline{R}}{J} \ge \overline{\tau }_{c } \end{aligned}$$
(9)

where A, I and J are the bond area, moment and polar moment of inertia of each bond. The bond radius \(\bar{R}\) is evaluated as the minimum radius of the two bonding particles, i.e. \(\overline{R}=\overline{\lambda }\,\min (r^1,r^2)\), with the bond radius multiplier \(\overline{\lambda }\) set to unity in this study. For simplicity, a homogeneous particle radius \(r=0.1\) mm is used in this study. Therefore, the bond radius is constant throughout the rock body. The bond forces, \(\overline{F}^n\) and \(\overline{F}^s\), and moments, \(\overline{M}^n\) and \(\overline{M}^s\), are incrementally updated for each time step according to

$$\begin{aligned} {\left\{ \begin{array}{ll} \Delta \overline{F}^{n}=\overline{k}^{n} A \Delta U^{n} &{} \\ \Delta \overline{F}^{s}=-\overline{k}^{s} A \Delta U^{s} &{} \\ \Delta \overline{M}^{\textrm{n}}=-\overline{k}^{s} J \Delta \theta ^{n} &{} \\ \Delta \overline{M}^{s}=-\overline{k}^{n} I \Delta \theta ^{s} &{} \end{array}\right. } \end{aligned}$$
(10)

where n and s represent the normal and shear direction, respectively. In Eq. (10), \(\overline{k}^{n}\) and \(\overline{k}^{s}\) correspond to the normal and shear bond stiffness, respectively, and \(\Delta U^{n/s}\) and \(\Delta \theta ^{n/s}\) correspond to the relative displacement and rotation. Furthermore, damping is present in the bonds through a local non-viscous damping force \(\varvec{F}^d_i\) given by

$$\begin{aligned} \varvec{F}^d_i=-\overline{\alpha } |\varvec{F}_i| {{\,\textrm{sign}\,}}(\varvec{v}_{i}) \end{aligned}$$
(11)

where \(\overline{\alpha }\) is the bond damping coefficient. Finally, the amount of particles being bonded together initially is governed through the interlocking range \(\overline{\beta }\), i.e. two spherical particles are bonded together if the distance between their perimeters is less than \(\overline{\beta }\).

To summarise, there are 14 parameters governing the DEM and BPM. Apart from the particle radius, these parameters consists of the 7 parameters associated with the DEM \(k_n\), \(k_s\),\(c_n\), \(c_s\), \(\mu _{p-p}\), \(\mu _{p-F}\) and \(c_{F}\), as well as the 6 parameters associated with the BPM, \(\overline{k}_n\), \(\overline{k}_s\), \(\overline{\sigma }_{c}\), \(\overline{\tau }_{c}\), \(\overline{\beta }\) and \(\overline{\alpha }\). The choice of each of these parameters is presented and motivated in Sect. 2.5.

2.4 Heterogeneous bonded particle model

Fig. 5
figure 5

A schematic view of six example grains and boundaries and their associated Weibull distributions

Full size image

The heterogeneity is introduced in two steps in this work - a geometrical heterogeneity due to randomised grain sizes and shapes as well as a mechanical heterogeneity due to statistically distributed bond strengths and stiffnesses. The grains are represented by randomised, irregular ellipsoidal subsets of bonded particles. The grain generation process starts with a body discretised with discrete particles of uniform size, such as the one in Fig. 4a consisting of 840,000 particles. Then, the following steps are performed:

  1. 1.

    A parent discrete element particle is randomly selected and excluded from future selection.

  2. 2.

    An ellipsoid, with its centre coinciding with the parent particle, is generated using the uniform distribution of the radii, i.e.

    $$\begin{aligned} R_x,R_y,R_z\sim U(R_{min},R_{max}) \end{aligned}$$
    (12)

    where \({R_x}\), \({R_y}\) and \({R_z}\) corresponds to the ellipsoidal radius in the x, y and z direction, respectively, and \(U(R_{min},R_{max})\) is the uniform distribution between \(R_{min}\) and \(R_{max}\).

  3. 3.

    A set is formed from all discrete element particles that lie inside this ellipsoid. This set corresponds to a rock grain. An important note is that the particles inside this set is excluded from future selection.

  4. 4.

    Steps 1 - 3 are repeated until there are no parent particles left for selection.

For the BDT sample used in this study, a uniform discrete element radius of \(r=0.1\) mm is used. Furthermore, by using an ellipsoidal subset radii of \((R_{min},R_{max})_{x,y,z}=(0.1, 3.085)\) mm, an average grain size of 2.4 mm is obtained, which corresponds to the average grain size of Kuru black diorite. An example of a generated sample can be seen in Fig. 4b. It should be noted here that this is an holistic approach to represent the grain structure of a rock material. Hence, the grain size distribution is not represented, only the average grain size. Since the fundamental shape of the discrete element subsets are ellipsoidal, some grains generated early in the process can be expected to exhibit that shape. However, as the generation process progresses, and more discrete elements are excluded from selection, more complex grain shapes are obtained (see examples in Fig. 4c).

The mechanical heterogeneity is achieved by statistically distributing the bond stiffnesses and strengths using three Weibull distributions [34]. The first distribution is used to distribute average bond strengths and stiffnesses to the grains through

$$\begin{aligned} f(\eta ,\eta _0,m)=\frac{{m}}{\eta _0}\big (\frac{\eta }{\eta _0}\big )^{{m}-1}e^{-(\eta /\eta _0)^{{m}}} \end{aligned}$$
(13)

Here, \(\eta \) corresponds to the parameters being distributed, \(\eta _0\) is the scale parameter and m is the shape parameter (or heterogeneity index). For a larger value of the heterogeneity index, the smaller the variation of the distributed parameter, and vice versa. The second distribution is used to distribute bond strengths and stiffnesses within each grain. For example, for grain i with bond parameters \(\eta _0^{g_i}\) given by Eq. 13,

$$\begin{aligned} f(\eta ^{g_i},\eta _0^{g_i},m^{g})=\frac{{m^{g}}}{\eta _0^{g_i}}\big (\frac{\eta ^{g_i}}{\eta _0^{g_i}}\big )^{{m^{g}}-1} e^{-(\eta ^{g_i}/\eta _0^{g_i})^{{m^{g}}}} \end{aligned}$$
(14)

where \(\eta ^{g,i}\) is the distributed parameters and \(m^{g}\) is the heterogeneity index within the grain. In order to represent grain cementing that is weaker than the cement itself, the bond strengths \(\overline{\sigma }^c_{c}\) and \(\overline{\tau }^c_{c}\) of the cement are obtained by scaling the prescribed strengths of the grains, i.e.

$$\begin{aligned} {\left\{ \begin{array}{ll} \overline{\sigma }^c_{c}=C_f\cdot \overline{\sigma }_{c} &{} \\ \overline{\tau }^c_{c}=C_f\cdot \overline{\tau }_{c} &{} \end{array}\right. } \end{aligned}$$
(15)

where \(C_f\) is the cement scaling factor. The third distribution is used for the cement strengths, \(\overline{\sigma }^c_{c}\) and \(\overline{\tau }^c_{c}\), between grains, given by

$$\begin{aligned} f(\eta ^c,\eta _0^c,m^c)=\frac{{m^c}}{\eta _0^c}\big (\frac{\eta ^c}{\eta _0^c}\big )^{{m^c}-1} e^{-(\eta ^c/\eta _0^c)^{{m^c}}} \end{aligned}$$
(16)

where \(\eta ^c\) is the distributed parameters, \(\eta _0^c\) is the scale parameter and \(m^c\) is the heterogeneity index for the cement. The schematics of the intragranular and cement distributions are shown in Fig. 5. A parametric study is conducted on the cement scale factor and the global heterogeneity index m. However, in order to limit the scope of this study, this is not done for the grain and cement heterogeneity indices \(m^{g_i}\) and \(m^c\). In order to obtain a rather large intragranular and cement heterogeneity [36], both of these are set to a value of five for all of the simulations.

2.5 Calibration

This section describes the calibration of the 14 parameters of the DEM and BPM, which can be an extensive process. A common approach is to use the BDT and UCT to capture the full elastic behaviour, tensile and compressive strength. In this work, however, a simplified approach, using only the BDT, is used [36]. Analytical expressions and literature values are used in order to decrease the amount of trial-and-error. The normal stiffness of the DEM has been suggested to be evaluated as [23]

$$\begin{aligned} k_n=4RE_c \end{aligned}$$
(17)

where \(E_c\) is the particle stiffness, set to 100 GPa in this study. For the shear stiffness of the DEM, elastic theory for the bulk and shear modulus is used [9],

$$\begin{aligned} k_s=\frac{k_n}{2(1+\nu )} \end{aligned}$$
(18)

By increasing the interlocking range, the amount of brittleness, i.e. the ratio between compresssive and tensile strength, can be increased. Scholtés et. al. [30] pointed out that the interlocking range should be chosen with care as to minimise the number of particles embedded between two bonded particles. With the uniform particle radius used in this study, a interlocking range of \(\overline{\beta }=0.2\) mm is used. This is the maximum distance that assures that no particle is embedded between two bonding particles. The BPM damping coefficient \(\overline{\alpha }\) is set to 0.7, which previously has been used to simulate Lac du Bonnet granite in quasi-static conditions [23]. All relevant parameters used in this study are summarised in Table 2. It is important to note that a constant loading rate (1400–1500 GPa \(s^-1\)) is investigated in this study. That is, strain rate dependency is not considered in this work.

A homogeneous version of the BDT, i.e. using the sample in Fig. 4a, is used to determine the bond stiffnesses and strengths. The normal bond stiffness \(\overline{k}_n\) and the ratio between shear and normal bond stiffness \({\overline{k_s}}/{\overline{k_n}}\) are determined so that the stiffness response and indirect tensile strain are similar to that of the dynamic experiments. It is of importance to set the mechanical properties so that the positions of crack initiations resemble the locations observed in the experiments. The shear strength is therefore chosen large enough to avoid severe shear damage or crack initiations at the contacts between sample and bars. Furthermore, the tensile bond strengths are determined as to yield the correct peak stress. A trial-and-error approach is used to achieve a response that corresponds to the experimental results. However, this process is accelerated by utilising linear dependency when applicable. For example, a linear dependency was found between the stiffness response and the normal bond stiffness, the indirect tensile strain rate and the ratio between shear and normal bond stiffness, as well as between the peak tensile stress and the normal bond strength.

A parametric study is conducted on the heterogeneity index for the average grain properties and cement scale factor (see Table 3). As heterogeneity and weaker grain boundaries are introduced into the sample, the measured indirect tensile strength decreases. For a higher heterogeneity and lower cement scale factor, the strength decreases even more. In order to match the simulated tensile strength for each case with the experimental tensile strength, the bond strengths were scaled individually using a scaling coefficient \(k_f\), which are listed in Table 3. For each of the cases in Table 3, five samples were generated and simulated, which results in a total of 45 samples. The samples were generated using Matlab (v.2022b) and the simulations were conducted using the commercial software LS-Dyna. Each simulation was conducted on a Xeon Gold 6248R 3GHz CPU with 48 cores.

Table 2 Input parameters for the simulation
Full size table
Table 3 The cement scaling factors of Eq. (15), shape parameters of Eq. (13) and strength scaling factors used for the different cases
Full size table

3 Results and discussion

In this section, the results are presented and discussed. First, results from the calibrated homogeneous BDT are presented, where force equilibrium and a good correspondence with the experiments are shown. Secondly, the results from the heterogeneous samples are shown. These results consist of the initiation and propagation of cracks in the BDT samples as well as the indirect tensile strength, strain and strain rate, and associated overloading. At the end of this section, the simulated crack paths are compared to high-speed images of an experimental test. Each of these samples was ensured to be in force equilibrium up until the point of fracture initiation. Each of these cases are compared to experimental data from five samples of Kuru black diorite. Each simulation had a total computational time of roughly 4 h.

3.1 Homogeneous Brazilian disc

Fig. 6
figure 6

The tensile stress at the centre of the disc and strain versus time for a homogeneous numerical dynamic Brazilian disc test

Full size image
Fig. 7
figure 7

The tensile stress at the centre of the disc and strain versus time for an experimental dynamic Brazilian disc test. Adapted from Wessling et al. [35]

Full size image
Fig. 8
figure 8

Force equilibrium evaluation from the split-Hopkinson simulation of a homogeneous dynamic Brazilian disc test

Full size image
Fig. 9
figure 9

Force equilibrium evaluation from the split-Hopkinson experiment of a dynamic Brazilian disc test. Adapted from Wessling et al. [35]

Full size image

In Fig. 6, the simulated stress and strain from the homogeneous BDT are shown, which is to be compared to the results from the experiments in Fig. 7 [35]. The homogeneous BDT is able to capture the tensile stress at the centre of the disc and the indirect tensile strain rate is similar to the experiment. The simulation also captures the two distinct strain rate regions seen in the experiments, the knee of which coincides well with the point of crack initiation. Furthermore, as discussed in Sect. 2.2, equilibrium in the sample needs to be ensured. Equilibrium in the sample can be evaluated by comparing the incident force \(P_1\) with the transmitted force \(P_2\), or equivalently the sum of the incident and reflected pulse with the transmitted pulse (see Eq. 2-3). This is shown for the homogeneous simulated BDT in Fig. 8, which can be compared to an example of experimental equilibrium in Fig. 9. From this it is clear that the sample is in equilibrium up until the point of fracture. Further, it is shown that the calibrated homogeneous bonded particle model is able to replicate the experiments well.

3.2 Heterogeneous Brazilian disc

This section covers the results from the 45 heterogeneous samples from all cases in Table 3. The results consist of points of crack initiation and propagation in the sample as well as the measured indirect tensile stress and strain, magnitude of overloading, crack propagation time and strain rate.

3.2.1 Crack initiation

Fig. 10
figure 10

Points of crack initiations for different cement scale factors \(C_f\) and heterogeneity indices m. The circular markers correspond to valid points of initiation while square markers correspond to invalid points of initiation

Full size image

In this study, only the initiation of the main axial splitting crack was documented. The crack initiation for each case is shown in Fig. 10. The majority of the crack initiations occurred at one of the flat sample surfaces. In most of the simulated samples, several cracks initiated at the same time, which is consistent with dynamic fracture of rock materials [32].

Out of all 45 samples, three samples exhibited main crack initiation from damage originating from the contact with the Hopkinson bars. These are denoted by square markers in Fig. 10. Therefore, all but three tests can be considered valid. Although, the accuracy of the indirect tensile strain measurement of the tests with initiation points further away from the centre might be debated. In contrast to when this modelling approach is subjected to quasi-static loading [36], no apparent pattern regarding initiation points and statistical parameters can be discerned from the initiation points. However, there are more initiation points on the incident side than on the transmitted side (30 vs 15 samples) in the SHPB configuration. This suggests a slight force unbalance in the sample. However, for all practical purposes, the equilibrium achieved in the simulations is sufficient (see Fig. 8).

3.2.2 Crack propagation

Fig. 11
figure 11

A representative example of the crack initiation and propagation in a typical heterogeneous sample, with heterogeneity index \(m=15\) and cement scale factor \(C_f=0.50\). Crack initiation (a), before the peak stress (b), at the peak stress (c) and the corresponding stress-time curve (d). It is important note is that only discrete elements that has lost at least 85 \(\%\) of their bonds with neighbouring particles are shown in (a)–(c). Furthermore, the cylinder cross-section is perfectly circular for simplicity in (a)–(c), while in reality it should be slightly compressed

Full size image

An example of the propagation of the main splitting crack can be seen in Fig. 11, with the associated tensile stress versus time curve. In Fig. 11a–c, only discrete elements with a damage parameter of at least 85 \(\%\) are shown. In other words, a discrete element is only shown if it has lost a minimum of 85 \(\%\) of its initial bonds with neighbouring particles. Here, it can be seen that one main crack initiates close to the centre of the disc, which then propagates horizontally towards the SHPB contact interfaces. At the same time as the main crack, two smaller cracks initiates as well. These smaller cracks does not propagate through the whole sample, rather, they stop as soon as the main crack is fully propagated. By comparing the cracks in Fig. 11a–c with the tensile stress curve in Fig. 11d, it is clear that crack initiation occurs well before the peak stress. This is known as the overloading effect that dynamic BDT exhibits [13, 35, 37]. This means that the peak stress in a dynamic BDT does not correspond to the actual tensile strength of the rock material. Instead, the peak stress corresponds to the point in time when the main splitting crack has fully propagated through the sample. From this, the crack propagation time, i.e. the time between initiation and a fully propagated splitting crack, in the sample can be estimated as the time between initiation and peak stress. For the particular example in Fig. 11, the crack propagation time was found to be 9 \(\mu \)s. Since all of the simulated samples follows the same trend, the crack propagation time was evaluated for all samples. Also, due to some weak grains at the contact point between the sample and SHPB bars, local crushing clusters can be seen (see Fig. 11).

Fig. 12
figure 12

A visualisation of crack propagation through a sample with heterogeneity index \(m=30\) and cement scale factor \(C_f=0.25\). Crack initiation at 76 \(\mu \)s (a), 81 \(\mu \)s (b), fully propagated main crack at 87 \(\mu \)s (c) and post peak at 91 \(\mu \)s (d). The colour overlay represents the bond damage, i.e. the amount of broken bonds for each particle, where no colour represents completely intact bonds

Full size image

An example of crack propagation in a heterogeneous sample with heterogeneity index \(m=30\) and cement scale factor \(C_f=0.25\) can be seen in Fig. 12. Due to the low grain cement strength, a large amount of minor damage is built up in the grain boundaries around the horizontal centre line of the sample, which can be seen at the instant of crack initiation in Fig. 12a. From this built up grain boundary damage, two larger cracks initiate and propagate through the sample. Another consequence of the low cement strength is that the cracks are propagating almost exclusively in the grain boundaries. At the point of crack initiation in Fig. 12a, the sample exhibits only a small amount of crushing at the contacts with the bars. However, the amount of contact crushing increases as the crack has propagated through the sample fully, i.e. at peak stress, see Fig. 12c. This is also when secondary cracks are initiated in grain boundaries around the perimeter of the sample. The fracture pattern after the peak stress can be seen in Fig. 12d, where more secondary perimeter cracks are visible. Using the same approach as is shown in Fig. 11, the crack propagation time is found to be 11 \(\mu \)s for this particular example.

Fig. 13
figure 13

A visualisation of crack propagation through a sample with heterogeneity index \(m=15\) and cement scale factor \(C_f=0.50\). Crack initiation at 78 \(\mu \)s (a), 82 \(\mu \)s (b), fully propagated main crack 88 \(\mu \)s (c) and post peak at 91 \(\mu \)s (d). The colour overlay represents the bond damage, i.e. the amount of broken bonds for each particle, where no colour represents completely intact bonds

Full size image

Another example of crack propagation in a heterogeneous sample can be seen in Fig. 13, here with a heterogeneity index \(m=15\) and cement scale factor \(C_f=0.50\). In contrast to the previous example in Fig. 12, the higher cement strength yields less built up minor damage in the sample before crack initiation. Nevertheless, the main splitting crack can be seen initiating from boundary damage close to the centre of the disc in Fig. 13a, b. Even with a higher cement strength, the crack still appears to prefer grain boundary propagation, although, the crack propagation path can be seen to cross some grains in Fig. 13c. Similarly to the previous case, the amount of crushing increases after peak stress (see Fig. 13d). Compared to the example in Fig. 12, more crushing occurs here due to weaker (darker) grains being located at the contact point between sample and the Hopkinson bars. However, compared to the previous sample, there is a smaller amount of secondary cracks initiating from the sample perimeter. It should be noted that, before the main crack initiation in the sample, a smaller crack has propagated from the incident bar face (see the left part of the sample in Fig. 13a, b). Due to this, the validity of the test can be debated, even though the main crack initiates from the sample centre. By measuring the time between crack initiation and peak stress, the crack propagation time is found to be 10 \(\mu \)s.

Fig. 14
figure 14

A visualisation of crack propagation through a sample with heterogeneity index \(m=5\) and cement scale factor \(C_f=0.75\). Crack initiation at 77 \(\mu \)s (a), 81 \(\mu \)s (b), fully propagated main crack at 84 \(\mu \)s (c) and post peak at 91 \(\mu \)s (d). The colour overlay represents the bond damage, i.e. the amount of broken bonds for each particle, where no colour represents completely intact bonds

Full size image

As a last example of crack propagation in a heterogeneous sample, a sample with heterogeneity index \(m=5\) and cement scale factor \(C_f=0.75\) is shown in Fig. 14. Compared to the previous two examples, there is only a small amount of built up damage in the grain boundaries due to the higher cement strength. Here, the main crack initiates in a grain boundary close to the transmitted bar (see Fig. 14a, b). As the main crack propagates through the sample, another crack initiates close to the centre part of the disc. These two cracks then coalesce to form the final crack path in Fig. 14. With this high cement strength, the crack propagation path crosses more grains compared to before. However, the crack still prefers the path of least resistance through the grain boundaries. Compared to the examples in Figs. 12, 13, the sample in Fig. 14 exhibits significantly less secondary cracks at the perimeter. The crack propagation time in this example is found to be 7 \(\mu \)s.

Fig. 15
figure 15

High-speed photograph (381,000 fps) of two dynamic Brazilian disc test on Kuru black diorite [35]. Example one (a)–(d) and example two (e), (f). Note that the surfaces are painted with paint for digital image correlation measurements and the actual rock surfaces are not visible

Full size image
Table 4 A summary of results from all of the heterogeneous cases summarised alongside their experimental counterparts
Full size table

Generally, the crack path is found to exhibit a more zig-zag pattern for lower cement strengths, which also has been observed using this model for quasi-static conditions [36]. Furthermore, the use of another Weibull distribution for the cement between different grains resulted in a more easily controlled grain boundary. Previously, the grain cementing between two grains was defined as a scaled down version of the strengths of one of the grains being cemented. This means that grains with a significantly low grain strength could result in grain boundaries with extremely low cementing strengths. By comparing the simulated heterogeneous samples in Figs. 12, 14 to the high-speed photography of the real crack propagation in Kuru black diorite, similarities are found. Although it is difficult to determine the exact crack initiation point in the experimental BDT, the example in Fig. 15a–d appears to show two distinct crack initiations that coalesce to form the final splitting crack. Also, secondary, smaller cracks can be seen around the main crack and close to the Hopkinson bars in both examples in Fig. 15. Furthermore, the crushed zones at interfaces between the samples and Hopkinson bars are visible, particularly in the example in Fig. 15e–h.

3.3 Measured results

Fig. 16
figure 16

The tensile stress at the centre of the disc, as predicted from Eq. (1), for all cases of heterogeneity indexes, \(m=\)30, 15, 5, and cement scale factors, \(C_f=\)0.25, 0.50, 0.75. Each simulated case consists of different samples and is compared to five experimental tests. The dotted lines correspond to ± the standard deviation

Full size image

All of the results from measurements on the simulated heterogeneous samples are summarised in Table 4 together with the experimental counterparts [35]. The indirect tensile stress (Eq. 1) versus time for each of the nine cases are shown and compared with data from five experimental measurements of Kuru black diorite in Fig. 16, where solid lines correspond to the average stress at the centre of the disc and the dotted lines denote the ± standard deviations. The peak stress of all of the simulated cases show good correspondence with the experimental data. The experimental stress response in Fig. 16 shows initial variations that the simulations does not exhibit. This is largely due to the fact that the input pulse for each experiment is different from each other (see Fig. 3), while the input pulses for the simulations are identical. However, as the stress increases, the variations are clear. Generally, the variations in peak stress increase with lower m, except for case 4 (\(m=15\), \(C_f=0.25\)) and case 6 (\(m=15\), \(C_f=0.75\)). This is expected since a lower heterogeneity index yields a distribution with larger variations. There is also a trend for an increase in peak strength variations for lower cement strengths, except for case 8 (\(m=5\), \(C_f=0.50\)). This is also an expected result since a lower cement strength yields a sample more sensitive to the grain shape heterogeneity.

For the overloading effect, the study shows a clear dependency on cement strength; the overloading effect decreases with an increase in cement strength. The reason for this is not clear to the authors, but it could be explained by the fact that samples with a lower cement strength will inhabit grain boundary particle bonds with extremely low strengths that will initiate a major crack at a significantly lower stress level. This stress level might be large enough to initiate a crack, but not large enough to make it propagate through the whole sample, which would also result in a lower peak stress. This theory is further supported by the crack propagation time in Table 4, where a clear decreasing trend in propagation time is seen when increasing the cement strength. As a consequence of higher overloading for lower cement strengths, the stress at the point of crack initiation is lower for the cases with lower cement strengths. This means that, even though the peak stresses are approximately the same, the actual indirect tensile strength is significantly different for cases with lower cement strength. The measured indirect tensile strain rate decreases with an increase in cement strength. This is probably due to the fact that samples with lower cement strengths exhibit a larger amount of induced damage before the main crack initiation. This induced damage degrades the sample material, which allows for larger vertical deformation.

By comparing the simulated results to the experimental counterparts in Table 4, it is clear that the model is able to capture the peak stress, overloading, actual tensile strength, strain rate and crack propagation time with reasonable accuracy. However, the experimental crack propagation time is only evaluated from one test.

4 Conclusion

In this study, a recently published heterogeneous bonded particle model was further developed, calibrated and used to study the effects of heterogeneity and cement strength on a dynamic Brazilian disc test. The further development of the model in this study includes the use of three Weibull distributions - one for distributing average mechanical properties for the grains, one for distributing the properties within each grain and one for the cement between different grains. The results were compared to experimental data in terms of crack initiation, propagation as well as indirect tensile stress, strain and strain rate.

  • The heterogeneous bonded particle model is able to capture key phenomena of the dynamic Brazilian disc test, such as the stochastic crack initiation and propagation as well as the magnitude and variations of measured quantities. Also, the results measured in the simulations correspond well with its experimental counterparts.

  • The heterogeneous samples are exhibiting the same overloading phenomenon as seen in experimental dynamic Brazilian disc tests. Furthermore, the importance of compensating for this effect was demonstrated in this study.

  • The results show that the cement strength is affecting several factors of the test, such as the crack propagation path, overloading magnitude and consequently the actual tensile strength. Furthermore, the indirect tensile strain rate is increasing for lower cement strengths, which is suggested to be due to more material degradation for weaker grain boundaries.

  • It was demonstrated how the peak stress corresponds to the point in time when the main splitting crack had propagated fully through the sample. This suggests that the crack propagation time can be estimated as the time between crack initiation and peak stress.

  • Future work using the heterogeneous bonded particle model will include a full calibration, using e.g. the Brazilian disc test in conjunction with the uniaxial compression test. Furthermore, since it is ensured that this model to be able to replicate dynamic rock fracture, this modelling approach will be used to model a full-scale rock drilling process. This would further require a loading rate dependent BPM in order to account for a wide range of strain rates.