Hybrid finite-discrete element modelling of rock fracture during conventional compressive and tensile strength tests under quasi-static and dynamic loading conditions

Abstract

A hybrid finite-discrete element method is proposed to model the rock fracture under various loading conditions. The key component of the hybrid method, i.e. transition from continuum to discontinuum through fracture and fragmentation, is introduced in detail. An empirical relationship between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments is implemented in the hybrid method to model the effect of loading rate. The hybrid method is calibrated by modelling the Brazilian tensile strength (BTS) test, uniaxial compressive strength (UCS) and notched Brazilian dis tests. Then the hybrid method is employed to model the dynamic rock fracture process in UCS and BTS tests. The proposed method has well modelled the dynamic rock fracture and fragmentation processes and captured the effect of loading rate on rock strengths. It is concluded that the hybrid finite-discrete element is a valuable tool to study the dynamic rock fracture as it takes the advantages of the continuum and discontinuum based method, and considers the effect of loading rate.

Keywords:
Hybrid finite-discrete element; transition from continuum to discontinuum; rock fracture; fracture toughness; loading rate

Graphical Abstract

1 INTRODUCTION

The fracture, failure and collapse of rock-like materials have been studied widely since it is significant to understand the fracture mechanism in civil engineering and mining engineering and several other files in which rock fractures play an important role (An and Liu 2014, Liu and An 2016, An, Liu et al. 2017a, An, Liu et al. 2017b, An, Hou et al. 2019, Fukuda, Mohammadnejad et al. 2019a). However, it is still far more from a better understanding of the rock mechanism, especially the rock under dynamic loading. It is imperative to study the dynamic rock mechanism as it not only can improve the efficiency of rock breaking and structure demolition but it also can help to prevent the collapse of geo-structure and increase their stabilities. With the fast development of computer technology, numerical methods have been widely developed to model the complex dynamic rock fracture process. (Zhu and Tang 2006, Liu, Kou et al. 2007, Stefanizzi, Barla et al. 2009, Mahabadi, Cottrell et al. 2010).In general, numerical method can be classified as continuum based method and discontinuum based method according to the hypothesis of the rock mass as continuous or discontinuous material. In addition, the coupled continuum-discontinuum method is a combination of the above two type methods by treating the materials as partly continuum or discontinuum but lack of transition from the continuum to discontinuum. However, for a realistic modelling of the rock fracture process, the numerical method should be able to model the rock mass as a continuous body and the interaction of the rock mass as a number of discrete bodies and the transition from a continuous body to a number of discontinuous bodies. Thus, the hybrid continuum-discontinuum method has recently attracted the attention of many engineers and researchers, which is able to model the transition from continuum to discontinuum through rock fracture and fragmentation (An, Liu et al. 2017a, Fukuda, Mohammadnejad et al. 2019a, 2019b). The hybrid finite-discrete element method (FDEM) proposed by Munjiza (2004) may be the most widely used hybrid continuum-discontinuum method. It takes the advantages of continuum method and discontinuum method and can model the transition process of rock from continuum to discontinuum trough rock fracture and fragmentation (An and Liu 2014, An, Liu et al. 2017a). As an open source, Y code is one of the main implementations of the hybrid FDEM originally proposed by Munjiza (2004). Some extension of the Y code includes Y-Geo (Mahabadi, Lisjak et al. 2012), Y-Flow (Lisjak, Mahabadi et al. 2018). The authors also have developed the Y-HFDEM IDE(Liu, Kang et al. 2015, An, Liu et al. 2017a).

In this paper, the hybrid finite-discrete element method is used to model the rock fracture process of the conventional compressive and tensile strength tests under quasi-static and dynamic loading. The obtained results of the quasi-static loading tests are compared with the experimental results and those well documented in the literatures to calibrate the proposed method and to show its abilities of modelling the transition from continuum to discontinuum through fracture and fragmentation. The modelling of the dynamic rock fracture process is used to show the ability of the proposed method to capture the rock characteristics on loading rates during the dynamic rock fracture process.

2 HYBRID FINITE-DISCRETE ELEMENT METHOD

The hybrid finite-discrete element mode is consist of one or a number of discrete bodies, which are then discretised into three-node finite elements. Among the finite element edges, four-node joint elements are embedded. The deformity of the discrete bodies is analysed by the finite elements while the fractures are modelled using the four-node joint elements depending on the calculated stress and strain of the discretised finite elements in the discrete bodies. Thus the hybrid finite-discrete element is superior to the continuum based finite element method and the discontinuum based discrete element method as it takes the advantages of the both two mentioned methods. The hybrid finite-discrete element method consists of the following components: contact detection, contact interaction between individual bodies, deformability and transition from continuum to discontinuum, temporal integration scheme. The contact detection component is used to detect the couples of the discrete bodies, which can be in contact, and eliminate the discrete bodies, which cannot be in contact. The transition from continuum to discontinuum is the key component, which makes the rock fracture process possible through a cohesive fracture model. The temporal integration scheme is applied to solve the equations of motion for the discretised bodies and to update the nodal coordinates at the time step.

In this section, the key component of the hybrid finite-discrete element method, i.e. the transition from continuum to discontinuum, is introduced in detail. In terms of modelling the rock failure process under various loading rates, the loading rate dependent should be taken into account. Therefore, an empirical relationship between the physical-mechanical properties of rock and the loading rate is implemented into the hybrid finite-discrete element method, which is also introduced in the section.

2.1 Transition from continuum to discontinuum through fracture and fragmentation

For a typical stress-strain curve of rock-like materials under dynamic loading, it can be divided into stress-hardening part (before the peak stress is reached) and stress-softening part, which presents decreasing stress with the increasing strain (Munjiza 2004). The stress-hardening part is implemented in a standard way through the constitutive law in the hybrid finite-discrete element method. The stress softening part of the stress-strain curve is related to the damage and fracture of the rock, i.e. transition from continuum to discontinuum. Formulations of the strain softening by means of stress and displacements are adopted to deal with the stress softening part of the stress-strain curve. As illustrated in Figure 1, the three nodes finite-element is used to model the deformity of the rock, i.e. the stress hardening part of the stress-strain curve, while the four-node joint element is used to model the fracture initiation and propagation depending on the calculated stress and strain of the discretised finite elements, i.e. the stress softening part of the stress-strain curve. In the hybrid finite-discrete element method, the fracture of the discrete bodies is assumed to occur at the surface of the finite element edges. The fracture initiation and propagation is implemented by the distortion of the joint elements, which involves a bonding stress. Figure 1 illustrates the distortion of the joint elements while the discrete body is under tensile, shear and the combined stress conditions, respectively. Thus, according to the stress condition or the distortion of the joint elements, the different fracture models, i.e. pure mode-I, pure mode-II and mixed-model I-II, can be modelled.

Figure 1
Hybrind finite-discrete element model under different stress conditions (red line represents edges of joint element).

At any point of the element edges, the separation of the adjacent finite element surface or the distortion of the joint element $\delta$ can be divided into normal direction and tangential direction as illustrated in Figure 1d.

$\delta ={\delta }_{n}n+{\delta }_{s}t$ (1)

where $n$ and $t$ are the unit vectors in the normal and tangential direction, respectively, of the surface at such a point, ${\delta }_n$and ${\delta }_s$ are the magnitudes of the components of$\delta$.

As illustrated in Figure 2, the relationship of the bonding stress and the separation of the adjacent finite element for pure mode-I fracture (tensile failure) indicate in the right side of the figure. For the pure mode-I fracture, as the tensile strength increases, the separation ${\delta }_n$ of the adjacent finite element in the form of the distortion of the joint element in normal direction (Figure 1 b) increases. As the separation ${\delta }_n$ reaches a critical displacement${\delta }_{np}$, prescribed by the tensile strength of the material ${\sigma }_n$ , the tensile failure occurs (i.e. mode-I fracture ) as illustrated in the right side of Figure 2. As the separation of the adjacent element continue increase, the bonding stress decrease gradually. While the separation exceeds the ultimate opening ${\delta }_{nu}$, determined by the mode-I fracture energy release rate $G_{fI}$, the adjacent finite elements are separated and the mode-I fracture completed. The mode-I fracture energy release rate $G_{fI}$ is equal to the area under the curve of the bonding stress and opening displacement as shown in the right side of Figure 2, and can be expressed by Equation 2.

$G_{fI}=\underset{{\delta }_{np}}{\overset{{\delta }_{nu}}{{\displaystyle \int }}}{\sigma }_n\left({\delta }_n\right)d{\delta }_n$ (2)

During the mode-I fracture process, the bonding stress in the normal direction can be calculated using Equation 3.

${\sigma }_n=\{\begin{array}{ccc}\left[2\cdot \frac{{\delta }_n}{{\delta }_{np}}-{\left(\frac{{\delta }_n}{{\delta }_{np}}\right)}^2\right]\cdot {\sigma }_t& if& 0\le {\delta }_n\le {\delta }_{np}\\ f\left(D\right)\cdot {\sigma }_t& if& {\delta }_{np}\le {\delta }_n\le {\delta }_{nu}\\ 0& if& {\delta }_n\ge {\delta }_{nu}\end{array}$ (3)

where D is a damage variable between 0 and 1,$f\left(D\right)$ is the damage function described in the mechanical damage model(Liu, Kou et al. 2004), and all other parameters have the same meaning as those introduced above.

Figure 2
Relationship between the bonding stress and opening/sliding displacement under tension and shear condition

For the pure mode-II fracture, the relationship of the bonding stress and sliding separation of the adjacent finite elements in the form of the joint elements distortion in the tangential direction (Figure 1c) is demonstrated in the left side of Figure 2. As the sliding displacement ${\delta }_s$ in tangential direction increases, the bonding stress increases, while the sliding displacement reaches a critical slip ${\delta }_{sp}$, described by the shear strength, the shear failure (pure mode-II fracture) occurs. As the sliding displacement continues to increase, i.e. ${\delta }_s>{\delta }_{sp}$, the bonding stress decreases gradually till residual stress according to a mechanical damage model. After the sliding displacement in the tangential direction is beyond the residual sliding displacement ${\delta }_{sr}$ determined by the model II fracture energy release rate$G_{fII}$, the shear failure has completed. However, a purely frictional resistance defined by Columb's model will still exist. The model II fracture energy release rate$G_{fII}$ can be described by the area under the curve of the bending stress and the sliding displacement and can be expressed in Equation 4

$G_{fII}=\underset{{\delta }_{sp}}{\overset{{\delta }_{sr}}{{\displaystyle \int }}}[\tau \left({\delta }_s\right)-{\tau }_r]d{\delta }_s$ (4)

During the pure mode-II fracture process, the bonding stress in the tangential direction can be calculated using Equation 5.

$\tau =\{\begin{array}{ccc}2\cdot \frac{{\delta }_s}{{\delta }_{sp}}\cdot {\sigma }_c& if& 0\le {\delta }_s\le {\delta }_{sp}\\ g\left(D\right)& if& {\delta }_{sp}\le {\delta }_s\le {\delta }_{sr}\\ {\sigma }_n\cdot \tan \left({\varnothing }_f\right)& if& {\delta }_s\ge {\delta }_{sr}\end{array}$ (5)

where $D$ is a damage variable between 0 and 1, $g\left(D\right)$ is damage functions described in the mechanical damage model (Liu, Kou et al. 2004), and ${\varnothing }_f$ is the joint residual friction angle.

For the rock fracture process under dynamic loading, the joint element may under the combination of the tensile stress and shear stress condition, which results in the separation of the finite elements in both normal and tangential direction, i.e. mixed-mode I-II fracture. Figure 3 illustrates the failure criteria for mixed-mode I-II. If Equation 6 is satisfied, the mixed-mode I-II fracture occurs. The fracture energy release rate $G_{fI-II}$ for the mixed-mode I-II equates to the shadow in Figure 3.

${\left(\frac{{\delta }_n-{\delta }_{np}}{{\delta }_{nu}-{\delta }_{np}}\right)}^2+{\left(\frac{{\delta }_s-{\delta }_{sp}}{{\delta }_{sr}-{\delta }_{sp}}\right)}^2\ge 1$ (6)

Figure 3
Failure criteria for Mixed I-II fracture mode.

2.2. Effect of loading rate on the dynamic behaviour of rock

The loading rate plays a significant role during the rock fracture process under dynamic loading (Zhang, Kou et al. 1999, Zhang, Kou et al. 2000, Zhao 2000). In the hybrid finite-discrete element modelling of the dynamic rock fracture process using various test techniques, the effect of loading rate is taken into account by through implementing an empirical relation between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments (Zhao 2000). Zhao(2000) conducted the uniaxial and triaxial compression, uniaxial tension and unconfined shear tests on Bukit Timah granite of Singapore and proposed a semi-log formula in Equation 7. The relationship estimates the dynamic rock strengths under different loading rates by comparing the dynamic loading rate and static loading rate.

${\sigma }_{cd}=A\cdot \log \left(\frac{{\dot{\sigma }}_{cd}}{{\dot{\sigma }}_c}\right)+{\sigma }_c$ (7)

where ${\sigma }_{cd}$ is the dynamic uniaxial compressive strength (MPa), ${\dot{\sigma }}_{cd}$ is the dynamic loading rate (MPa/s), ${\dot{\sigma }}_c$ is the quasi-static loading rate (approximately $5\times {10}^{-2}$ MPa/s), ${\sigma }_c$ is the uniaxial compressive strength at the quasi-static loading rate (MPa) and $A$ is a material parameter, which is 11.9 for the Bukit Timah granite (Zhao 2000).

3 EXPERIMENTAL STUDY OF THE ROCK FRACTURE UNDER QUASI-STATIC LOADING

Before modelling the rock fracture processes, UCS and BTS tests have been done in laboratory to obtain the rock fracture patterns and rock prosperity. Figure 4a shows the rock samples for uniaxial compressive tests, while Figure 4b [[Q5: Q5]] illustrates the rock sample under test. The diameters of the sample for UCS test is 50mm while the length is 100mm. Figure 4c depicts the rock samples with the diameter of 50mm and the height of 25mm for BTS test, while Figure 4d illustrates the rock samples under test. Figure 5 illustrates the rock fracture pattern for UCS test (Figure 5a) and BTS test (Figure 5b), which will then be compared with the modelled result to calibrate the proposed method. Table 1 lists the rock parameters obtained by a serial of laboratory tests, which are then used as input values for modelling rock fracture process by the proposed method. The details of the experimental tests can be found in the first author dissertation(Huaming 2018).

Figure 4
Rock sample for UCS test and BTS test: (a) Rock specimen for UCS test; (b) Rock specimen under UCS test; (c) Rock specimen for BTS test; (d) Rock specimen under BTS test.

Figure 5
Rock fracture patterns in UCS and BTS tests: (a) UCS tests fracture pattern; (b) BTS tests fracture pattern

Table 1
Rock parameters obtained from laboratory tests

Properties	Symbols	Values	Units
Young′s modulus	$E$	59	$GPa$
Poisson′s ration	$\nu$	0.23	$N/A$
Density	$\rho$	2595	$Kg{m}^{-3}$
Tensile strength	${\sigma }_t$	10.86	$MPa$
Compressive Strength	${\sigma }_c$	103	$MPa$
Internal friction coefficient	$\varnothing$	30	$°$(Degree)
Surface friction coefficient	$u$	0.1	$N/A$

4 NUMERICAL MODELLING OF ROCK FRACTURE UNDER QUASI-STATIC LOADING

In this section, a serials of rock fracture tests, i.e. uniaxial compressive strength (UCS) test, Brazilian tensile strength (BTS) test and notched Brazilian disc (NBD) tests, under quasi-static loading are modelled to calibrate the proposed method and demonstrate the capability of the hybrid method in modelling the rock fracture initiation and propagation.

4.1 Numerical modelling of quasi-static rock fracture process during UCS tests

The UCS test is simplified as a plane stress problems and only the vertical section is taken into account. As illustrated in Figure 6a, the geometrical model is built according to the geometry suggested by ISRM for rock characterization, testing and monitoring 2007-2014 (Ulusay 2015), i.e. 50mm for the width and 100 for the length of the sample, while the numerical model is discretized by triangle elements. The materials properties of the rock specimen follow those in Table 1. For hybrid finite-discrete element method, rock fracture energy release rates are critical parameters. Those parameters are not directly obtained in our laboratory tests. However, it can be estimated by the following equations:

$G_{fI}=\frac{K_{Ic}{}^2}{E}$ (8)

$G_{fII}=\frac{\left(1-v^2\right)K_{IIc}{}^2}{E}$ (9)

where $G_{fI}$ and $G_{fII}$ are the mode-I and mode-II fracture energy release rates, $K_{Ic}$ and $K_{IIc}$ are the mode-I and mode-II fracture toughness. For mode-I fracture toughness, it can be calculated according to Equation 10, which is used the tensile strength to estimate the rock fracture toughness(Zhang 2002).

${\sigma }_t=6.88K_{Ic}$ (10)

The mode-II rock fracture toughness is assumed to be 1.7 times larger than the mode-I rock fracture toughness (Ingraffea 1981, Whittaker, Singh et al. 1992, Al-Shayea, Khan et al. 2000, Khan and Al-Shayea 2000). Thus all the rock parameters for hybrid modelling are obtained on the basis of Equation 8 to 10.

The material properties of the loading plates follow those of standard steels, while the penalty terms of the rock specimen and loading plates are assumed to be equal to the elastic modulus of the rock specimen and half of that the steel, respectively.

During the test, a constant displacement increment of 0.1m/s is applied on the top plate in the vertical direction, while the bottom plate is fixed both in the vertical and horizontal directions. It should be noted that the loading rate is much higher than those used in the laboratory test. Nevertheless, it should not influence the rock properties significantly as these are still much smaller than the certain threshold according to the laboratory test (Zhang 2001).

Figure 6
Geometrical and numerical model of UCS test: (a) Geometrical model; (b) Numerical model

Figure 7a depicts the evolution of the minor principal stress while Figure 7b shows the fracture initiation and propagation during the UCS test under quasi-static load, i.e. 0.1m/s. Figure 8 records the relationships between forces from the two plates and the displacements. In Figure 7a, the magnitudes of the minor principal stress can be referred to the legend on the top of Figure 7a. Following the solid mechanics conventions, the compressive stresses are taken as negative, while the tensile stresses are regarded as positive. In Figure 7b the tensile failure is marked using red colour while the shear failure and the boundaries are denoted using a blue colour.

It can be seen from Figure 7a-A that as the top-loading plate contacts the specimen, a high-stress concentration is induced and propagate from the top of the specimen to the bottom. As the top plate further moves, the stress continues to propagate and reaches the bottom of the specimen (Figure 7a-B). Then the stress reflected from the contact of the bottom of the specimen and the fixed bottom plate (Figure 7a-C). After the stress travels a few times between the top and bottom plates, the relatively uniform stress distribution is formed (Figure 7a-E and F).

During this period, the force-loading increases dramatically and reaches its peak when the palate has moved 23$\mu$m (Figure 8E and F). Meanwhile, small cracks initiate along an approximate 45-degree plane (Figure 7a-E and F). As the top plate continues to move, the force-loadings drop dramatically from the peak to the bottom, which indicates that the specimen loses its ability of carrying loads. Then stress fluctuates at a very low level (Figure 8a-H, I and L). During this process, more cracks are produced along an approximate 45-degree plane which are mainly shear cracks, while cracks along the vertical direction are observed which indices part tensile cracks and part shear cracks. The final modelled result (Figure 7b-L) agrees well with the experimental result (Figure 5a). The main cracks in the both the experimental and numerical results are vertical cracks along the loading line and inclined cracks along approximately 45 degree plane. Thus, the Hybrid finite-discrete element method can well reproduce the rock fracture process.

The modelled results reproduce the brittle rock failure process under loading and agree with laboratory results (Wawersik and Fairhurst 1970): a nearly linear loading increase region (Fig.8-AE), a big jump from the peak value to bottom (Figure 8-FI) and a post-failure region (Figure 8-IJ). According to the stress-displacement curve, the compressive strength of the specimen is approximately 106MPa, which is only a little bit larger than the input value (103MPa).

Figure 7
Modelled failure process during uniaxial compression test under the loading rate of constant displacement of 0.1m/s: (a) Evolution of minor principal stress in the vertical direction; (b) Crack initiation and propagation (red: shear failure; Blue: tensile failure)

Figure 8
Stress-loading displacement curve during the uniaxial compression test under a loading rate of 0.1m/s

4.2 Modelling of quasi-static rock fracture process during BTS test

BTS test is developed to determine the tensile strength of brittle materials indirectly using cylindrical specimen since conventional methods of tensile test present difficulties to cope with low tensile resistance materials, such as concrete, rock, and rock-like materials. Since Akazawa (1943) and Carneiro (1943) independently developed the BTS test at almost the same time, the BTS test has gained his popularity for calculating the tensile strength. Then Hondros (1959) gave a complete stress solution for disc BTS test valid for both plane stress and plane strain condition.

In this section, hybrid finite-discrete element method is employed to model the quasi-static and dynamic rock failure processes in BTS test to calculate the tensile strength under quasi-static loading and to capture the fracture initiation and propagation process. Figure 9 illustrates the geometrical and numerical models. The numerical model is built according to the geometry suggested for rock characterization, testing and monitoring 2007-2014 by ISRM(Ulusay 2015), i.e. the diameter of the specimen is 50mm as shown in Figure 9a. The specimen is placed between two plates which move towards each other with a constant displacement increment during BTS test. The loads from the plates are assumed to be radially applied over a short strip of the circumference with a radian of$2\alpha$ as shown in Figure9a. Then the numerical model is discretized by triangular elements as shown in Figure 9b while the material properties follow those in the UCS test.

Thus, taking into account the specimen thickness $t$, the tensile strength can be calculated by Equation 11 (Timoshenko and Goodier 1970):

${\sigma }_t=\frac{2P}{\pi Dt}$ (11)

where the ${\sigma }_t$ is the tensile strength, $P$ is the applied load, $R$ is the diameter.

The complete stress solution along the load diameter given by Hondros (1959) can be described as follows.

${\sigma }_x=\frac{P}{\pi Rt\alpha }\left\{\frac{\left[1-{\left(\frac{r}{R}\right)}^2\right]\mathrm{Sin}2\alpha }{1-2{\left(\frac{r}{R}\right)}^2\mathrm{Cos}2\alpha +{\left(\frac{r}{R}\right)}^4}-{\tan }^{-1}\left[\frac{1+{\left(\frac{r}{R}\right)}^2}{1-{\left(\frac{r}{R}\right)}^2}\tan \left(\alpha \right)\right]\right\}$ (12)

${\sigma }_y=-\frac{P}{\pi Rt\alpha }\left\{\frac{\left[1-{\left(\frac{r}{R}\right)}^2\right]\mathrm{Sin}2\alpha }{1-2{\left(\frac{r}{R}\right)}^2\mathrm{Cos}2\alpha +{\left(\frac{r}{R}\right)}^4}+{\tan }^{-1}\left[\frac{1+{\left(\frac{r}{R}\right)}^2}{1-{\left(\frac{r}{R}\right)}^2}\tan \alpha \right]\right\}$ (13)

where$r$ is the distance from the entre of the disc, 2$\alpha$ is radian of contact between plates and the specimen, ${\sigma }_{xx}$ and ${\sigma }_{yy}$ are stresses along the horizontal and vertical directions respectively.

Figure 9
Geometrical and numerical models for the Brazilian tensile disc test: $P$ is the applied load, $R$ is the disc radius,$r$ is the distance from the center of the disc, $t$ is the disc thickness, 2$\alpha$ is the angular distance of load arc, ${\sigma }_{xx}$ and ${\sigma }_{yy}$ are stresses along the horizontal and vertical directions respectively.

During modelling the quasi-static rock failure process in BTS test, a constant displacement increment of 0.1m/s is applied on both the top and bottom plates in the vertical direction while the two plates are fixed in the horizontal direction. Figure 10a visually illustrates the temporal and spatial rock failure process for BTS test under the loading rate of 0.1m/s while Figure 11 demonstrates the corresponding force-loading displacement curves.

As the plates contact the specimen, high-stress concentrations are induced immediately at the contact areas (Figure 10a-A). Then the stresses transmit into the specimen and mainly propagate along the loading line (Figure 10a-BC). As the loading plates continue move, a rather uniform stress distribution along the loading diameter is formed (Figure 10a-CD). During this period, the Force-loadings increase almost linearly and peak at point D as shown in Figure 11a-A-D and Figure 11b-A-D. Then a tensile crack is induced (Figure 10b-E) and propagate along the loading diameter (Figure 10b-F).It should be note that, in the BTS test, the tensile crack is marked using the red color while the shear crack and boundaries are denoted using the blue color. As the tensile crack propagates toward the loading areas, shear cracks are induced due to the compressive concentration at the loading vicinities (Figure 10b-J). During this period, the forces drop dramatically from the peak points to the bottom of the curves (Figure 11a-D-F and Figure 11b-D-F), which means that the specimen gradually losses its load bearing capacity. Finally more tensile and shear cracks are induced and the specimen is split into two halves due to the tension in the horizontal direction (Figure 10b-K).

The modelled fore-loading displacement curves of the BTS recorded in Figure 11 present the typical behavior of brittle rock under compression: a compressive deformation region (AB), a linear-elastic deformation region (BC), a non-linear deformation region (CD), a strain-softening deformation region (DEF) and a residual deformation region (FG). According to the maximum load $P_{Max}$ at point D and the Equation 10, the tensile strength of the specimen can be calculated as follows.

${\sigma }_t=\frac{2P}{\pi Dt}=\frac{2\times 1.01\times {10}^6}{3.14\times 0.050\times 1}=12.87MPa$ (14)

The modelled result is larger than the obtained result from the laboratory test (10.86Mpa), and the variance is (12.86-10.86)/12.86=15.5. The variance is caused by the adopting a larger loading rate (0.1m/s) than that in the laboratory test (around 0.05mm/s).

Figure 12 compares the modelled rock fracture pattern (Figure 12a) with the experimental result (Figure 12b). It can be seen that the modelled result agrees well with the experimental result as they both fracture along the vertical diameter and also some fragments are produced at the top loading vicinities. In addition, Li and Wong (2013), and Hobbs (1964) indicate that a typical failure pattern of BTS test includes tensile failure along the loading diameter and shear failure at the top and bottom loading areas. Thus, the hybrid method can well reproduce the rock failure process in the BTS test under static loading.

Figure 13 depicts the theoretical and numerical stress distribution along the longing diameter, which are normalized by 2P/πDt for easy comparison at the time when stress equilibrium reaches. The theoretical stresses distribution are calculated according to Equation 12 and 13 (where $2\alpha$ = 0.1047, i.e.$6^0$, t=1m, D=0.054m) (Hondros 1959).It can be seen that the modelled stress distributions along the loading diameter agree well with the theoretical solution as the two curves are almost overlapped each other. However, at the loading points, the stress curves are not agreeable well, which is caused by the shear cracks occurring at the loading points.

Figure 10
Evaluation of minor principal stress and the fracture initiation and propagation for BTS test under the loading rate of 0.1m/s: (a) Evolution of minor principal stress; (b) Fracture initiation and propagation.

Figure 11
Force-loading displacement curve for during BTS test unde a quasi-static load (0.1m/s).(a) Top force-loading displacement curve; (b) Bottom force-loading displacement curve

Figure 12
Comparison of the modelled and experimental fracture pattern for BTS test: (a) modelled result; (b) experimental results.

Figure 13
Comparison of the numerical and analytical stress along the loading diameter for BTS test under the loading rate of 0.1m/s

4.3 Modelling of quasi-static rock fracture process during NBD tests

In this section, the NBD test is used to model rock fracture during the NBD tests under quasi-static loading. The objective of this section is to investigate the fracture initiation and propagation of the pure mode-I, pure mode-II and the mixed-mode I-II fracture processes and obtained the corresponding fracture toughness.

Figure 14 illustrates the geometrical model for the NBD test. The geometrical model is the same as the BTS test expect that NBD test model has a pre-fabricated notch with various orientation angles. By changing the angles of the notch with respect to the direction of the loading, the mode-I, mode-II and the mix-mode I-II can be obtained. The analytical solution for the pure mode-I and pure mode-II stress intensity factors (SIFs) can be calculated according to the following equations

Figure 14
Geometrical mode of notched Brazilian disc: a is the half crack length; R is the diameter; and $\theta$ is the orientation angle

$K_I=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_1$ (15)

$K_{II}=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_2$ (16)

where $K_I$ is the model I stress intensity factor, $K_{II}$ is the mode-II stress intensity factor, R is the radius of the Brazilian disc, B is the thickness of the disc, P is the peak load at failure, and $N_1$ and $N_2$ are non-dimensional coefficients and can be defined as follows.

$N_1=1-4{\sin }^2\theta +4{\sin }^2\theta \left(1-4{\cos }^2\theta \right){\left(\frac{a}{R}\right)}^2$ (17)

$N_2=\left[2+\left(8{\cos }^2\theta -5\right){\left(\begin{array}{c}a\\ R\end{array}\right)}^2\right]\sin 2\theta$ (18)

where $\theta$ is the orientation angle.

$N_1$ and $N_2$ are related to the half crack (pre-fabricated notch) $a$, the diameter $R$ and the crack orientation angle, $\theta$. For the specific value of a/R=0.3, the variation of the $N_1$ and $N_2$ with the change of orientation angles can be described in Figure 15.

Figure 15
Stress intensity factor(SIF) in diametral compression of the NBD at a=15.4mm ,R=54mm,a/R = 15.4/54(after Whittaker et al., 1992)(Barry, Raghu et al. 1992)

On the basis of Equation 17 and Equation 18, the oriental angles for the pure mode-I and pure mode-II fractures can be obtained by letting $N_2=0$ and $N_1=0$, respectively. For the fixed a/R=0.3, the pure mode-I fracture can be obtained at the oriental angle of $0^0$, while pure mode-II fracture can be obtained at the oriental angle of ${27.5}^0$. In addition, the oriental angle of ${45}^0$ is selected to obtain the mixed-mode I-II fracture.

Figure 16 illustrates the geometrical and numerical models of NBD test for modelling of pure mode-I, pure mode-II and mixed-mode I-II fractures. The geometrical model is the same as that for the BTS test, except a notch of 15mm in length is pre-fabricated in the specimen as shown in Figure 16a. The numerical model is discretized using finite element and the material parameters are the same as used for UCS and BTS tests. During the test, a constant displacement increment of 0.1m/s is applied on the top and bottom loading plates in the vertical direction while the two plates are fixed in the horizontal direction.

Figure 16
Geometrical and numerical model for diametral compression of the notched Brazilian Dis tests: (a) Geometrical model; (b) numerical model; Angles for model A, B, and C are$0^0$,${27.5}^{0}and{45}^0$ respectively.

Figure 17 illustrates the rock fracture initiation and propagation process of NBD test for modelling the pure mode-I fracture. As shown in Figure 17, during the diametral compression of the top and bottom plates, the cracks firstly initiate from the two tips of the pre-fabricated notch (Figure 17B). Then they propagate along the loading diameter to the loading top and bottom points (Figure 17C). As the two loading plates continue to move, the stresses concentrate on the top and bottom loading contacts. As the stress concentration meet the shear failure crania, shear failures occurs (Figure 17D). Finally, the disc is split into two halves along the loading diameter and some fragments are produced at the top and bottom loading areas (Figure 17F).

Figure 18 shows the numerical and experimental results in diametral compression of NBD test under pure mode-I loading from the literatures. The hybrid finite-discrete element method modelled fracture pattern (Figure 17 F) show a good agreement with the experimental and numerical results in literatures: the cracks initiate around the tips of the pre-fabricated notch and then propagate almost parallel to the loading diameter. Finally, the cracks approach to the loading points and some fragments might be produced depending on the loading duration and the loading rate. Thus the proposed hybrid method can well simulate the pure mode-I fracture initiation and the subsequent crack propagation during the NBD test under diametral loading.

The obtained force-loading displacement curve of NBD test in pure mode-I loading is recorded in Figure 19, which is similar to those recorded in the UCS or BTS test. The alphabets in Figure 19 correspond to those in Figure 17. The Region AB of the force-loading displacement curve is considered as nearly linear elastic region, and the deformation can be recoverable before the cracks occur. Region BD is regarded as a nonlinear elastic region and the fractures propagate steadily. As for the Region CD, it is an unstable fracture propagation region. The last Region DEF, is the post failure region. During the region, the force-loading decreased sharply and the specimen lose the bearing capability completely.

Figure 17
Rock fracture initiation and propagation for pure mode-I NBD test under quasi-static loading (red color represents tensile failure while the blue color indicates the shear failure, pre-fabricated notch and the boundaries)

Figure 18
Numerical and experimental results obtained in the diametral compression of NBD tests in pure mode-I loading: i) Simulated by the $R-T^{2D}$ code(Liu, Kou et al. 2007) ii) Experimental result (Jia, Castro-Montero et al. 1996); iii) Experimental and numerical result(Chen, Pan et al. 1998)

According to the peak load at point D, the pure mode-I fracture toughness can be calculated on the basis of the aforementioned method. Equation 19 and Equation 20 are used to calculate the mode-I fracture toughness.

$N_1=1-4{\sin }^{2}0+4{\sin }^{2}0\left(1-4{\cos }^{2}0\right){\left(\frac{0.0075}{0.025}\right)}^2=1$ (19)

$K_I=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_1=\frac{0.95\times {10}^6\times \sqrt{7.5\times {10}^{-3}}}{\sqrt{3.14\times }0.025\times 1}\times 1=1.857MPa\sqrt{m}$ (20)

Figure 19
Force loading-displacement curve for pure mode-I NBN test

Figure 20 demonstrates the rock fracture process for the pure mode-II NBD test. As can be seen that from Figure 20, the fractures initiate from the tips of the pre-fabricated notch (Figure 20B) and propagate to the top and bottom loading points (Figure 20B). Finally, the initiated fractures approached the loading points in zigzag paths (Figure 20D). In addition, two fractures initiating from the rim of the specimen and propagating along the diameter through the pre-fabricated notch are also observed.

Figure 21 illustrates the numerical and experimental results for the pure mode-II NBD test in literatures. The hybrid modelled result (Figure 20D) shows a good agreement with those in the literatures (Figure 21). The fracture for those tests all initiate from the tips of the loading points and approach to the loading points. Figure 22 shows the force-loading displacement curve for pure mode-II NBD test, which is similar to the UCS, BTS tests. The curve represents typical brittle failure characteristics with linear elastic region AB, nonlinear elastic region BC and post-failure region CE. According to the peak force at point C, the pure mode-II fracture toughness can be calculated as follows.

$N_2=\left[2+\left(8{\cos }^2\theta -5\right){\left(\frac{a}{R}\right)}^2\right]\sin 2\theta =\left[2+\left(8{\cos }^2{27.5}^0-5\right){\left(\frac{0.0075}{0.027}\right)}^2\right]\sin (2\times {27.5}^0)=1.72405$ (21)

$0)K_{II}=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_2=\frac{0.87\times {10}^6\times \sqrt{0.0075}}{\sqrt{\pi 3.14}\times 0.025\times 1}\times 1.724052=2.93MPa\sqrt{m}$ (22)

Figure 20
Rock fracture initiation and propagation for pure mode-II NBD test under quasi-static loading

Figure 21
Numerical and experimental results obtained in the diametral compression of NBD tests in pure mode-II loading: i) Simulated by the $R-T^{2D}$ code (Liu, Kou et al. 2007) ii) Experimental result (Jia, Castro-Montero et al. 1996); iii) Experimental and numerical result (Chen, Pan et al. 1998)

Figure 22
Force loading-displacement curve for pure mode-II NBN test

Figure 23 shows the rock fracture process for the mixed-mode I- II NBD test while Figure 24 shows the corresponding force-loading displacement curve. The fracture initiation and propagation process and the final fracture pattern are similar to that for the pure mode-II test. The fracture initiates from the tips of the pre-fabricated notch, propagate in a zigzag path and finally approach to the top and bottom loading points. Figure 24 show the force-loading displacement curve for mixed-mode I-II NBD test. The curve shows a typical failure process of brittle materials. According to the aforementioned theory and the peal forces at point B of the force-loading displacement cure (Figure 24), the mixed-mode fracture toughness can be calculated as follows.

$\begin{array}{l}{N}_1=1-4{\sin }^2\theta +4{\sin }^2\theta \left(1-4{\cos }^2\theta \right)\times {\left(\frac{a}{R}\right)}^2\\ =1-4{\sin }^2{45}^0+4{\sin }^2{45}^0\times \left(1-4{\cos }^2{45}^0\right){\left(\frac{0.0075}{0.025}\right)}^2\\ =-1.16\end{array}$ (23)

$\begin{array}{l}{N}_2=\left[2+\left(8{\cos }^2\theta -5\right){\left(\frac{a}{R}\right)}^2\right]\sin 2\theta \\ =\left[2+\left(8{\cos }^2{45}^0-5\right)\times {\left(\frac{0.0075}{0.025}\right)}^2\right]\sin (2\times {45}^0)\\ =1.92\end{array}$ (24)

Figure 23
Rock fracture initiation and propagation for pure mixed-mode I-II NBD test under quasi-static loading (red color represents tensile failure while the blue color indicates the shear failure, pre-fabricated notch and the boundaries)

Figure 24
Force loading-displacement curve for mixed-mode I-II NBN test

$K_{I-Mix}=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_1=\frac{0.8\times {10}^6\times \sqrt{0.0075}}{\sqrt{3.14}\times 0.025\times 1}\times 1.16=1.81MPa\sqrt{m}$ (25)

$K_{II-Mix}=\frac{P\sqrt{a}}{\sqrt{\pi }RB}{N}_2=\frac{0.8\times {10}^6\times \sqrt{0.0075}}{\sqrt{3.14}\times 0.025\times 1}\times 1.92=3MPa\sqrt{m}$ (26)

The fracture toughness obtained from the mixed-mode I-II modelling is 1.69 for the model I fracture and 2.82 for the mode-II fracture. In comparison with the values obtained from pure mode-I fracture modelling (1.78$MPa\sqrt{m}$) and pure mode-II fracture modelling (2.71$MPa\sqrt{m}$), the variations can calculated as follows:

$\frac{K_I-K_{I-Mix}}{K_I}=\frac{1.857-1.81}{1.857}=2.5\%$ (27)

$\frac{K_{II}-K_{II-Mix}}{K_{II}}\frac{2.93-3}{2.93}=-2.3\%$ (28)

The variations 2.5% for mode-I fracture toughness and 2.3% for model II fracture toughness are relatively small if considering the different test methods. Thus the hybrid finite element method can well capture the rock characteristics during the rock fracture process and the hybrid method is a stable and valuable numerical method for rock fracture modelling.

5 NUMERICAL MODELLING OF ROCK FRACTURE UNDER DYNAMIC LOADING

In this section, the rock fracture processes in UCS test and BTS test under various dynamic loading rates, i.e. 1m/s, 5m/s, and 10m/s, are modelled by the hybrid finite-discrete element method. The stress propagation, the fracture initiation and propagation are modelled. The force-loading displacement curves are obtained. The purpose of this section is to demonstrate the capabilities of the hybrid finite-discrete element method in modelling the dynamic rock fracture possess and capturing characteristics of the dynamic behaviour of brittle materials.

5.1 Modelling dynamic rock fracture during UCS test

5.1.1 Modelling rock fracture during UCS test under loading rate of 1m/s

Figure 25a visually illustrates the temporal and spatial stress propagation and the crack initiation and propagation for the UCS test under the loading rate of 1m/s, while the corresponding force-loading displacement curve is plotted in Figure 26 in which the alphabets correspond to those in Figure 25.

Figure 25
Modelled failure process during uniaxial compression test under the loading rate of constant displacement of 1m/s: (a) Evolution of minor principal stress in the vertical direction; (b) Crack initiation and propagation (red: shear failure; Blue: tensile failure)

As for the stress propagation, it is similar to that for the UCS test under quasi-static loading (Figure 7a, under the loading rate of 0.1m/s). As the top plate contacts the specimen, the compressive stress is induced immediately and propagates along the vertical direction of the specimen (Figure 25a-A). As the stress reaches the bottom of the specimen and contacts the bottom plate, the stress is reflected and propagate toward the top plate (Figure 25a-C). Then the stress propagates between the two loading plates a few times, and stress equilibrium is nearly achieved at 37 $\mu s$ (Figure 25a-E). During this period, the stress increases steadily and reaches its peak value (Figure 26-ADEF). Thus the maximum shear stress, i.e. shear strength, has reached. Meanwhile, shear fractures occur and propagate along an inclined angle of approximately 45 degrees to the horizontal direction (Figure 25b-EF). As the top-loading plate continues to move, more cracks are produced and those cracks make the specimen into pieces (Figure 25b-G-L).

As for the force-loading displace curve, the stress drops to the bottom immediately after it reaches its peak, which means the specimen has lost its ability to carry loads (Figure 26-FGL). Overall, the crack initiation and propagation process for UCS test under the loading rate of 1m/s is much similar to that under the loading rate of 0.1 m/s, since the rock property will not be significantly changed before the loading rate reaches a certain threshold level(Zhang 2001).

For the UCS test, under the loading rate of 1m/s, the force-loading displacement still indicates a brittle rock failure process under loading and agree with laboratory results (Wawersik and Fairhurst 1970). It includes a nearly linear loading increase region (Figure 26-AE), a big jump from the peak value to bottom (Figure 26-FI) and a post-failure region (Figure 26-IJ). According to the peak load of the force-loading displacement curve, the compressive strength of the UCS test can be calculated as follows.

Figure 26
Force-loading displacement curve during the uniaxial compression test under a loading rate of 1m/s

${\sigma }_{c1}=\frac{F}{A}=\frac{11.08MN}{0.05m\times 1m}=221.6MPa$ (29)

Where ${\sigma }_{c1}$ is the compressive strength (under the loading rate of 1m/s), and A is the area of the specimen cross-section. The calculated is much higher than the input value or the obtained value under the loading rate of 0.1m/s, which is mainly caused due to the much higher loading rate.

5.1.2 Modelling dynamic rock fracture during UCS test under loading rate of 5m/s, 10m/s and 50m/s

Figure 27 visually illustrates the temporal and spatial rock failure process during UCS test under the loading rate of 5m/s, 10m/s and 50m/s.

For the UCS test under the loading rate of 5m/s, the rock failure process is similar to that under the loading rate of 0.1m/s and 1m/s. A shear crack is firstly produced at the top part of the specimen (Figure 22a-B), then an almost symmetrical shear crack initiates at the bottom part of the specimen (Figure 27a-C). Thus the specimen is manly divided into three pieces by the two almost symmetrical shear cracks (Figure 27a-DE). As the force-loading continues to apply on the top part of the specimen, it forces the three pieces to slide along the formed two shear cracks (Figure 27a-F-I). Finally, the specimen is mainly comprised of three pieces and some small fragments (Figure 27a-J).

For the loading rates of 10m/s and 50m/s, the fractures are mainly produced at the top parts of the specimens(Figure 27b-B and Figure 27c-B) due to that the high-stress concentration at the top part of the specimen. Then the cracks propagate from the top parts of the specimen to the lower part. The propagation processes are similar to the stress propagation as illustrated in Figure 7a. Due to the high loading rate and the high-stress concentration at the top part of the specimen, many fragments are produced at the top vicinity of the specimens (Figure 27b-D-J and Figure 27c-D-J). The final failure pattern consists of fine fragments at the top part of the specimen and relatively bigger fragments at the lower part of the specimen as illustrated in Figure 27b-J and Figure 27c-J.

By comparing the failure processes of UCS test under the loading rates various loading rates, it is found that: Cracks mainly occur at the center of the specimen first for lower loading rates, e.g. 0.1m/s, 1m/s, 5m/s, while the cracks mainly concentrate on the top vicinity of the specimen for higher loading rates, i.e. 10m/s, 50m/s. Besides, the size of the fragments highly depends on the loading rates, i.e. the higher loading rate, the smaller the size of the fragments

Figure 28 demonstrates the force-loading displacement curves for the UCS test under the loading rate of 5m/s, 10m/s and 50m. Those curves also can demonstrate the brittle materials failure process and divided into three main regions, i.e. a nearly linear loading increase region (before the peak of the curve), a big jump from the peak value to the bottom line and post-failure region (the last region)

The peak forces of the obtained force-loading displacement curves can be used to calculate to the dynamic compressive strengths at different loading rates according to the following equations.

${\sigma }_{c5}=\frac{F}{A}=\frac{13.54MN}{0.05m\times 1m}=270.2MPa$ (30)

${\sigma }_{c10}=\frac{F}{A}=\frac{20.58MN}{0.05m\times 1m}=411.6MPa$ (31)

${\sigma }_{c50}=\frac{F}{A}=\frac{61.8MN}{0.05m\times 1m}=1236MPa$ (32)

According to the calculated results of the compressive strengths in the UCS test, the increase in the rate of loading is accompanied by an increase in compressive strength of rock. The modelled results agree well with those documented in literatures that peak load is depended mainly on the rate of loading, as it increased with increasing loading rate (Sukontasukkul, Nimityongskul et al. 2004). In addition, more damage can be found for the specimen suffer higher damage than those subjected to static loading (Sukontasukkul, Nimityongskul et al. 2004).

However, it should be noted that the compressive strengths obtained during the UCS tests under the loading rate of 10m/s and 50m/s might not be accurate. Due to high loading rates, i.e. 10m/s, 50m/s, the stresses are mainly concentrated on loading vicinity and cracks occur before the stress researches an equilibrium state. However, for the Equations 30-32, only when the stress in the specimen reaches equilibrium, the equation is valid.

5.2 Modelling of Dynamic Rock fracture During BTS test

5.2.1 Modelling dynamic rock fracture during BTS test under loading rate of 1m/s

Figure 29 visually illustrates the temporal and spatial rock failure process for BTS test under the loading rate of 1m/s while Figure 30 demonstrates the corresponding force-loading displacement curves. As can be seen from Figure 29a that the stress propagation process is similar to that for the BTS test under quasi-static loading (0.1m/s). Figure 30 depicts the fore-loading displacement curves for the BTS test and the alphabets in the curves correspond to that in Figure 29. The stresses initiate from the top and bottom loading contacts (Figure 29a-A) and propagate along the loading diameter (Figure 29a-B). Then the relative stress equilibrium achieves after the force-loading applied on the loading plates for a while (Figure 29a-C). During this period, on fracture is produced (Figure 29b-AB) and the force-loadings from the two loading plates increase almost linearly (Figure 30a-AB and Figure 30b-AB) and reach their peak points (Figure 30a-C and Figure 30b-C), which mean that the specimen is no longer able to bear more loads. As can be seen that a tensile facture initaites at the centre of the dis and along the loading diameter. Then the fracture propagates towards the loading points (Figure 29b-D). Meanwhile, the force-loadings decrease nearly linearly (Figure 30a-CD and 30b-CD). As the loading plates continue to move, the fracture reaches the loading points and the specimen completely loses its ability to carry loads (Figure 29b-E ). Due to the stress concentration at the two loading points, shear cracks occur as the stress reaches the shear strength of the specimen (Figure 29b-E). Finally, the cracks are mainly along the loading diameter and are comprised of tensile cracks connecting the two loading points and some shear cracks at the loading vicinities (Figure 29b-F).

Figure 27
Modelled failure process during uniaxial compression test under the loading rate of (a) 5m/s, (b) 10m/s and (c) 50m/s, respectively.

According to the maximum load of Figure 30a at point D and Equation 10, the tensile strength of the specimen can be calculated as follows.

${\sigma }_t=\frac{2P}{\pi Dt}=\frac{2\times 1.85\times {10}^6}{3.14\times 0.05\times 1}=23.57MPa$ (33)

Compared with the rock failure process for the BTS test under the loading rate of 0.1m/s, the rock fracture behaviour including the force-loading displacements are quite similar to that under static loading. However, the tensile strength is more than twice larger than the input value, which is induced by the higher loading rate.

5.2.2 Modelling dynamic rock fracture during BTS test under loading rate of 5m/s, 10m/s and 50m/s

To gain insight into the rock behaviour under dynamic loads, three loading rates, i.e. 5m/s, 10m/s, and 50m/s, are applied on the specimens, respectively. Figure 31 visually shows the temporal and spatial distribution of the initiated cracks in the BTS test under different loading rates, while the corresponding force-loading displacement curves are plotted in Figure 32 in which the alphabets correspond to those in Figure 31.

In the case of the rock dynamic failure process under the loading rate of 5m/s, shear cracks, instead of tensile cracks, first occur on the top and bottom ends of the specimen (Figure 31a-B), due to the stress concentrations at the two ends of the disc. This phenomenon is confirmed by literatures (Fairhurst 1964, Hudson, Brown et al. 1972) which indicate that that stress concentration at the vicinity of loading plates occasionally lead early shear failure in the rock. In addition, the relatively hard materials do not permit smooth transition of stress wave from the loaded rim to the unloaded portion of the specimens resulting in a fracture initiation from the perimeter of the disk (Li and Wong 2013).

Figure 28
Force-loading displacement curve during the uniaxial compression test under the loading rates of (a) 5m/s, (b) 10m/s and (c) 50m/s, respectively.

As the two loading plates continue to move, the initiated shear cracks at the vicinity of the loading points keep on propagating, meanwhile, a tensile crack initiates at the centre of the disc along the loading diameter (Figure 31a-C). At the same time, the corresponding force loadings-displacement climb to the summits (Figure 32a-I-C and Figure32 b-I-C), which indicates that the specimen is about to lose its capability of carrying loads. Then the tensile crack at the centre of the loading line propagates towards the two loading points while more shear cracks imitate at the loading vicinities (Figure 31a-DE). As can be seen that on the force-loading curves, the forces drops from the peaks to the bottom (Figure 32a/b-I-CF), which means that the specimen losses its bearing capability. At the same time, the tensile crack at the centre of the loading line reaches the ends of the specimen while more shear cracks coupled with tensile cracks initiating from the loading vicinities propagates towards the opposite ends (Figure 31 a-EF).

While the loading plates further move, the specimen is split into two halves and more shear cracks are produced mainly distributing on the loading vicinities (Figure 31 a-FG).

Therefore, the modelled results agree well with laboratory test (Mahabadi, Cottrell et al. 2010, Zhou, Zou et al. 2013) and numerical simulation results (Zhu and Tang 2006), i.e. a tensile crack initiating at the centre of loading line splitting the disc into two halves, shear cracks appearing at the loading vicinity, more cracks occur during the post-failure period.

The modelled force-loading displacement curve of the BTS recorded in Figure32 also present the typical behavior of brittle rock under compression. According to the maximum load $P_{Max}$ at point C of Figure 32i-a, the tensile strength of the specimen at the loading rate of 5m/s can be calculated as follows.

${\sigma }_{t5}=\frac{2P}{\pi Dt}=\frac{2\times 3.21\times {10}^6}{3.14\times 0.05\times 1}=40.89MPa$ (34)

The modelled tensile strength ($40.89MPa$) is much higher than that under static loading (12.89Mpa), which demonstrates the loading rates significantly influences the strength of the rock specimen.

As can be seen from Figure 31a that the failure process in the BTS test under the loading rate of 5m/s is quite similar to that under the loading rate of 0.1 and 1m/s except that shear cracks occur first rather than tensile cracks and more shear cracks initiate at the loading areas. The force-loading displacement curve for the BTS test under the three loading rates, i.e. 0.1m/s,1m/s and 5m/s, are similar which all represent the typical failure process of brittle materials.

In the case of the rock failure process of BTS test under the loading rates of 10 m/s, the crack initiation and propagation process is much similar to that under the loading rates of 5m/s, i.e. shear cracks initiating first instead of tensile cracks (Figure 31b-B and C) and more shear cracks concentrating on the loading areas (Figure 31b-G). The recorded force-loading displacement cure (Figure 32-II)for the BTS test under the loading rate of 10m/s demonstrates the typical brittle materials failure process again as described for the BTS test under 1m/s and 5m/s. The maximum load $P_{Max}$ from the top loading plate (Figure 32ii-a) can use be to calculate the tensile strength at the loading rate of 10m/s.

${\sigma }_{t10}=\frac{2P}{\pi Dt}=\frac{2\times 4.18\times {10}^6}{3.14\times 0.05\times 1}=53.25MPa$ (35)

For the BTS test under the loading rate of 50m/s, the specimen is dominated by the shear cracks. As can be seen form Fig. 26c-B and C, the cracks initiate from the two loading ends and propagate to the opposite ends. .As the loading plates further move, shear cracks from the two loading area meet and began to interacts (Figure 31c-EF). Finally, with the plates moving towards each other, the loading areas are crushed into fragments and the specimen is dominated by shear cracks (Figure 31c-G).

In the case of force-loading displacement curve (Figure 32 i-c and ii-c), the force goes up and down widely with the force decline gradually, which maybe because of the repeating process of cracks initiating coupled with the closure of cracks and nucleating fragments. The maximum load $P_{Max}$ from the top loading plate at Point B in Figure 22iii-b) is used to calculate the tensile strength at the loading rate of 50m/s

${\sigma }_{t50}=\frac{2P}{\pi Dt}=\frac{2\times 10.67\times {10}^6}{3.14\times 0.05\times 1}=145MPa$ (36)

However, it should be noted that, due to the intensive fluctuation of the force during the BTS test under the loading rate of 50m/s and many shear cracks occurring on the loading areas before tensile cracks initiating, it is hard to know how the force loading contributes to the tensile crack at the center of the loading line. Moreover, the high loading rate, i.e. 50m/s, make it difficult to form a uniform tensile stress distribution along the loading diameter, which plays the role of spilt the disc along the loading line. Therefore, the calculated tensile strength on the loading rate of 50m/s might not be accurate and only can be used for reference.

In summary, with the loading rate increases, the shear crack initiates first rather than the tensile cracks due to the stress concentration at the loading areas. Correspondingly, a higher loading rate contributes to more shear cracks because of the stress concentration on the loading areas. Tensile strength highly depends on the loading rates, i.e. the higher loading rate, the higher tensile strength. Instead of separating the dis into two halves along the loading line, high loading rates, e.g. 50m/s, can crush the specimen into fragments directly. For the BTS test under various loading rates of this study, the force-loading displacement curves represent a typical brittle material failure process as that in quasi-static rock failure processes and this agrees well with those documented in the literatures (Mahabadi, Cottrell et al. 2010, Zhou, Zou et al. 2013). Residual deformation region of the force-loading displacement curve fluctuates significantly if a higher loading rate (i.e. 10m/s) is applied on the loading plates while the entire curve fluctuates intensively if the loading rates is high enough (e.g. 50m/s).

Figure 29
Evaluation of minor principal stress and the fracture initiation and propagation for BTS test under the loading rate of 1m/s: (a) Evolution of minor principal stress; (b) Fracture initiation and propagation

Figure 30
Force-loading displacement curve for during BTS test unde the loading rate of 1m/s :(a) Top force-loading displacement curve; (b) Bottom force-loading displacement curve

Figure 31
Dynamic rock fracture during BTS test under loading rate of 5m/s, 10m/s and 50m/s.

Figure 32
Force-loading displacement curves during BTS test under the loading rate of 5m/s, 10m/s and 50m/s, respectively: (a) Top force-loading displacement curve; (b) Bottom force-loading displacement curve

6 Discussion

6.1 Influence of the loading rate on the rock behavior

In this research, the UCS test, BTS test and NBD test adopt a lower loading rate, i.e. 0.1m/s, to calibrate the proposed method in modelling the rock fracture under quasi-static loading rate. However, it should be noted that the loading rate is much higher than those used in the laboratory (around 0.01mm/s) for the static test. The reason to use a high loading rate is to significantly decrease computational time, as the increase of the loading rate can dramatically decrease the computational time.

According to the hybrid finite-discrete element modelling results of the BTS test under different loading rates, the loading rate significantly influences the rock behaviours. Figure 33 demonstrates the compression of the force-loading displacement curves for the BTS test under various loading rates. Peak force of the curves increases with the loading rate, which indicates the tensile strength increases with the loading rates. Figure 34 illustrate the relationship of the hybrid modelled tensile strength with the various loading rates, i.e. 0.05 m/s、0.1m/s、0.25 m/s、0.5 m/s、0.75 m/s、1m/s、2.5m/s、5 m/s、10 m/s. For the loading rate at a very low level, e.g. 0.05m/s, 0.1m/s, the strength is not significantly influenced by the loading rate, while wit the increase of the loading rate, the rock strength is obviously affected by the loading rate. According to the hybrid finite-discrete element modelling results of the UCS test in section 4.1 and Section 5.1, the compressive strength is influenced by the loading rate in the same way as that for the tensile strength test.

Dynamic strength increasing factor (DIF) is the ratio between the dynamic and static strengths, which is then used here for easy comparison of hybrid finite-discrete element method modelled result with those obtain in literatures. Figure 35 illustrates the comparison of DIF for the hybrid modelling results with those obtained in the literature. As can be seen from Figure 35 that the DIFs obtained both in the literatures and from hybrid modeling increase slightly with the loading rate increases when the loading rate is lower than a threshold level. As the loading rate become larger than the threshold, all the DIFs including modelled results in this research increase dramatically with the increasing of the loading rate. By the comparison, the hybrid finite-discrete element method modelled results under various dynamic loading rates agree well with those obtained in literatures and the hybrid method can well captures the effect of the loading rate on the dynamic properties through the implementation of the equation describing relationship between the static and dynamic loading on the rock strength.

Figure 33
Comparison of the force-loading displacements of the BTS test curves under various loading rates

Figure 34
Relationship of tensile strength and the loading rate

Figure 35
Dynamic strength increasing factors obtained from hybrid modelling and their compassion with the existing results from literatures (Parts of them are taken from Asprone et al. 2009 and Cho et al. 2003).

6.2 effect of the meshes

A typical stress-strain curve can be divided in to stress Hardening part (before peak stress) and stress softening part (after peak stress). The strain hardening part of the stress-strain curve presents no difficulties when implemented in the combined finite-discrete element method, which is therefore implemented in a standard way through the constitutive law(Munjiza 2004). However, the strain softening part of the stress-strain curve is connected with the mesh size and mesh orientation (Munjiza 2004). To deal with this problems, the strain softening is described in the terms of stress and displacement in the hybrid finite-discrete element method. However, as the transition from continuum to discontinuum through fracture and fragmentation is implemented through the separation of the crack elements or joint elements which are placed among the finite elements, the mesh size and mesh orientation still effect the modelling results especially when an inappropriate mesh orientation is adopted. To illustrate the effect of meshes, an UCS test with structural mesh for the model is simulated by the hybrid finite-discrete element method, which is compared with the modelling results with the free mesh.

Figure 36
Numerical model for UCS test with structural mesh

Figure 37
Hybrid finite-discrete element modelling of rock fracture during UCS test with structural mesh for the models

The geometrical model for the UCS test is the same as illustrated in Figure 6a while the numerical model is meshed using structural mesh as show in Figure 36. Figure 37 illustrates the rock fracture process for the UCS under the loading rate of 1m/s and 5m/s modelled by the hybrid finite discrete element method. For the UCS test under the loading rate of 1m/s, the shear crack first initiates from the left-bottom Conner of the specimen (Figure 37a-A). Then it propagates to the right side of the specimen and separates the specimen into two pieces (Figure 37a-B). Further loading causes the two pieces slide along the formed shear cracks (Figure 37a-CD).For the UCS test under the loading rate of 5m/s, a shear crack first appears at the right-top Conner (Figure 37b-B). With the loading continuing, more cracks are produced parallel to the first shear crack (Figure 37a-D). Final fracture pattern (Figure 37a-E) includes both shear and tensile cracks.

For the modelled results using structural mesh (Figure 37a with loading of 1m/s and Figure 37b with loading rate of 5m/s) and free mesh (Figure 25 with loading rate 1m/s and Figure 27a with loading rate of 5m/s), the fracture patterns generally agree with the experimental results as they both have included cracks or/and vertical cracks. However, the fracture patterns for the two kind models are not the same as the adopted different mesh orientations. For the hybrid finite-discrete element model, the fractures occur along the element boundaries. Thus the mesh size and mesh orientation significantly affects the final rock fracture patterns.

6 CONCLUSION

In this study, the hybrid finite-discrete element method is implemented to study the dynamic rock fracture and fragmentation under various loading rate. The hybrid finite discrete element method inherits the advantages of the finite element method in describing elastic deformations and the capabilities of the discrete element method in capturing interactions and fracturing processes of solids. The transition from continuum to discontinuum through fracture and fragmentation makes the hybrid method superior to the traditional continuum-based finite element method and discontinuum-based discrete element method. To model the rock fracture and fragmentation under dynamic loading, the effect of the loading rate is implemented by an empirical relationship between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments. Then the hybrid method is implemented to mode the quasi-static rock fracture in UCS, BTS and NBD test to calibrate the proposed method and demonstrate its ability in modelling the rock fracture process. The hybrid method has modelled the rock fracture process under static loading in BTS test and UCS test. The obtained rock fracture patterns show a good agreement with the experimental results and the modelled force-loading displacement curves indicate a typical failure process of brittle materials. The modelling of the NBD test with the various oriental angles for pre-fabricated notches indicates that the hybrid method not only can well model the rock fracture initiation and propagation but also can capture the fracture characteristics. After that, the hybrid method is used to model the dynamic rock fracture process of UCS and BTS test under various loading rates. The rock fracture and fragmentation processes have obtained. It is found that more fractures are produced compared with that under quasi-static loading and some fragments can be observed at the loading vicinities. The modelled strength is increasing with the loading rate which agrees with those in the documents. The research can conclude that:

• The hybrid finite-discrete element method has well modelled the transition of rock from continuum to discontinuum through fracture and fragmentation in the UCS and BTS tests under quasi-static and dynamic loading.

• The hybrid finite-discrete element method has well modelled the pure Mode-I, pure Mode-II and Mixed-model I-II fracture initiation and propagation in the NBD tests. The successful modelling of the transition from continuum to discontinuum makes the hybrid finite-discrete element method superior to the traditional continuum-based finite element method and discontinuum-based discrete element method.

• The force loading-displacement curves obtained from modelling dynamic rock fracture process in the UCS and BTS tests indicate that the hybrid finite-discrete element can well capture the dynamic fracture behaviour under various dynamic loads through implementing an empirical relation between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments.

• The hybrid finite-discrete element method is a valuable tool for studying the dynamic rock fracture and fragmentation since it can reproduce the fracture initiation and propagation process under various loading rates.

References

• Akazawa, T. (1943). New test method for evaluating internal stress due to compression of concrete (the splitting tension test)(part 1). J Jpn Soc Civ Eng 29: 777-787.
• Al-Shayea, N., K. Khan and S. Abduljauwad (2000). Effects of confining pressure and temperature on mixed-mode (I-II) fracture toughness of a limestone rock. International Journal of Rock Mechanics and Mining Sciences 37:4: 629-643.
• An, H., S. Hou and L. Liu (2019). Experimental and Numerical Study of the Concrete Stress and Fracture Propagation Processes by Blast. Engineering Letters 27:4: 669-675.
• An, H. and H. Liu (2014). Hybrid finite-discrete element modelling of dynamic fracture of rock and resultant fragment arching by rock blast. ISRM International Symposium-8th Asian Rock Mechanics Symposium, International Society for Rock Mechanics and Rock Engineering.
• An, H., H. Liu, H. Han, X. Zheng and X. Wang (2017a). Hybrid finite-discrete element modelling of dynamic fracture and resultant fragment casting and muck-piling by rock blast. Computers and Geotechnics 81: 322-345.
• An, H., H. Liu, X. Wang, J. Shi and H. Han (2017b). HYBRID CONTINUUM-DISCONTINUUM MODELLING OF ROCK FRACUTRE PROCESS IN BRAZILIAN TENSILE STRENGTH TEST. CIVIL ENGINEERING JOURNAL-STAVEBNI OBZOR 3.
• Asprone, D., E. Cadoni and A. Prota (2009). Experimental analysis on tensile dynamic behavior of existing concrete under high strain rates. ACI Structural Journal 106:1: 106.
• Barry, N. W., N. S. Raghu and S. Gexin (1992). Rock fracture mechanics principles design and applications, Amsterdam: ELSEVIER.
• Carneiro, F. (1943). A new method to determine the tensile strength of concrete. Proceedings of the 5th meeting of the Brazilian Association for Technical Rules, 3d. Section.
• Chen, C. S., E. Pan and B. Amadei (1998). Fracture mechanics analysis of cracked discs of anisotropic rock using the boundary element method. International Journal of Rock Mechanics and Mining Sciences 35:2: 195-218.
• Cho, S. H., Y. Ogata and K. Kaneko (2003). Strain-rate dependency of the dynamic tensile strength of rock. International Journal of Rock Mechanics and Mining Sciences 40:5: 763-777.
• Fairhurst, C. (1964). On the validity of the ‘Brazilian’test for brittle materials. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Elsevier.
• Fukuda, D., M. Mohammadnejad, H. Liu, S. Dehkhoda, A. Chan, S. H. Cho, G. J. Min, H. Han, J. i. Kodama and Y. Fujii (2019a). Development of a GPGPU‐parallelized hybrid finite‐discrete element method for modeling rock fracture. International Journal for Numerical and Analytical Methods in Geomechanics 43:10: 1797-1824.
• Fukuda, D., M. Mohammadnejad, H. Liu, Q. Zhang, J. Zhao, S. Dehkhoda, A. Chan, J.-i. Kodama and Y. Fujii (2019b). Development of a 3D hybrid finite-discrete element simulator based on GPGPU-parallelized computation for modelling rock fracturing under quasi-static and dynamic loading conditions. Rock Mechanics and Rock Engineering: 1-34.
• Hobbs, D. (1964). The tensile strength of rocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Elsevier.
• Hondros, G. (1959). The evaluation of Poisson’s ratio and the modulus of materials of a low tensile resistance by the Brazilian (indirect tensile) test with particular reference to concrete. Australian Journal of Applied Science 10:3: 243-268.
• Huaming, A. (2018). Hybrid finite-discrete element model for simulation fracture in rocks under dynamic loading PhD thesis, University of Science and Technology Beijing.
• Hudson, J., E. Brown and F. Rummel (1972). The controlled failure of rock discs and rings loaded in diametral compression. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Elsevier.
• Ingraffea, A. R. (1981). Mixed-mode fracture initiation in Indiana limestone and Westerly granite. The 22nd US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association.
• Jia, Z., A. Castro-Montero and S. P. Shah (1996). Observation of mixed mode fracture with center notched disk specimens. Cement & Concrete Research 26:1: 125-137.
• Khan, K. and N. Al-Shayea (2000). Effect of specimen geometry and testing method on mixed mode I-II fracture toughness of a limestone rock from Saudi Arabia. Rock mechanics and rock engineering 33:3: 179-206.
• Li, D. and L. N. Y. Wong (2013). The Brazilian disc test for rock mechanics applications: review and new insights. Rock mechanics and rock engineering 46:2: 269-287.
• Lisjak, A., O. Mahabadi, L. He, B. Tatone, P. Kaifosh, S. Haque and G. Grasselli (2018). Acceleration of a 2D/3D finite-discrete element code for geomechanical simulations using General Purpose GPU computing. Computers and Geotechnics 100: 84-96.
• Liu, H. and H. An (2016). Hybrid finite-discrete modelling of dynamic fracture of rock under impact loads and resultant fragment grinding. International Conference on Geomechanics, Geo-energy and Geo-resources (IC3G 2016).
• Liu, H., Y. Kang and P. Lin (2015). Hybrid finite-discrete element modeling of geomaterials fracture and fragment muck-piling. International Journal of Geotechnical Engineering 9:2: 115-131.
• Liu, H., S. Kou, P.-A. Lindqvist and C. Tang (2004). Numerical studies on the failure process and associated microseismicity in rock under triaxial compression. Tectonophysics 384:1: 149-174.
• Liu, H., S. Kou, P.-A. Lindqvist and C. Tang (2007). Numerical modelling of the heterogeneous rock fracture process using various test techniques. Rock mechanics and rock engineering 40:2: 107-144.
• Mahabadi, O., B. Cottrell and G. Grasselli (2010). An example of realistic modelling of rock dynamics problems: FEM/DEM simulation of dynamic Brazilian test on Barre granite. Rock mechanics and rock engineering 43:6: 707-716.
• Mahabadi, O., A. Lisjak, A. Munjiza and G. Grasselli (2012). Y-Geo: new combined finite-discrete element numerical code for geomechanical applications. International Journal of Geomechanics 12:6: 676-688.
• Munjiza, A. (2004). The Combined Finite-Discrete Element Method, Wiley Online Library.
• Stefanizzi, S., G. Barla, P. Kaiser and G. Grasselli (2009). Numerical modeling of standard rock mechanics laboratory tests using a finite/discrete element approach. Proceedings of the 3rd CANUS Rock Mechanics Symposium, Toronto.
• Sukontasukkul, P., P. Nimityongskul and S. Mindess (2004). Effect of loading rate on damage of concrete. Cement and Concrete Research 34:11: 2127-2134.
• Timoshenko, S. and J. Goodier (1970). Theory of Elasticity, McGraw-Hill, New York, NY.
• Ulusay, R. (2015). The ISRM suggested methods for rock characterization, testing and monitoring: 2007-2014, Springer.
• Wawersik, W. and C. Fairhurst (1970). A study of brittle rock fracture in laboratory compression experiments. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Elsevier.
• Whittaker, B. N., R. N. Singh and G. Sun (1992). Rock fracture mechanics: principles, design, and applications, Elsevier Amsterdam.
• Zhang, Z. (2001). Laboratory studies of dynamic rock fracture and in-situ measurements of cutter forces for a boring machine.
• Zhang, Z. (2002). An empirical relation between mode I fracture toughness and the tensile strength of rock. International Journal of Rock Mechanics and Mining Sciences 39:3: 401-406.
• Zhang, Z., S. Kou, L. Jiang and P.-A. Lindqvist (2000). Effects of loading rate on rock fracture: fracture characteristics and energy partitioning. International Journal of Rock Mechanics and Mining Sciences 37:5: 745-762.
• Zhang, Z., S. Kou, J. Yu, Y. Yu, L. Jiang and P.-A. Lindqvist (1999). Effects of loading rate on rock fracture. International Journal of Rock Mechanics and Mining Sciences 36:5: 597-611.
• Zhao, J. (2000). Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock. International Journal of Rock Mechanics and Mining Sciences 37:7: 1115-1121.
• Zhou, Z.-l., Y. Zou, X.-b. Li and Y.-h. Jiang (2013). Stress evolution and failure process of Brazilian disc under impact. Journal of Central South University 20: 172-177.
• Zhu, W. and C. Tang (2006). Numerical simulation of Brazilian disk rock failure under static and dynamic loading. International Journal of Rock Mechanics and Mining Sciences 43:2: 236-252.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  18 Sept 2020
• Date of issue
  2020

History

• Received
  19 May 2020
• Reviewed
  23 July 2020
• Accepted
  03 Aug 2020
• Published
  07 Aug 2020