# FDEM Modelling of Rock Fracture Process during Three-Point Bending Test under Quasistatic and Dynamic Loading Conditions

Authors: Huaming An; Yushan Song; Hongyuan Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

A hybrid finite-discrete element method (FDEM) is proposed to model rock fracture initiation and propagation during a three-point bending test under quasistatic and dynamic loading conditions. Three fracture models have been implemented in the FDEM to model the transition from continuum to discontinuum through fracture and fragmentation. The loading rate effect on rock behaviour has been taken into account by the implementation of the relationship between the static and dynamic rock strengths derived from dynamic rock fracture experiments. The Brazilian tensile strength test has been modelled to calibrate the FDEM. The FDEM can well model the stress and fracture propagation and well show the stress distribution along the vertical diameter of the disc during the Brazilian tensile strength test. Then, FDEM is implemented to study the rock fracture process during three-point bending tests under quasistatic and dynamic loading conditions. The FDEM has well modelled the stress and fracture propagation and can obtain reasonable fracture toughness. After that, the effects of the loading rate on the rock strength and rock fracture toughness are discussed, and the mesh size and mesh orientation on the fracture patterns are also discussed. It is concluded that the FDEM can well model the rock fracture process by the implementation of the three fracture models. The FDEM can capture the loading rate effect on rock strength and rock fracture toughness. The FDEM is a valuable tool for studying the rock behaviour on the dynamic loading although the proposed method is sensitive to the mesh size and mesh orientation.


---

## Body

## 1. Introduction

The rock fracture mechanism plays a significant role in the field of civil and mining engineering as well as other fields, such as geothermal, hydraulic, and oil and gas engineering [1]. It is imperative to study the rock fracture mechanism not only for improving the efficiency of rock breaking and structure demolition but also for preventing the collapse of geo-structures and increasing their stabilities. Many test techniques are employed to study rock fracture mechanisms, e.g., uniaxial compressive strength test, Brazilian tensile strength, notched Brazilian disc test, and three- and four-point bending tests. Those test techniques are used to obtain the general rock properties, e.g., rock strengths and rock toughness. Besides the experiential method, the numerical techniques have provided a better way to understand the rock fracture mechanics [2–7]. Based on the continuous or discontinuous material hypothesis, the numerical techniques include three main types, i.e., continuum-based methods, discontinuum-based methods, and hybrid or combined continuum-discontinuum-based methods. In the case that the discontinuities of rock-like materials can be ignored, the continuum-based method can be well used to solve engineering problems. There are many continuum-based methods available in the literature such as the finite element method (FEM), finite difference method (FDM), boundary element method (BEM), scaled boundary finite element method (SBFEM), and extended finite element method (XFEM) [8].However, in some cases, the discontinuities of the rock-like materials have to be taken into account since the original existing fractures are comparable to the interest area [9]. Under this scenario, the rock mass is assumed to be the assemblies of discrete bodies and the discontinuum-based method can well deal with the interaction of the discrete bodies and can well model the material behaviour under loads. At present, the representative discontinuum-based methods include the distinct element method (DEM), the lattice model (LM) method, molecular dynamics (MD), and discontinuous deformation analysis (DDA). Morris et al. and Mohammadnejad et al. gave a comprehensive review of those methods in computational fracture mechanics of rock. To more realistically model the rock behaviour, many coupled methods, hybrid/combined methods, and multiscale coupled methods have been developed [10]. The hybrid methods include the hybrid boundary element method-finite element method (BEM-FEM), discrete element method-finite element method (DEM-FEM), and discrete element method-boundary element method (DEM-BEM) [11].In this study, a hybrid finite-discrete element method (FDEM) has been proposed to model the rock fracture initiation and propagation process during a three-point bending test under various loading rates. Three fracture models, i.e., pure mode-I, pure mode-II, and mixed mode I-II, are proposed to model the transition from continuum to discontinuum through fracture and fragmentation. A relationship between the static strength and the dynamic strength obtained through experimental tests has been implemented in the FDEM to characterize the loading rate effect on rock behaviours. The loading rate effect is considered through the implementation of a relationship between the static strength and the dynamic strength for modelling the dynamic rock fracture process. The purpose of this study is to illustrate the abilities of the proposed method in modelling the transition from continuum to discontinuum through fracture and fragmentation and to demonstrate the capabilities of the proposed method in capturing the effect of loading rate on rock behaviour.


## 2. Fundamental Principles of Hybrid Finite-Discrete Element Method

The hybrid finite-discrete element method was proposed by Munjiza [12], and he also developed the Y2D library to implement the method. Although the Y2D library is a robust and efficient open source and capable of modelling material continuum/discontinuum behaviour, it is very difficult for other researchers to use due to no graphical user interface. To overcome this shortcoming, Mahabadi et al. [13] developed the Y-GUI, which is a graphical user interface and preprocessor for the hybrid FEM/DEM. Y-GUI can be conveniently used to set up and check models, but it cannot be used for graphically displaying the analyzed results. Munjiza et al. [14] proposed the Virtual Geoscience Workbench (VGW) to simplify the use of Y2D/3D library. For VGW, the GID or GMESH is used to visually prepare the numerical models while MAYAVI is employed to graphically display the calculated results. Thus, the users have to learn a few periphery software applications to prepare the model and display the modelled results. The authors [5, 7, 15–17] proposed Y2D/3D IDE to simply use the Y2D/3D library and visually display the calculated results. In this research, the Y2D/3D IDE is used for modelling the rock failure processes during three-point bending tests under various loading rates. The Y2D/3D IDE is implemented by Liu et al. [15] based on their previous enriched finite element codes RFPA-RT2D [18] and TunGeo3D [19] and the open-source combined finite-discrete element libraries Y2D and Y3D originally developed by Munjiza [12] and Xiang et al. [20], respectively. The Y2D/3D IDE can generate the hybrid FDEM models and set up the initial conditions, physical properties, contact properties, boundary conditions, fracture criteria, and explosive charges. It is also capable of tracing errors visually, displaying the modelled stresses, displacement, velocity, force, damage, fracture, and fragmentation in real-time graphs.The hybrid finite-discrete element method (FDEM) takes advantage of the continuum-based finite element method and discontinuum-based discrete element method. The FDEM not only can model the damage of the brittle materials before fractures occur but also can model the interactions of the fractures and fragments. The transition from continuum to discontinuum is the key issue for the proposed method. Munjiza et al. [12] proposed a combined single and smeared crack model for modelling mode-I fracture of concrete only. The authors extend the model for being able to model three fracture modes, i.e., mode-I, mode-II, and mixed mode, which is introduced in detail. In addition, for modelling the dynamic behaviour of rock-like materials, the loading rate plays a significant role especially when the rock is experiencing strong dynamic loading, e.g., rock blasting. Thus, the effect of the loading implemented in the hybrid finite-discrete element method is also introduced in this section.


### 2.1. Governing Equation

The hybrid finite-discrete element model can have a single discrete body or a number of discrete bodies with general shapes and sizes, each of which is then represented by a single discrete element [15]. Each discrete body is then discretized into three-node triangular finite elements. The central difference explicit time integration scheme is applied in the hybrid finite-discrete method to integrate the equations of motion of either the initially discrete body or the discrete elements formed by the fracture and fragmentation algorithm. The generalized governing equation for the motion of the discrete bodies can be expressed as follows [12, 21]:(1)M∂XX∂t2+C∂X∂t=F,where M and C are the discrete body mass and damping diagonal matrix, respectively, X is the vector of nodal displacements, and F is the node force vector.The damping diagonal matrixC can be expressed as follows, which accounts for the energy dissipation due to the nonlinear material behaviour [12, 21].(2)C=μI,where μ is a constant viscous damping coefficient and I is the identity matrix.


### 2.2. Contact Detection and Interaction

In the FDEM models, thousands or even millions of discrete elements are involved, and it is time-consuming to process contact interaction when there is no contact. Thus, it is essential to detect those couples close to each other and eliminate couples of discrete elements that are far from each other. Contact detection algorithms in the FDEM are employed to detect close couples and eliminate those discrete bodies that are impossible to contact. There are many algorithms available for automatic contact definition and detection in literature [12, 22], such as buffer zone, binary tree, no binary search, space decomposition, and alternating digital tree. After the couples of discrete bodies are detected, the contact interaction algorithms are then employed to enforce the contact constraint between discrete bodies in contact. In the FDEM used in this article, the penalty method is implemented to calculate the contact forces in the tangential and normal directions for the contacted bodies or elements overlapping in space. Figure 1 illustrates two discrete bodies overlapping each other. One is called contactor while the other is called target. For a penetration of area dA of the contactor into a target, the infinitesimal contact force can be given as equation (3). The total contact force can be obtained by equation (4).(3)df=gradφcPc−gradφtPtdA,(4)fc=∫S=Et∩​EcgradφcPc−gradφtPtdA,where df is the infinitesimal contact force due to the infinitesimal overlap dA. φc and φt are potential functions that can be expressed as equations (5) and (6) after the target and contact are discretized into n and m finite element as shown in Figure 2.(5)φt=∑i=1nφti,(6)φc=∑j=1mφcj.Figure 1
Contact force due to penetration of contactor into target (after Munjiza [12]). (a) Contact force due to penetration. (b) The target and contactor discretized into finite elements.
(a)(b)The total contact force can be expressed as a summation over the finite elements.(7)fc=∑i=1n∑j=1m∫Et∩​EcgradφcjPc−gradφtiPtdA.


### 2.3. Constitutive Rock Fracture Model

A constitutive model indicates the stress-strain behaviour of rock mass, which has been incorporated in rock mass behaviour study by numerical modelling. In FDEM modelling, the stress-strain curve is divided into two parts, i.e., strain-hardening part and the strain-softening part, as illustrated in Figure2. The strain-hardening part of the stress-strain curve is implemented through the elastoplastic constitutive law, which is commonly used in continuum methods, such as the finite element method. The strain-softening part of the stress-strain curve indicates the localization of strain, loss of elasticity of the governing equation, and ill-posed problems. To deal with this problem, a fracture model that describes the relationship between stress and displacement is implemented in the hybrid finite-discrete element method. This enables the hybrid finite-discrete element to model the transition from continuum to discontinuum through fracture and fragmentation.Figure 2
Typical stress-strain of geomaterials under loading (after Liu et al. [15]). fp is the peak stress, εp is the strain at peak stress, and εu is the ultimate strain.


### 2.4. Transition from Continuum to Discontinuum

The hybrid finite-discrete element model can have a single discrete body or a number of discrete bodies with general shapes and sizes, each of which is then represented by a single discrete element [15]. Figure 3 illustrates the numerical models under various stress conditions, i.e., tensile, shear, and mixed tensile and shear stress conditions. The numerical models are meshed using the three-node finite elements, while four-node joint elements are embedded among the boundaries of the finite elements, as illustrated in Figure 3. Three-node finite elements are used to analyze the deformation of the numerical models, while the four-node joint elements are employed to model the fracture initiation and propagation, i.e., the transition from continuum to discontinuum.Figure 3
FDEM models under different stress conditions (red line represents edges of joint element). (a) No stress. (b) Under tension condition. (c) Under shear condition. (d) Under both tension and shear conditions.
(a)(b)(c)(d)According to the stress conditions, a fracture could occur in three modes, i.e., pure mode-I, pure mode-II, and mixed mode I-II fractures. The crack initiation and propagation process, i.e., separating the finite elements, is implemented through distortion of the four-node joint elements, which involves bonding stresses. Figure4 illustrates the relationship between the bonding stress and the opening or sliding displacement of adjacent finite elements for pure mode-I or pure mode-II fracture. As illustrated in Figure 3(d), stress at any point of the finite element edge can be divided into the normal direction and tangential direction, which results in the separation of the adjacent finite element surface or the distortion of the joint element in the normal and tangential directions which can be expressed as follows:(8)δ=δnn+δst,where n and t are the unit vectors in the normal and tangential directions, respectively, of the surface at such a point and δn and δs are the magnitudes of the components of δ.Figure 4
Relationship between the bonding stress and opening/sliding displacement under tension and shear conditions (failure criteria for mode-I and mode-II).The right part of Figure4 illustrates the relationship between the bonding stresses and the displacement of the adjacent finite elements in the normal direction for tensile failure, i.e., pure mode-I fracture. As the separation of the finite element surface δn reaches a critical value δnp determined by the tensile strength, σt, the tensile failure occurs. During this period, i.e., 0≤δn≤δnp, the bonding stress in the normal direction can be given by the following equations:(9)σn=2δnδnp−δnδnp2σt.As the separation of the finite element surfaceδn continue increase, i.e., δnp≤δn≤δnu, the bonding stress gradually decreases and can be expressed by the following equation:(10)σn=fDσt,where D is a damage variable between 0 and 1 and fD is the damage function described in the mechanical damage model [18].At the separation,δn≥δnu, the bending stress becomes zero and the fracture is assumed to propagate. The mode-I fracture process is governed by the mode-I fracture energy release rate GfI, which is equal to the area under the curve of the bonding stress and opening displacement as shown on the right side of Figure 4, and can be expressed by the following equation:(11)GfI=∫δnpδnuσnδndδn.In summary, during the model fracture and propagation process, the bonding stress in the normal direction can be expressed using the following equation:(12)σn=2δnδnp−δnδnp2σt,if0≤δn≤δnp,fDσt,ifδnp≤δn≤δnu,0,ifδn≥δnu,where all other parameters have the same meaning as those introduced above.The left part of Figure4 shows the relationship of the bonding stress and sliding separation of the adjacent finite elements for the shear failure process, i.e., pure mode-II fracture. As can be seen in Figure 4, the bonding stress in the tangential direction, i.e., shear stress, increases with the increase of the sliding displacement δs. Before the sliding displacement reaches a critical value δsp described by the shear strength σc, the shear stress can be calculated according to the following equation:(13)τ=2δsδspσc.When the critical valueδsp is reached, the fracture occurs. After that, the bonding stress in the tangential direction decreases with the increase of the sliding displacement till a residual stress δsr according to a mechanical damage model. The value can be calculated as follows:(14)τ=gD.After that, bonding stress becomes a pure frictional resistance defined by Columb’s model and can be expressed as follows:(15)τ=σntan∅f.The shear failure process is governed by pure mode-II fracture energy release rateGfII which can be described by the area under the curve of the bonding stress and the sliding displacement and can be expressed as follows:(16)GfII=∫δspδsrτδs−τrdδs.In summary, the bonding stress for pure mode-II fracture can be calculated using the following equation:(17)τ=2δsδspσc,if0≤δs≤δsp,gD,ifδsp≤δs≤δsr,σntan∅f,ifδs≥δsr,where D is a damage variable between 0 and 1, gD denotes damage functions described in the mechanical damage model [18], and ∅f is the joint residual friction angle.Pure mode-I and pure mode-II fractures are rarely produced. In most cases, a fracture is produced due to the combination of both the shear and tensile stresses, which results in mixed mode I-II fracture. For this situation, if equation (18) is satisfied, the mixed mode I-II fracture occurs. The fracture energy release rate GfI−II for the mixed mode I-II equates to the shadow in Figure 5.(18)δn−δnpδnu−δnp2+δs−δspδsr−δsp2≥1.Figure 5
Failure criteria for mixed I-II fracture mode.


### 2.5. Effect of the Loading Rate on the Dynamic Behaviour of Rock

The effect of the loading rate significantly influences the dynamic rock behaviour. The dynamic rock strength and rock fracture toughness are quite different from those under static loading. In the hybrid finite element method, thefracture propagation process is governed by the fracture energy release rateinstead of rock strength as mentioned in Section2.1. However, the relationship between rock fracture energy release rate and the loading rate is not easy to be obtained directly. Many experiments have been carried out to obtain the influence of the loading rate on rock strength [23–25]. Zhao conducted the uniaxial and triaxial compression, uniaxial tension, and unconfined shear tests on Bukit Timah granite of Singapore and proposed a semilog formula:(19)σcd=Alogσ˙cdσ˙c+σc,where σcd is the dynamic uniaxial compressive strength (MPa), σ˙cd is the dynamic loading rate (MPa/s), σ˙c is the quasistatic loading rate (approximately 5×10−2 MPa/s), σc is the uniaxial compressive strength at the quasistatic loading rate (MPa), and A is a material parameter, which is 11.9 for the Bukit Timah granite [23].In the hybrid finite-discrete element method, the dynamic fracture energy release rate is assumed to increase with the increase of loading rate according to equation (19). Thus, the hybrid finite-discrete element method’s modelled results can reflect the effect of loading rate on dynamic tensile strength, shear strength, and mode-I and mode-II fracture energy release rates.


## 2.1. Governing Equation

The hybrid finite-discrete element model can have a single discrete body or a number of discrete bodies with general shapes and sizes, each of which is then represented by a single discrete element [15]. Each discrete body is then discretized into three-node triangular finite elements. The central difference explicit time integration scheme is applied in the hybrid finite-discrete method to integrate the equations of motion of either the initially discrete body or the discrete elements formed by the fracture and fragmentation algorithm. The generalized governing equation for the motion of the discrete bodies can be expressed as follows [12, 21]:(1)M∂XX∂t2+C∂X∂t=F,where M and C are the discrete body mass and damping diagonal matrix, respectively, X is the vector of nodal displacements, and F is the node force vector.The damping diagonal matrixC can be expressed as follows, which accounts for the energy dissipation due to the nonlinear material behaviour [12, 21].(2)C=μI,where μ is a constant viscous damping coefficient and I is the identity matrix.


## 2.2. Contact Detection and Interaction

In the FDEM models, thousands or even millions of discrete elements are involved, and it is time-consuming to process contact interaction when there is no contact. Thus, it is essential to detect those couples close to each other and eliminate couples of discrete elements that are far from each other. Contact detection algorithms in the FDEM are employed to detect close couples and eliminate those discrete bodies that are impossible to contact. There are many algorithms available for automatic contact definition and detection in literature [12, 22], such as buffer zone, binary tree, no binary search, space decomposition, and alternating digital tree. After the couples of discrete bodies are detected, the contact interaction algorithms are then employed to enforce the contact constraint between discrete bodies in contact. In the FDEM used in this article, the penalty method is implemented to calculate the contact forces in the tangential and normal directions for the contacted bodies or elements overlapping in space. Figure 1 illustrates two discrete bodies overlapping each other. One is called contactor while the other is called target. For a penetration of area dA of the contactor into a target, the infinitesimal contact force can be given as equation (3). The total contact force can be obtained by equation (4).(3)df=gradφcPc−gradφtPtdA,(4)fc=∫S=Et∩​EcgradφcPc−gradφtPtdA,where df is the infinitesimal contact force due to the infinitesimal overlap dA. φc and φt are potential functions that can be expressed as equations (5) and (6) after the target and contact are discretized into n and m finite element as shown in Figure 2.(5)φt=∑i=1nφti,(6)φc=∑j=1mφcj.Figure 1
Contact force due to penetration of contactor into target (after Munjiza [12]). (a) Contact force due to penetration. (b) The target and contactor discretized into finite elements.
(a)(b)The total contact force can be expressed as a summation over the finite elements.(7)fc=∑i=1n∑j=1m∫Et∩​EcgradφcjPc−gradφtiPtdA.


## 2.3. Constitutive Rock Fracture Model

A constitutive model indicates the stress-strain behaviour of rock mass, which has been incorporated in rock mass behaviour study by numerical modelling. In FDEM modelling, the stress-strain curve is divided into two parts, i.e., strain-hardening part and the strain-softening part, as illustrated in Figure2. The strain-hardening part of the stress-strain curve is implemented through the elastoplastic constitutive law, which is commonly used in continuum methods, such as the finite element method. The strain-softening part of the stress-strain curve indicates the localization of strain, loss of elasticity of the governing equation, and ill-posed problems. To deal with this problem, a fracture model that describes the relationship between stress and displacement is implemented in the hybrid finite-discrete element method. This enables the hybrid finite-discrete element to model the transition from continuum to discontinuum through fracture and fragmentation.Figure 2
Typical stress-strain of geomaterials under loading (after Liu et al. [15]). fp is the peak stress, εp is the strain at peak stress, and εu is the ultimate strain.


## 2.4. Transition from Continuum to Discontinuum

The hybrid finite-discrete element model can have a single discrete body or a number of discrete bodies with general shapes and sizes, each of which is then represented by a single discrete element [15]. Figure 3 illustrates the numerical models under various stress conditions, i.e., tensile, shear, and mixed tensile and shear stress conditions. The numerical models are meshed using the three-node finite elements, while four-node joint elements are embedded among the boundaries of the finite elements, as illustrated in Figure 3. Three-node finite elements are used to analyze the deformation of the numerical models, while the four-node joint elements are employed to model the fracture initiation and propagation, i.e., the transition from continuum to discontinuum.Figure 3
FDEM models under different stress conditions (red line represents edges of joint element). (a) No stress. (b) Under tension condition. (c) Under shear condition. (d) Under both tension and shear conditions.
(a)(b)(c)(d)According to the stress conditions, a fracture could occur in three modes, i.e., pure mode-I, pure mode-II, and mixed mode I-II fractures. The crack initiation and propagation process, i.e., separating the finite elements, is implemented through distortion of the four-node joint elements, which involves bonding stresses. Figure4 illustrates the relationship between the bonding stress and the opening or sliding displacement of adjacent finite elements for pure mode-I or pure mode-II fracture. As illustrated in Figure 3(d), stress at any point of the finite element edge can be divided into the normal direction and tangential direction, which results in the separation of the adjacent finite element surface or the distortion of the joint element in the normal and tangential directions which can be expressed as follows:(8)δ=δnn+δst,where n and t are the unit vectors in the normal and tangential directions, respectively, of the surface at such a point and δn and δs are the magnitudes of the components of δ.Figure 4
Relationship between the bonding stress and opening/sliding displacement under tension and shear conditions (failure criteria for mode-I and mode-II).The right part of Figure4 illustrates the relationship between the bonding stresses and the displacement of the adjacent finite elements in the normal direction for tensile failure, i.e., pure mode-I fracture. As the separation of the finite element surface δn reaches a critical value δnp determined by the tensile strength, σt, the tensile failure occurs. During this period, i.e., 0≤δn≤δnp, the bonding stress in the normal direction can be given by the following equations:(9)σn=2δnδnp−δnδnp2σt.As the separation of the finite element surfaceδn continue increase, i.e., δnp≤δn≤δnu, the bonding stress gradually decreases and can be expressed by the following equation:(10)σn=fDσt,where D is a damage variable between 0 and 1 and fD is the damage function described in the mechanical damage model [18].At the separation,δn≥δnu, the bending stress becomes zero and the fracture is assumed to propagate. The mode-I fracture process is governed by the mode-I fracture energy release rate GfI, which is equal to the area under the curve of the bonding stress and opening displacement as shown on the right side of Figure 4, and can be expressed by the following equation:(11)GfI=∫δnpδnuσnδndδn.In summary, during the model fracture and propagation process, the bonding stress in the normal direction can be expressed using the following equation:(12)σn=2δnδnp−δnδnp2σt,if0≤δn≤δnp,fDσt,ifδnp≤δn≤δnu,0,ifδn≥δnu,where all other parameters have the same meaning as those introduced above.The left part of Figure4 shows the relationship of the bonding stress and sliding separation of the adjacent finite elements for the shear failure process, i.e., pure mode-II fracture. As can be seen in Figure 4, the bonding stress in the tangential direction, i.e., shear stress, increases with the increase of the sliding displacement δs. Before the sliding displacement reaches a critical value δsp described by the shear strength σc, the shear stress can be calculated according to the following equation:(13)τ=2δsδspσc.When the critical valueδsp is reached, the fracture occurs. After that, the bonding stress in the tangential direction decreases with the increase of the sliding displacement till a residual stress δsr according to a mechanical damage model. The value can be calculated as follows:(14)τ=gD.After that, bonding stress becomes a pure frictional resistance defined by Columb’s model and can be expressed as follows:(15)τ=σntan∅f.The shear failure process is governed by pure mode-II fracture energy release rateGfII which can be described by the area under the curve of the bonding stress and the sliding displacement and can be expressed as follows:(16)GfII=∫δspδsrτδs−τrdδs.In summary, the bonding stress for pure mode-II fracture can be calculated using the following equation:(17)τ=2δsδspσc,if0≤δs≤δsp,gD,ifδsp≤δs≤δsr,σntan∅f,ifδs≥δsr,where D is a damage variable between 0 and 1, gD denotes damage functions described in the mechanical damage model [18], and ∅f is the joint residual friction angle.Pure mode-I and pure mode-II fractures are rarely produced. In most cases, a fracture is produced due to the combination of both the shear and tensile stresses, which results in mixed mode I-II fracture. For this situation, if equation (18) is satisfied, the mixed mode I-II fracture occurs. The fracture energy release rate GfI−II for the mixed mode I-II equates to the shadow in Figure 5.(18)δn−δnpδnu−δnp2+δs−δspδsr−δsp2≥1.Figure 5
Failure criteria for mixed I-II fracture mode.


## 2.5. Effect of the Loading Rate on the Dynamic Behaviour of Rock

The effect of the loading rate significantly influences the dynamic rock behaviour. The dynamic rock strength and rock fracture toughness are quite different from those under static loading. In the hybrid finite element method, thefracture propagation process is governed by the fracture energy release rateinstead of rock strength as mentioned in Section2.1. However, the relationship between rock fracture energy release rate and the loading rate is not easy to be obtained directly. Many experiments have been carried out to obtain the influence of the loading rate on rock strength [23–25]. Zhao conducted the uniaxial and triaxial compression, uniaxial tension, and unconfined shear tests on Bukit Timah granite of Singapore and proposed a semilog formula:(19)σcd=Alogσ˙cdσ˙c+σc,where σcd is the dynamic uniaxial compressive strength (MPa), σ˙cd is the dynamic loading rate (MPa/s), σ˙c is the quasistatic loading rate (approximately 5×10−2 MPa/s), σc is the uniaxial compressive strength at the quasistatic loading rate (MPa), and A is a material parameter, which is 11.9 for the Bukit Timah granite [23].In the hybrid finite-discrete element method, the dynamic fracture energy release rate is assumed to increase with the increase of loading rate according to equation (19). Thus, the hybrid finite-discrete element method’s modelled results can reflect the effect of loading rate on dynamic tensile strength, shear strength, and mode-I and mode-II fracture energy release rates.


## 3. Calibration of the FDEM

In this section, the FDEM is employed to model the Brazilian tensile strength (BTS) test. The modelled results are then compared with the experiential result to calibrate the proposed method.Figure6 illustrates the geometrical model and the numerical model for the BTS test. The model is built following the ISRM suggested method [26], i.e., the diameter is 50 m. The BTS test is simplified as a plane strain problem, and only the vertical direction is taken into account. As can be seen in Figure 6(a), the rock sample is placed between two loading plates, which will move at a constant displacement increment during the BTS test. Figure 6(b) shows the numerical model for the BTS test. As can be seen in Figure 6(b), the numerical model is discretized by triangle finite elements. The material properties used for the modelling can be found in Table 1.Figure 6
Geometrical and numerical models for the Brazilian tensile strength 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α is the angular distance of load arc, and σxx and σyy are stresses along the horizontal and vertical directions, respectively. (a) Geometrical model. (b) Numerical model.
(a)(b)Table 1
Input parameters of the rock samples for the hybrid finite-discrete element model.
SymbolsPropertiesValuesUnitsEYoung’s modulus60GPavPoisson’s ratio0.26N/AρDensity2600kgm−3σtTensile strength20MPaσcCompressive strength200MPaθInternal friction coefficient30°CuSurface friction coefficient0.1N/AGfIMode-I fracture energy release50Nm−1GfIMode-II fracture energy release250Nm−1Figure7(a) shows the stress propagation, while Figure 7(b) illustrates the rock fracture initiation and propagation. Figure 8 depicts the force-displacement curves during the test. As the top and bottom loading plates contact the specimen, the stresses from the top and bottom of the specimen are produced and propagate towards the center of the disc as illustrated in (A) in Figure 7(a). As the stresses from the two loading areas propagate, a relatively uniform stress distribution along the vertical diameter is formed. Meanwhile, a tensile crack at the center of the specimen along the loading diameter is produced, as demonstrated in (C) in Figure 7(b). During this period, the forces initiated from the loading plates increase almost linearly and reach their peaks ((A)–(D) in Figures 8(a) and 8(b)). As the loading plates continue to move, the crack propagates along the loading diameter and reaches the two contacts between the plates and the specimen ((E) in Figure 7(b)). Due to the strong compressive stress at the contacts, shear cracks are produced at the top and bottom loading vicinities ((F) and (G) in Figure 7(b)). During this period, the forces from the two loading plates drop dramatically as the specimen lost its bearing loads (D–F in Figures 8(a) and 8(b)).Figure 7
FDEM modelling of Brazilian tensile strength test under the loading rate of 1 m/s: (a) evolution of minor principal stress; (b) fracture initiation and propagation; (A)d = 1 μm; (B) d = 20 μm; (C) d = 33.5 μm; (D) d = 39 μm; (E) d = 41 μm; (F) d = 58 μm; (G) d = 166 µm.Figure 8
Force loading-displacement curves obtained during hybrid modelling of the BTS test for a constant displacement increment of 1 m/s on plates: (a) top plate force loading-displacement curve; (b) bottom plate force loading-displacement curve.
(a)(b)The force-displacement curves in Figure8 indicate a typical behaviour of brittle rock under compression. The curves include a compressive deformation region (AB), a linear-elastic deformation region (BC), a nonlinear deformation region (CD), a strain-softening deformation region (DEF), and a residual deformation region (FG). On the basis of the maximum load PMax at point D, the tensile strength of the rock specimen can be calculated as follows:(20)σt1=2PMaxπDt=2×2.86983.14×0.054×1=34MPa.It can be seen that the modelled dynamic tensile strength (34MPa) is much higher than the input static tensile strength (20MPa) due to the effects of loading rate.To better calibrate the hybrid finite-discrete element method, the modelled BTS result (Figure9(a)) is compared with the experimental result (Figure 9(b)) and the typical rock failure pattern in literature (Figure 9(c)) [27, 28]. As can be seen in Figures 9(a) and 9(b), the modelled fracture pattern is quite similar to the experimental result, as they both have a long fracture along the loading diameter and small fragments in the top and bottom loading areas. Figure 9(c) gives the typical fracture pattern of the BTS test under static loading. It has a long fracture along the vertical diameter and two main damage areas on the top and bottom loading vicinities. Thus, the modelled results well agree with the typical fracture patterns of BTS tests in the literature.Figure 9
Comparison of the modelled result with the experimental fracture pattern and typical failure pattern for BTS test: (a) modelled result; (b) experimental result from An et al. [17]; (c) typical rock failure pattern of BTS after Li and Wong [27] and Hobbs [28].
(a)(b)(c)Figure10 compares the fracture propagation process in the BTS test with those well documented in the literature [29]. As shown in Figure 10(a), a typical fracture propagation in the BTS test includes primary tensile cracking along the loading diameter, secondary cracking from the sides parallel to the primary cracks, and tertiary cracking due to shear failure at the contact areas between the loading plates and disc. As can be seen in Figure 10(b), the FDEM well reproduces primary cracks ((A) in Figure 10(b)) and secondary cracks ((B) in Figure 10(b)). In addition, the FDEM also can well model the tertiary cracks ((C) in Figure 10(b)). If the smaller mesh size is adopted in this BTS model, the tertiary cracks might be clearer.Figure 10
Comparison of the crack propagation at different stages between hybrid modelling results and the typical fracture propagation for BTS test. (a) Typical fracture propagation in BTS test [29]. (b) Hybrid modelling results for BTS test at different stages. (A) Primary tensile cracking. (B) Secondary cracking. (C) Tertiary cracking.
(a)(b)The BTS test is developed to obtain the tensile strength since the stress along the loading diameter is almost uniform and mainly distributed in the horizontal direction. The complete stress distribution along the loading diameter for the BTS test is given by Hondros [30] as follows:(21)σx=PπRtα1−r/R2Sin2α1−2r/R2Cos2α+r/R4−tan−11+r/R21−r/R2tanα,σy=−PπRtα1−r/R2Sin2α1−2r/R2Cos2α+r/R4+tan−11+r/R21−r/R2tanα.As illustrated in Figure11, the theoretical and numerical stress distributions along the loading diameter are compared. For better comparison, the stresses are normalized by 2P/πDt for both the numerical and theoretical results. Figure 11 indicates that the stress in the horizontal direction is almost a constant value, except the stress at the top and bottom of the specimen due to stress concentration. For the stress in the vertical direction, it increases from the center of the specimen to the ends of the loading diameter. By comparison, the modelled stress distribution shows a good agreement with the theoretical result. Thus, the hybrid finite-discrete method element can well model the rock failure process in the BTS test and capture the main characters of rock-like material under loading.Figure 11
Comparison of normalized stresses along the vertical diameter between the theoretical (Hondros’ solution) and numerical results.


## 4. FDEM Modelling of Stress and Fracture Propagation during Three-Point Bending Test

In this section, the hybrid finite element method is used to model the rock fracture process during the three-point bending (3PB) test. Figure12(a) depicts the geometrical modes for the symmetrical three-point bending test. As shown in Figure 12(a), a notch is prefabricated at the center of the low part. The lower two rolls are fixed in both horizontal and vertical directions. During the test, the top-loading roll is moving downward at a constant displacement increment.Figure 12
Geometrical and numerical models for 3PB test. (a) Geometrical model. (b) Numerical model.
(a)(b)3PB test is a standard procedure for the determination of the mode-I (tensile) fracture toughnessKIC, which can be calculated according to Brown and Srawley [31] as follows:(22)KIC=PMaxLaBD22.9−4.6aD+21.8aD2−37.6aD3+38.7aD4,where KIC is the mode-I fracture toughness, PMax is the peak load, L is the distance between the supporting points, a is the length of the prefabricated notch, B is the thickness of the rectangular beam, and D is the width. All those parameters are shown in Figure 12(a).Figure12(b) illustrates the numerical model, which is discretized using finite elements. It can be seen from Figure 12(a) that the rectangular beam used for hybrid finite-discrete modelling is simplified as plane strain problems and only the vertical sections are considered.


### 4.1. Under Quasistatic Loading Rate Condition

For modelling the pure mode-I fracture during the 3PB test under quasistatic loading, the constant displacement of 0.1 m/s for the top-loading roll is adopted. Figure13 illustrates the stress propagation during the 3PB test, while Figure 14 depicts the crack initiation and propagation for the 3PBT under the quasistatic loading condition. The corresponding force-loading displacement curve, force-loading crack mount opening displacement (CMOD) curve, and CMOD time curve are recorded in Figure 15. The alphabets in Figure 15 correspond to those in Figure 14.Figure 13
Stress propagation during FDEM modelling three-point bending test under quasistatic loading.Figure 14
FDEM modelling rock fracture process in 3PB test under quasistatic loading.Figure 15
Force-related curves for 3PB test under quasistatic loading condition: (a) force-loading displacement curve; (b) force loading-CMOD curve; (c) CMOD-time curve.
(a)(b)(c)As the top-loading roll contacts the beam, the compressive stress at the loading vicinity is produced immediately (Figure13 (at 2.5 µs)). Then, it propagates downwards (Figure 13 (at 27.5 µs)). The compressive stress can be observed at the loading vicinity due to the stress concentration, while tensile stresses can be found at the tip of the prefabricated notch (Figure 13 (at 27.5 µs)). As the top-loading roll continues to move downwards, stresses from the two bottom contact points are produced (Figure 13 (at 35 µs and 100 µs)). Then, stresses are mainly produced from the three contact points. The tensile stress concentrates at the tip of the prefabricated notch (Figure 13 (from 120 µs to 500 µs)).As the rigid roll on the top of the beam moves downwards, the stress is increasing (A and B in Figure15(a)) while the prefabricated notch starts to open (A and B in Figures 15(b) and 15(c)). It should be noted that the prefabricated cracks in the hybrid finite-discrete model are set by deleting corresponding joint elements. Therefore, there is no opening for a prefabricated crack initially. That is the reason why the force loading-CMOD curve and CMOD-time curve initiate from zero. Before reaching the peak force, no new crack is induced and the increasing rate of the CMOD-time curve is relatively low (A and B in Figure 15(c)), which means the CMOD is gradually increasing although no crack occurs at the tips of the prefabricated notch. A tensile crack first initiates from the tip of the prefabricated crack ((B) in Figure 14) when the force reaches its peak (B in Figure 15(a)). While the rigid roll continues to move downwards, the force drops dramatically (B and C in Figure 15(a)) and the crack continues to propagate upwards ((C) in Figure 14). Meanwhile, although the force drops rapidly, the CMOD continues to increase (B and C in Figures 15(b) and 15(c)) and the increasing rate of CMOD becomes much higher compared with the increasing rate of CMOD before crack occurs. Finally, the crack reaches the top surface of the beam ((D) in Figure 14) while the force drops from the peak force (B in Figure 15(b)) to the bottom (D in Figure 15(b)) and the CMOD achieves its maximum (D in Figure 15(c)).The modelled crack initiation and propagation for the 3PB test agree with the literature that the cracks initiate at the tip of the prefabricated crack and propagate forward to the top central loading supporting points [32]. Additionally, the recorded force loading-CMOD curve agrees with Zhou et al. [33] who did the 3PB test using limestone with range loading rates from 0.0001 mm/s to 0.1 m/s. In their test [33], the prefabricated crack started to open as the loading force increased and a crack firstly was induced at the tip of the prefabricated crack when the loading force increased to its maximum .According to the peak force at point B in Figures15(a) and 15(b), the pure mode-I fracture toughness can be calculated as follows:(23)KIC=PMaxLaBD22.9−4.6aD+21.8aD2−37.6aD3+38.7aD4=0.134×106×3.33D×0.4DBD2×3.13=3.8Pam.Liu et al. [32] modelled the notched Brazilian disc test and the three-point bending test to determine the fracture toughness of the granite using RT2D. The obtained results for mode-I fracture toughness are 1.74 MPa m and 3.74 Pa m for notched Brazilian disc test and three-point bending test, respectively. The obtained mode-I fracture toughness is significantly influenced by the test techniques. The relationship of the mode-I fracture toughness and the energy release rate can be expressed as follows:(24)GfI=KIc2E.According to equation (18), the mode-I fracture toughness can be obtained.(25)KIc=GfI×E=50Nm×60GPa=1.732MPam.The calculated result (1.732) is comparable to modelled result (1.74) by Liu et al. using the Brazilian disc test method, while the hybrid modelled result (3.78) is comparable to RT2D modelled result (3.74) by the same test technique, i.e., 3PB test. Additionally, the hybrid finite-discrete element method is used to model the notched Brazilian disc test to obtain the mode-I fracture toughness for granite in our previous research [5]. The obtained result is 1.857MPam, which is comparable to Liu et al.’s obtained result, i.e., 1.732 MPam.Thus, the hybrid finite-discrete element method not only can model the fracture initiation and propagation process but also can capture the reasonable fracture toughness and reflect the influence of the test techniques.


### 4.2. Under Dynamic Loading Rate Condition

#### 4.2.1. Under the Loading Rate of 1 m/s and 5 m/s

Figure16 visually illustrates the crack initiation and propagation of the 3PB test under the loading rates of 1 m/s and 5 m/s, respectively. It can be seen that the failure processes for the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are much similar. As indicated by Zhang, the rock properties are not significantly changed until a certain threshold of loading rate is achieved. The alphabets in Figure 17 correspond to those in Figure 16 which records the force loading-displacement curves, the force loading-CMOD curves, and the force loading-time curves.Figure 16
Hybrid finite-discrete element method modelled results of 3PBT under the loading rate of (a) 1 m/s and (b) 5 m/s.
(a)(b)Figure 17
(a) Force loading-displacement curve, (b) force loading-CMOD curve, and (c) CMOD-time curve for 3PB test under the loading rates of (i) 1 m/s (figures in the left) and (ii) 5 m/s (figures in the right), respectively.
(a)(b)(c)(A) in Figures16(a) and 16(b) depicts the 3PB test with the prefabricated crack in the bottom center of the beam. As for the 3PB test under the loading rate of 1 m/s, with the rigid roll on the top of the beam moving downwards, the force increases rapidly and peaks at point B (A and B in (i) in Figure 17(a)) which corresponds to the initiation of the crack ((B) in Figure 16(a)).The force increases linearly (A and B in (i) in Figure 17(b)) which corresponds to elastically opening process of the prefabricated crack, and the CMOD has a slight increment (A and B in (i) in Figure 17(c)) before the force reaches its maximum. As the loading roll continues to move downwards, the crack propagates upwards ((C) and (D) in Figure 16(a)) while the force drops dramatically (B–D in (i) in Figure 17(b)) from its maximum. Meanwhile, the CMOD-time curve climbs up considerably (B–D in (i) in Figure 17(c)). Finally, the crack reaches the top surface of the beam and the beam is split into two parts ((E) in Figure 16(a)). Correspondingly, the force drops to its minimum (E in (i) in Figure 17(b)) while the CMOD grows to its maximum (E in (i) in Figure 17(c)).The modelled crack initiation and propagation for the 3PB test agree with the literature that the cracks initiate at the tip of the prefabricated crack and propagate to the top central loading supporting points [32]. Additionally, the recorded force loading-CMOD curve agrees with Zhou et al. [33] who did the 3PB test using limestone with loading rates ranging from 0.0001 mm/s to 0.1 m/s. In their test [33], the prefabricated crack was opening with the force increasing and a crack firstly was induced from the tip of the prefabricated crack when the force reached its maximum. Figure 13 indicates the same conclusion.In terms of the rock failure process of the 3PB test under the loading rate of 5 m/s, the crack initiation and propagation are almost the same as those under the loading rates of 1 m/s and 0.1 m/s. Initially, while the loading roll contacts the beam, the force is induced and almost increases linearly ((A) and (B) in (ii) in Figures17(a) and 17(b)). Correspondingly, the prefabricated crack is elastically opened while the force loading is increasing ((A) and (B) in (ii) in Figure 17(a). However, the CMOD opens slowly before a newly formed crack is produced ((A) and (B) in (ii) in Figures 17(b) and 17(c)). Then, a crack initiates at the tip of the prefabricated crack ((B) in Figure  16(b)) when the force reaches point B (B in (ii) in Figure 17(a)) rather than at the maximum force (C in (ii) in Figure 17(a)). While the loading roll further moves downwards, the force achieves its peak (C in (ii) in Figures 17(a) and 17(b)) and the crack continues to propagate forwards ((C) in Figure 16(b)). It can be seen from (ii) in Figure 17(c) that the CMOD rises considerably after passing point B which is considered as a turning point for expressing the material properties objectively [34]. After force achieves its maximum, it began to drop considerably (C and D in (ii) in Figures 17(a) and 17(b)) and the crack further propagates upwards ((D) in Figure 16(b)). Finally, while the crack reaches the top surface, the force drops to the bottom and the CMOD achieves its maximum.The model rock failure processes in the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are similar in terms of crack initiation and propagation and agree with the literature [32]. The three curves for the 3PB test under the loading rates of 1 m/s and 5 m/s show good agreements with those well documented in the literature [33, 34]. It should be noted that when force loading increases, the crack initiates before the force reaches its peak. Therefore, the hybrid finite-discrete element method shows its ability to model the rock failure process in the 3PB test in terms of crack initiation and propagation. The modelled failure processes under the three loading rates, i.e., 0.1 m/s, 1 m/s, and 5 m/s, show good agreements with those well documented in the literature. The corresponding curves, i.e., force-loading displacement curve, force-loading CMOD curve, and the CMOD time curve, also agree well with those in literature.As when the forces reach Point B (B in (i) in Figure17(a), B in (ii) in Figure 17(a)), the cracks start to open at the tips of the prefabricated notches for the 3PB test under the loading rate of 1 m/s and 5 m/s, respectively; the forces at Point B (B in (i) in Figure 17(a), B in (ii) in Figure 17(a)) can be used to calculate the fracture toughness for 3PB test under the loading rate of 1 m/s and 5 m/s, respectively. The fracture toughness for 3PB test under the loading rates of 1 m/s and 5 m/s can be calculated as follows:(26)KIC1=0.697×106×3.33×D0.4×D1×D2×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.70.44=19.76MPam,KIC5=4.563×106×3.33×540.4×54×10−31×542×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.7×0.44=129.32MPam.


### 4.3. Under the Loading Rate of 10 m/s

Figure18 visually depicts the fracture initiation and propagation process of the 3PB test under the loading rate of 10 m/s while the alphabets correspond to those in the force loading-displacement curve, the force loading-CMOD curve, and the CMOD-time curve. As the loading rigid roll moves downwards, it can be seen in Figure 18 that a tensile crack at the tip of the prefabricated notch and a shear crack at the vicinity of the loading area occur at the same time instead of only a tensile crack initiated at the tip of the preexisting crack, while the force loading increases dramatically and peaks at point B as shown in Figure 19(a). Due to the much higher loading rate (10 m/s), strong compressive stress concentrates at the loading vicinity which causes shear cracksat the loading area. Then, the force drops rapidly which indicates that the beam loses its bearing compatibility. During the postfailure process (B–D in Figure 19(a)), the tensile crack initiated from the tip of the prefabricated notch continues to propagate upwards while more shear cracks and tensile cracks are produced at the vicinity of the loading area, which propagates downwards. During the residual deformation process (D and H in Figure 19(a)), more shear cracks and tensile cracks are produced and some tensile cracks initiate from the tips of the shear cracks and propagate towards the lower rigid rolls as shown in (E) in Figure 18 (the two long red cracks). Meanwhile, some tensile cracks initiate from the tips of the shear cracks and propagate towards the top surface of the beam as shown in (E) and (F) in Figure 18 (the red cracks at the loading vicinity). Finally, the tensile cracks initiated from the tip of the prefabricated notch connect with the cracks induced around the loading vicinity while the loading area is crushed into fragments due to the strong stress concentration as shown in (H) in Figure 18.Figure 18
Crack initiation and propagation of 3PB test under the constant displacement increment of 10 m/s.Figure 19
Force loading-displacement curve for 3PB test constant displacement of 10 m/s. (a) Force loading-displacement curve. (b) Force loading-CMOD curve. (c) CMOD-time curve.
(a)(b)(c)Figures19(b) and 19(c) depict the prefabricated crack opening process and the corresponding relationship and time, respectively. The prefabricated notch is elastically opened before the crack initiates at the tip of the prefabricated crack (A and B in Figure 19(b)) while the CMOD increases slowly (A and B in Figure 19(c)). Then, the force peaks at point B when a crack is induced at the tip of the prefabricated crack. After that, the force drops dramatically (B, C, and D in Figure 19(b)) and finally, the force fluctuates at a very low range (E–H in Figure 19(b)). However, for the CMOD, it increased rapidly after passing point B (B in Figure 19(c)) and finally it reaches its maximum (H in Figure 19(c)).For the 3PB test under the loading rate of 10 m/s, although many shear cracks occur at the loading vicinity due to the stress concentration, the mode-I rock failure process agrees with the literature [34] in terms of crack initiation and propagation, the relationship of crack opening and the loading, and the crack propagation speed.According to the peak force on the force-loading displacement curve at point B (Figure19(a)), the mode-I fracture toughness under the loading of 10 m/s can be calculated as follows:(27)KIC10=9.5×106×3.33×540.4×54×10−31×542×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.7×0.44=269.30MPm.


### 4.4. Under the Loading Rate of 50 m/s

Figure20 depicts the fracture initiation and propagation of the 3PB test under the constant increment displacement of 50 m/s, while Figure 21 recorded the force loading-displacement curve in which the alphabets correspond to those in Figure 20.Figure 20
Crack initiation and propagation of 3PB test under the constant displacement increment of 50 m/s.Figure 21
Force loading-displacement curve for 3PB test constant displacement of 50 m/s.(A) in Figure20 shows the initial beam with the prefabricated notch at the low center of the beam. As the rigid roll contacts the beam, shear cracks initiate immediately due to the stress concentration at the contact area ((B) in Figure 20). Then, the induced cracks propagate down while more cracks are produced including both shear cracks and tensile cracks ((C) in Figure 20). While more produced cracks propagate towards the lower surface of the beam, a crack initiates from the tips of the prefabricated notch ((D) in Figure 20). It seems that the induced crack from the tip of the prefabricated notch plays little role in the rock failure process of the beam. The induced cracks due to the strong stress concentration reach the prefabricated notch and continue to propagate downwards ((E) in Figure 20). With the rigid roll moving downwards, more tensile cracks are produced, which are mainly distributed at the left and right sides of the beam ((F) and (G) in Figure 20). Finally, the loading area is crushed into fragments while more cracks are produced which are distributed along the loading area.Figure21 demonstrates that the force increases rapidly (A and B in Figure 21) as the rigid roll contacts the beam due to the high loading rate (50 m/s). Then, the force drops gradually (B–F in Figure 21) and finally the beam completely loses its ability to carry loads (G and H in Figure 21). Therefore, although the beam is under dynamic loads, it demonstrates the typical brittle material failure process. Since the prefabricated notch does not play a significant role in the rock failure process under the loading rate of 50 m/s, the peak force on the force-loading displacement curve is not used to calculate the fracture toughness under the dynamic load, i.e., 50 m/s.


## 4.1. Under Quasistatic Loading Rate Condition

For modelling the pure mode-I fracture during the 3PB test under quasistatic loading, the constant displacement of 0.1 m/s for the top-loading roll is adopted. Figure13 illustrates the stress propagation during the 3PB test, while Figure 14 depicts the crack initiation and propagation for the 3PBT under the quasistatic loading condition. The corresponding force-loading displacement curve, force-loading crack mount opening displacement (CMOD) curve, and CMOD time curve are recorded in Figure 15. The alphabets in Figure 15 correspond to those in Figure 14.Figure 13
Stress propagation during FDEM modelling three-point bending test under quasistatic loading.Figure 14
FDEM modelling rock fracture process in 3PB test under quasistatic loading.Figure 15
Force-related curves for 3PB test under quasistatic loading condition: (a) force-loading displacement curve; (b) force loading-CMOD curve; (c) CMOD-time curve.
(a)(b)(c)As the top-loading roll contacts the beam, the compressive stress at the loading vicinity is produced immediately (Figure13 (at 2.5 µs)). Then, it propagates downwards (Figure 13 (at 27.5 µs)). The compressive stress can be observed at the loading vicinity due to the stress concentration, while tensile stresses can be found at the tip of the prefabricated notch (Figure 13 (at 27.5 µs)). As the top-loading roll continues to move downwards, stresses from the two bottom contact points are produced (Figure 13 (at 35 µs and 100 µs)). Then, stresses are mainly produced from the three contact points. The tensile stress concentrates at the tip of the prefabricated notch (Figure 13 (from 120 µs to 500 µs)).As the rigid roll on the top of the beam moves downwards, the stress is increasing (A and B in Figure15(a)) while the prefabricated notch starts to open (A and B in Figures 15(b) and 15(c)). It should be noted that the prefabricated cracks in the hybrid finite-discrete model are set by deleting corresponding joint elements. Therefore, there is no opening for a prefabricated crack initially. That is the reason why the force loading-CMOD curve and CMOD-time curve initiate from zero. Before reaching the peak force, no new crack is induced and the increasing rate of the CMOD-time curve is relatively low (A and B in Figure 15(c)), which means the CMOD is gradually increasing although no crack occurs at the tips of the prefabricated notch. A tensile crack first initiates from the tip of the prefabricated crack ((B) in Figure 14) when the force reaches its peak (B in Figure 15(a)). While the rigid roll continues to move downwards, the force drops dramatically (B and C in Figure 15(a)) and the crack continues to propagate upwards ((C) in Figure 14). Meanwhile, although the force drops rapidly, the CMOD continues to increase (B and C in Figures 15(b) and 15(c)) and the increasing rate of CMOD becomes much higher compared with the increasing rate of CMOD before crack occurs. Finally, the crack reaches the top surface of the beam ((D) in Figure 14) while the force drops from the peak force (B in Figure 15(b)) to the bottom (D in Figure 15(b)) and the CMOD achieves its maximum (D in Figure 15(c)).The modelled crack initiation and propagation for the 3PB test agree with the literature that the cracks initiate at the tip of the prefabricated crack and propagate forward to the top central loading supporting points [32]. Additionally, the recorded force loading-CMOD curve agrees with Zhou et al. [33] who did the 3PB test using limestone with range loading rates from 0.0001 mm/s to 0.1 m/s. In their test [33], the prefabricated crack started to open as the loading force increased and a crack firstly was induced at the tip of the prefabricated crack when the loading force increased to its maximum .According to the peak force at point B in Figures15(a) and 15(b), the pure mode-I fracture toughness can be calculated as follows:(23)KIC=PMaxLaBD22.9−4.6aD+21.8aD2−37.6aD3+38.7aD4=0.134×106×3.33D×0.4DBD2×3.13=3.8Pam.Liu et al. [32] modelled the notched Brazilian disc test and the three-point bending test to determine the fracture toughness of the granite using RT2D. The obtained results for mode-I fracture toughness are 1.74 MPa m and 3.74 Pa m for notched Brazilian disc test and three-point bending test, respectively. The obtained mode-I fracture toughness is significantly influenced by the test techniques. The relationship of the mode-I fracture toughness and the energy release rate can be expressed as follows:(24)GfI=KIc2E.According to equation (18), the mode-I fracture toughness can be obtained.(25)KIc=GfI×E=50Nm×60GPa=1.732MPam.The calculated result (1.732) is comparable to modelled result (1.74) by Liu et al. using the Brazilian disc test method, while the hybrid modelled result (3.78) is comparable to RT2D modelled result (3.74) by the same test technique, i.e., 3PB test. Additionally, the hybrid finite-discrete element method is used to model the notched Brazilian disc test to obtain the mode-I fracture toughness for granite in our previous research [5]. The obtained result is 1.857MPam, which is comparable to Liu et al.’s obtained result, i.e., 1.732 MPam.Thus, the hybrid finite-discrete element method not only can model the fracture initiation and propagation process but also can capture the reasonable fracture toughness and reflect the influence of the test techniques.


## 4.2. Under Dynamic Loading Rate Condition

### 4.2.1. Under the Loading Rate of 1 m/s and 5 m/s

Figure16 visually illustrates the crack initiation and propagation of the 3PB test under the loading rates of 1 m/s and 5 m/s, respectively. It can be seen that the failure processes for the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are much similar. As indicated by Zhang, the rock properties are not significantly changed until a certain threshold of loading rate is achieved. The alphabets in Figure 17 correspond to those in Figure 16 which records the force loading-displacement curves, the force loading-CMOD curves, and the force loading-time curves.Figure 16
Hybrid finite-discrete element method modelled results of 3PBT under the loading rate of (a) 1 m/s and (b) 5 m/s.
(a)(b)Figure 17
(a) Force loading-displacement curve, (b) force loading-CMOD curve, and (c) CMOD-time curve for 3PB test under the loading rates of (i) 1 m/s (figures in the left) and (ii) 5 m/s (figures in the right), respectively.
(a)(b)(c)(A) in Figures16(a) and 16(b) depicts the 3PB test with the prefabricated crack in the bottom center of the beam. As for the 3PB test under the loading rate of 1 m/s, with the rigid roll on the top of the beam moving downwards, the force increases rapidly and peaks at point B (A and B in (i) in Figure 17(a)) which corresponds to the initiation of the crack ((B) in Figure 16(a)).The force increases linearly (A and B in (i) in Figure 17(b)) which corresponds to elastically opening process of the prefabricated crack, and the CMOD has a slight increment (A and B in (i) in Figure 17(c)) before the force reaches its maximum. As the loading roll continues to move downwards, the crack propagates upwards ((C) and (D) in Figure 16(a)) while the force drops dramatically (B–D in (i) in Figure 17(b)) from its maximum. Meanwhile, the CMOD-time curve climbs up considerably (B–D in (i) in Figure 17(c)). Finally, the crack reaches the top surface of the beam and the beam is split into two parts ((E) in Figure 16(a)). Correspondingly, the force drops to its minimum (E in (i) in Figure 17(b)) while the CMOD grows to its maximum (E in (i) in Figure 17(c)).The modelled crack initiation and propagation for the 3PB test agree with the literature that the cracks initiate at the tip of the prefabricated crack and propagate to the top central loading supporting points [32]. Additionally, the recorded force loading-CMOD curve agrees with Zhou et al. [33] who did the 3PB test using limestone with loading rates ranging from 0.0001 mm/s to 0.1 m/s. In their test [33], the prefabricated crack was opening with the force increasing and a crack firstly was induced from the tip of the prefabricated crack when the force reached its maximum. Figure 13 indicates the same conclusion.In terms of the rock failure process of the 3PB test under the loading rate of 5 m/s, the crack initiation and propagation are almost the same as those under the loading rates of 1 m/s and 0.1 m/s. Initially, while the loading roll contacts the beam, the force is induced and almost increases linearly ((A) and (B) in (ii) in Figures17(a) and 17(b)). Correspondingly, the prefabricated crack is elastically opened while the force loading is increasing ((A) and (B) in (ii) in Figure 17(a). However, the CMOD opens slowly before a newly formed crack is produced ((A) and (B) in (ii) in Figures 17(b) and 17(c)). Then, a crack initiates at the tip of the prefabricated crack ((B) in Figure  16(b)) when the force reaches point B (B in (ii) in Figure 17(a)) rather than at the maximum force (C in (ii) in Figure 17(a)). While the loading roll further moves downwards, the force achieves its peak (C in (ii) in Figures 17(a) and 17(b)) and the crack continues to propagate forwards ((C) in Figure 16(b)). It can be seen from (ii) in Figure 17(c) that the CMOD rises considerably after passing point B which is considered as a turning point for expressing the material properties objectively [34]. After force achieves its maximum, it began to drop considerably (C and D in (ii) in Figures 17(a) and 17(b)) and the crack further propagates upwards ((D) in Figure 16(b)). Finally, while the crack reaches the top surface, the force drops to the bottom and the CMOD achieves its maximum.The model rock failure processes in the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are similar in terms of crack initiation and propagation and agree with the literature [32]. The three curves for the 3PB test under the loading rates of 1 m/s and 5 m/s show good agreements with those well documented in the literature [33, 34]. It should be noted that when force loading increases, the crack initiates before the force reaches its peak. Therefore, the hybrid finite-discrete element method shows its ability to model the rock failure process in the 3PB test in terms of crack initiation and propagation. The modelled failure processes under the three loading rates, i.e., 0.1 m/s, 1 m/s, and 5 m/s, show good agreements with those well documented in the literature. The corresponding curves, i.e., force-loading displacement curve, force-loading CMOD curve, and the CMOD time curve, also agree well with those in literature.As when the forces reach Point B (B in (i) in Figure17(a), B in (ii) in Figure 17(a)), the cracks start to open at the tips of the prefabricated notches for the 3PB test under the loading rate of 1 m/s and 5 m/s, respectively; the forces at Point B (B in (i) in Figure 17(a), B in (ii) in Figure 17(a)) can be used to calculate the fracture toughness for 3PB test under the loading rate of 1 m/s and 5 m/s, respectively. The fracture toughness for 3PB test under the loading rates of 1 m/s and 5 m/s can be calculated as follows:(26)KIC1=0.697×106×3.33×D0.4×D1×D2×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.70.44=19.76MPam,KIC5=4.563×106×3.33×540.4×54×10−31×542×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.7×0.44=129.32MPam.


## 4.2.1. Under the Loading Rate of 1 m/s and 5 m/s

Figure16 visually illustrates the crack initiation and propagation of the 3PB test under the loading rates of 1 m/s and 5 m/s, respectively. It can be seen that the failure processes for the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are much similar. As indicated by Zhang, the rock properties are not significantly changed until a certain threshold of loading rate is achieved. The alphabets in Figure 17 correspond to those in Figure 16 which records the force loading-displacement curves, the force loading-CMOD curves, and the force loading-time curves.Figure 16
Hybrid finite-discrete element method modelled results of 3PBT under the loading rate of (a) 1 m/s and (b) 5 m/s.
(a)(b)Figure 17
(a) Force loading-displacement curve, (b) force loading-CMOD curve, and (c) CMOD-time curve for 3PB test under the loading rates of (i) 1 m/s (figures in the left) and (ii) 5 m/s (figures in the right), respectively.
(a)(b)(c)(A) in Figures16(a) and 16(b) depicts the 3PB test with the prefabricated crack in the bottom center of the beam. As for the 3PB test under the loading rate of 1 m/s, with the rigid roll on the top of the beam moving downwards, the force increases rapidly and peaks at point B (A and B in (i) in Figure 17(a)) which corresponds to the initiation of the crack ((B) in Figure 16(a)).The force increases linearly (A and B in (i) in Figure 17(b)) which corresponds to elastically opening process of the prefabricated crack, and the CMOD has a slight increment (A and B in (i) in Figure 17(c)) before the force reaches its maximum. As the loading roll continues to move downwards, the crack propagates upwards ((C) and (D) in Figure 16(a)) while the force drops dramatically (B–D in (i) in Figure 17(b)) from its maximum. Meanwhile, the CMOD-time curve climbs up considerably (B–D in (i) in Figure 17(c)). Finally, the crack reaches the top surface of the beam and the beam is split into two parts ((E) in Figure 16(a)). Correspondingly, the force drops to its minimum (E in (i) in Figure 17(b)) while the CMOD grows to its maximum (E in (i) in Figure 17(c)).The modelled crack initiation and propagation for the 3PB test agree with the literature that the cracks initiate at the tip of the prefabricated crack and propagate to the top central loading supporting points [32]. Additionally, the recorded force loading-CMOD curve agrees with Zhou et al. [33] who did the 3PB test using limestone with loading rates ranging from 0.0001 mm/s to 0.1 m/s. In their test [33], the prefabricated crack was opening with the force increasing and a crack firstly was induced from the tip of the prefabricated crack when the force reached its maximum. Figure 13 indicates the same conclusion.In terms of the rock failure process of the 3PB test under the loading rate of 5 m/s, the crack initiation and propagation are almost the same as those under the loading rates of 1 m/s and 0.1 m/s. Initially, while the loading roll contacts the beam, the force is induced and almost increases linearly ((A) and (B) in (ii) in Figures17(a) and 17(b)). Correspondingly, the prefabricated crack is elastically opened while the force loading is increasing ((A) and (B) in (ii) in Figure 17(a). However, the CMOD opens slowly before a newly formed crack is produced ((A) and (B) in (ii) in Figures 17(b) and 17(c)). Then, a crack initiates at the tip of the prefabricated crack ((B) in Figure  16(b)) when the force reaches point B (B in (ii) in Figure 17(a)) rather than at the maximum force (C in (ii) in Figure 17(a)). While the loading roll further moves downwards, the force achieves its peak (C in (ii) in Figures 17(a) and 17(b)) and the crack continues to propagate forwards ((C) in Figure 16(b)). It can be seen from (ii) in Figure 17(c) that the CMOD rises considerably after passing point B which is considered as a turning point for expressing the material properties objectively [34]. After force achieves its maximum, it began to drop considerably (C and D in (ii) in Figures 17(a) and 17(b)) and the crack further propagates upwards ((D) in Figure 16(b)). Finally, while the crack reaches the top surface, the force drops to the bottom and the CMOD achieves its maximum.The model rock failure processes in the 3PB test under the loading rates of 0.1 m/s, 1 m/s, and 5 m/s are similar in terms of crack initiation and propagation and agree with the literature [32]. The three curves for the 3PB test under the loading rates of 1 m/s and 5 m/s show good agreements with those well documented in the literature [33, 34]. It should be noted that when force loading increases, the crack initiates before the force reaches its peak. Therefore, the hybrid finite-discrete element method shows its ability to model the rock failure process in the 3PB test in terms of crack initiation and propagation. The modelled failure processes under the three loading rates, i.e., 0.1 m/s, 1 m/s, and 5 m/s, show good agreements with those well documented in the literature. The corresponding curves, i.e., force-loading displacement curve, force-loading CMOD curve, and the CMOD time curve, also agree well with those in literature.As when the forces reach Point B (B in (i) in Figure17(a), B in (ii) in Figure 17(a)), the cracks start to open at the tips of the prefabricated notches for the 3PB test under the loading rate of 1 m/s and 5 m/s, respectively; the forces at Point B (B in (i) in Figure 17(a), B in (ii) in Figure 17(a)) can be used to calculate the fracture toughness for 3PB test under the loading rate of 1 m/s and 5 m/s, respectively. The fracture toughness for 3PB test under the loading rates of 1 m/s and 5 m/s can be calculated as follows:(26)KIC1=0.697×106×3.33×D0.4×D1×D2×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.70.44=19.76MPam,KIC5=4.563×106×3.33×540.4×54×10−31×542×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.7×0.44=129.32MPam.


## 4.3. Under the Loading Rate of 10 m/s

Figure18 visually depicts the fracture initiation and propagation process of the 3PB test under the loading rate of 10 m/s while the alphabets correspond to those in the force loading-displacement curve, the force loading-CMOD curve, and the CMOD-time curve. As the loading rigid roll moves downwards, it can be seen in Figure 18 that a tensile crack at the tip of the prefabricated notch and a shear crack at the vicinity of the loading area occur at the same time instead of only a tensile crack initiated at the tip of the preexisting crack, while the force loading increases dramatically and peaks at point B as shown in Figure 19(a). Due to the much higher loading rate (10 m/s), strong compressive stress concentrates at the loading vicinity which causes shear cracksat the loading area. Then, the force drops rapidly which indicates that the beam loses its bearing compatibility. During the postfailure process (B–D in Figure 19(a)), the tensile crack initiated from the tip of the prefabricated notch continues to propagate upwards while more shear cracks and tensile cracks are produced at the vicinity of the loading area, which propagates downwards. During the residual deformation process (D and H in Figure 19(a)), more shear cracks and tensile cracks are produced and some tensile cracks initiate from the tips of the shear cracks and propagate towards the lower rigid rolls as shown in (E) in Figure 18 (the two long red cracks). Meanwhile, some tensile cracks initiate from the tips of the shear cracks and propagate towards the top surface of the beam as shown in (E) and (F) in Figure 18 (the red cracks at the loading vicinity). Finally, the tensile cracks initiated from the tip of the prefabricated notch connect with the cracks induced around the loading vicinity while the loading area is crushed into fragments due to the strong stress concentration as shown in (H) in Figure 18.Figure 18
Crack initiation and propagation of 3PB test under the constant displacement increment of 10 m/s.Figure 19
Force loading-displacement curve for 3PB test constant displacement of 10 m/s. (a) Force loading-displacement curve. (b) Force loading-CMOD curve. (c) CMOD-time curve.
(a)(b)(c)Figures19(b) and 19(c) depict the prefabricated crack opening process and the corresponding relationship and time, respectively. The prefabricated notch is elastically opened before the crack initiates at the tip of the prefabricated crack (A and B in Figure 19(b)) while the CMOD increases slowly (A and B in Figure 19(c)). Then, the force peaks at point B when a crack is induced at the tip of the prefabricated crack. After that, the force drops dramatically (B, C, and D in Figure 19(b)) and finally, the force fluctuates at a very low range (E–H in Figure 19(b)). However, for the CMOD, it increased rapidly after passing point B (B in Figure 19(c)) and finally it reaches its maximum (H in Figure 19(c)).For the 3PB test under the loading rate of 10 m/s, although many shear cracks occur at the loading vicinity due to the stress concentration, the mode-I rock failure process agrees with the literature [34] in terms of crack initiation and propagation, the relationship of crack opening and the loading, and the crack propagation speed.According to the peak force on the force-loading displacement curve at point B (Figure19(a)), the mode-I fracture toughness under the loading of 10 m/s can be calculated as follows:(27)KIC10=9.5×106×3.33×540.4×54×10−31×542×2.9−4.6×0.4+21.8×0.42−37.6×0.43+38.7×0.44=269.30MPm.


## 4.4. Under the Loading Rate of 50 m/s

Figure20 depicts the fracture initiation and propagation of the 3PB test under the constant increment displacement of 50 m/s, while Figure 21 recorded the force loading-displacement curve in which the alphabets correspond to those in Figure 20.Figure 20
Crack initiation and propagation of 3PB test under the constant displacement increment of 50 m/s.Figure 21
Force loading-displacement curve for 3PB test constant displacement of 50 m/s.(A) in Figure20 shows the initial beam with the prefabricated notch at the low center of the beam. As the rigid roll contacts the beam, shear cracks initiate immediately due to the stress concentration at the contact area ((B) in Figure 20). Then, the induced cracks propagate down while more cracks are produced including both shear cracks and tensile cracks ((C) in Figure 20). While more produced cracks propagate towards the lower surface of the beam, a crack initiates from the tips of the prefabricated notch ((D) in Figure 20). It seems that the induced crack from the tip of the prefabricated notch plays little role in the rock failure process of the beam. The induced cracks due to the strong stress concentration reach the prefabricated notch and continue to propagate downwards ((E) in Figure 20). With the rigid roll moving downwards, more tensile cracks are produced, which are mainly distributed at the left and right sides of the beam ((F) and (G) in Figure 20). Finally, the loading area is crushed into fragments while more cracks are produced which are distributed along the loading area.Figure21 demonstrates that the force increases rapidly (A and B in Figure 21) as the rigid roll contacts the beam due to the high loading rate (50 m/s). Then, the force drops gradually (B–F in Figure 21) and finally the beam completely loses its ability to carry loads (G and H in Figure 21). Therefore, although the beam is under dynamic loads, it demonstrates the typical brittle material failure process. Since the prefabricated notch does not play a significant role in the rock failure process under the loading rate of 50 m/s, the peak force on the force-loading displacement curve is not used to calculate the fracture toughness under the dynamic load, i.e., 50 m/s.


## 5. Discussion

### 5.1. Influence of the Loading Rate on Rock Behaviour

In our previous research, the loading rate on the rock strength has been studied [5]. In this section, the effects of the loading rate on the rock strength and rock fracture toughness are studied. Figure 22 depicts the influence of the loading rate on the tensile strength. The vertical axis indicates the tensile strength while the horizontal axis shows the logarithm value of the strain rate. As illustrated in Figure 22, the tensile strength increases slightly with the strain rate when the loading rate is lower than a threshold level, i.e., logarithm value equates to 1. After that, the tensile strength increases dramatically with the strain rate. The FDEM modelled results show a good agreement with those well documented in the literature [35, 36]. In those research studies, there is a threshold separating the loading rate into static and dynamic conditions. In the static loading condition, the strength increases slowly with the increase of loading rate, while in the dynamic condition, the strength increases dramatically with the increase of loading rate.Figure 22
Effect of loading rate on the rock strength (the data are taken from Huaming et al. [5]).Figure23 illustrates the effect of loading rate on rock fracture toughness. It can be seen that the fracture toughness remains almost constant before 1 for the logarithm value of the strain rate. After that, it increases dramatically with the increase of the strain rate. This shows the same character for the loading effect on rock strength as shown in Figure 22. Additionally, the effect of the loading rate on fracture toughness shows a good agreement with Zhang et al. [25], who did a series of fracture toughness experimental tests by means of wedge loading applied to short-rock fracture specimens. The experimental results show that fracture toughness nearly remains constant before the loading rate reaches a certain value, and after the certain value, it increases with the loading rate.Figure 23
Loading effect on the rock fracture toughness.


### 5.2. Influence of the Mesh Size and Mesh Orientation

For modelling the rock fracture using the hybrid finite-discrete element method, the modelling results are sensitive to the mesh size and mesh orientation. As mentioned before, the strain-hardening part of the typical stress-strain curve is implemented in the hybrid method through the constitutive laws [12], while the strain-softening part is implemented through the separation of the crack elements or joint elements which are placed among the finite elements. Thus, the fractures can only occur along the mesh boundaries and the mesh size and orientation can significantly influence the modelling results. In addition, the central difference explicit integration scheme is implemented in the hybrid method to obtain the solutions for forces, stresses, and displacements. The time steps can be obtained by the following equation:(28)tcr≈lminλ+2μ/ρ,where lmin is the smallest elemental length, λ and μ are Lamé constants, and ρ is the density of geomaterials.Thus, the mesh size can greatly influence the time steps and can further influence the computational time. For saving the computational resources, relatively bigger mesh size is adopted in this research.To illustrate the influence of mesh orientation on the rock behaviour, the free mesh and the structural mesh are adopted for the uniaxial compressive test model as shown in Figure24. The top-loading plate moves at the constant displacement of 1 m/s, while the bottom plate is fixed in both horizontal and vertical directions. The modelled results for the two types of meshed models generally agree with the experimental results as they both have included tensile cracks. However, the fracture propagation paths in the twomodelled results are both limited by the mesh orientations as the fracturesonly can occur along the element boundaries. Thus, the mesh orientation significantly affects the rock fracture patterns.Figure 24
Influence of the mesh orientation on the failure patterns. (a) Free mesh model and the corresponding failure pattern. (b) Structural mesh model and the corresponding failure pattern.
(a)(b)


## 5.1. Influence of the Loading Rate on Rock Behaviour

In our previous research, the loading rate on the rock strength has been studied [5]. In this section, the effects of the loading rate on the rock strength and rock fracture toughness are studied. Figure 22 depicts the influence of the loading rate on the tensile strength. The vertical axis indicates the tensile strength while the horizontal axis shows the logarithm value of the strain rate. As illustrated in Figure 22, the tensile strength increases slightly with the strain rate when the loading rate is lower than a threshold level, i.e., logarithm value equates to 1. After that, the tensile strength increases dramatically with the strain rate. The FDEM modelled results show a good agreement with those well documented in the literature [35, 36]. In those research studies, there is a threshold separating the loading rate into static and dynamic conditions. In the static loading condition, the strength increases slowly with the increase of loading rate, while in the dynamic condition, the strength increases dramatically with the increase of loading rate.Figure 22
Effect of loading rate on the rock strength (the data are taken from Huaming et al. [5]).Figure23 illustrates the effect of loading rate on rock fracture toughness. It can be seen that the fracture toughness remains almost constant before 1 for the logarithm value of the strain rate. After that, it increases dramatically with the increase of the strain rate. This shows the same character for the loading effect on rock strength as shown in Figure 22. Additionally, the effect of the loading rate on fracture toughness shows a good agreement with Zhang et al. [25], who did a series of fracture toughness experimental tests by means of wedge loading applied to short-rock fracture specimens. The experimental results show that fracture toughness nearly remains constant before the loading rate reaches a certain value, and after the certain value, it increases with the loading rate.Figure 23
Loading effect on the rock fracture toughness.


## 5.2. Influence of the Mesh Size and Mesh Orientation

For modelling the rock fracture using the hybrid finite-discrete element method, the modelling results are sensitive to the mesh size and mesh orientation. As mentioned before, the strain-hardening part of the typical stress-strain curve is implemented in the hybrid method through the constitutive laws [12], while the strain-softening part is implemented through the separation of the crack elements or joint elements which are placed among the finite elements. Thus, the fractures can only occur along the mesh boundaries and the mesh size and orientation can significantly influence the modelling results. In addition, the central difference explicit integration scheme is implemented in the hybrid method to obtain the solutions for forces, stresses, and displacements. The time steps can be obtained by the following equation:(28)tcr≈lminλ+2μ/ρ,where lmin is the smallest elemental length, λ and μ are Lamé constants, and ρ is the density of geomaterials.Thus, the mesh size can greatly influence the time steps and can further influence the computational time. For saving the computational resources, relatively bigger mesh size is adopted in this research.To illustrate the influence of mesh orientation on the rock behaviour, the free mesh and the structural mesh are adopted for the uniaxial compressive test model as shown in Figure24. The top-loading plate moves at the constant displacement of 1 m/s, while the bottom plate is fixed in both horizontal and vertical directions. The modelled results for the two types of meshed models generally agree with the experimental results as they both have included tensile cracks. However, the fracture propagation paths in the twomodelled results are both limited by the mesh orientations as the fracturesonly can occur along the element boundaries. Thus, the mesh orientation significantly affects the rock fracture patterns.Figure 24
Influence of the mesh orientation on the failure patterns. (a) Free mesh model and the corresponding failure pattern. (b) Structural mesh model and the corresponding failure pattern.
(a)(b)


## 6. Conclusion

The hybrid finite-discrete element method (FDEM) is proposed to model the rock fracture initiation and propagation process during a three-point bending test under quasistatic and dynamic loading conditions. Three fracture modes are implemented in the FDEM to model the transition from continuum to discontinuum through fracture and fragmentation, which make the FDEM superior to the traditional continuum-based finite element method and discontinuum-based discrete element method. The FDEM takes advantage of finite element method in describing elastic deformations and the capabilities of the discrete element method in capturing interactions and fracturing processes of solids. Additionally, the loading rate effect on the rock behaviour has been taken into account by the implementation of an empirical relationship between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments. The FDEM is calibrated by modelling the Brazilian tensile strength test and comparing the modelled results with those well documented in the literature. Then, the FDEM is implemented to model the three-point bending test under quasistatic and dynamic loading conditions. The FDEM has well modelled the stress and fracture propagation and can well capture the loading effect on the rock behaviours. Finally, the effect of loading rate and mesh size and mesh orientation are discussed. It is concluded that(1)
The FDEM has well modelled the rock fracture process during the three-point bending test, and this indicates that the FDEM can well model the transition of rock from continuum to discontinuum through fracture and fragmentation by implementation of three fracture modes.(2)
The FDEM has well captured the effect of the loading rate on the rock strength, rock fracture toughness, and rock fracture behaviour by implementing an empirical relationship between the static strengths and the dynamic strengths derived from the dynamic rock fracture experiments.(3)
The FDEM is a valuable tool to study the rock behaviour since the FDEM takes advantage 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.



---

*Source: 5566992-2021-08-05.xml*