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

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.


Introduction
e 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 fourpoint bending tests. ose 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][3][4][5][6][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.
ere are many continuumbased 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 rocklike 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]. e hybrid methods include the hybrid boundary element method-finite element method (BEM-FEM), discrete element methodfinite 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. ree 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. e 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. e 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.

Fundamental Principles of Hybrid Finite-Discrete Element Method
e 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. us, the users have to learn a few periphery software applications to prepare the model and display the modelled results. e authors [5,7,[15][16][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. e 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 opensource combined finite-discrete element libraries Y2D and Y3D originally developed by Munjiza [12] and Xiang et al. [20], respectively. e 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. e hybrid finite-discrete element method (FDEM) takes advantage of the continuum-based finite element method and discontinuum-based discrete element method. e 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. e 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. e 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. us, the effect of the loading implemented in the hybrid finite-discrete element method is also introduced in this section.

Governing Equation.
e 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. e 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. e generalized governing equation for the motion of the discrete bodies can be expressed as follows [12,21]: 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. e damping diagonal matrix C can be expressed as follows, which accounts for the energy dissipation due to the nonlinear material behaviour [12,21].
where μ is a constant viscous damping coefficient and I is the identity matrix.
2 Shock and Vibration

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. us, 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. ere 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). e total contact force can be obtained by equation (4).
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.
e total contact force can be expressed as a summation over the finite elements.
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 Figure 2. e 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. e 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. is enables the hybrid finite-discrete element to model the transition from continuum to discontinuum through fracture and fragmentation.

Transition from Continuum to
Discontinuum. e 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. e numerical models are meshed using the three-node finite elements, while fournode joint elements are embedded among the boundaries of the finite elements, as illustrated in Figure 3. ree-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.
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. e 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. Figure 4 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: 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 δ. e right part of Figure 4 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: 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:

Shock and Vibration 3
Target body E t Figure 1: Contact force due to penetration of contactor into target (after Munjiza [12]). (a) Contact force due to penetration. (b) e target and contactor discretized into finite elements.   [15]). f p is the peak stress, ε p is the strain at peak stress, and ε u is the ultimate strain.
ree-node constant-strain elastic element where D is a damage variable between 0 and 1 and f(D) 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. e mode-I fracture process is governed by the mode-I fracture energy release rate G fI , 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: In summary, during the model fracture and propagation process, the bonding stress in the normal direction can be expressed using the following equation: where all other parameters have the same meaning as those introduced above. e left part of Figure 4 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: 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. e value can be calculated as follows: After that, bonding stress becomes a pure frictional resistance defined by Columb's model and can be expressed as follows: e shear failure process is governed by pure mode-II fracture energy release rate G fII which can be described by the area under the curve of the bonding stress and the sliding displacement and can be expressed as follows: In summary, the bonding stress for pure mode-II fracture can be calculated using the following equation:

Shock and Vibration
where D is a damage variable between 0 and 1, g(D) 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. e fracture energy release rate G fI− II for the mixed mode I-II equates to the shadow in Figure 5.

Effect of the Loading Rate on the Dynamic Behaviour of
Rock . e effect of the loading rate significantly influences the dynamic rock behaviour. e 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 Section 2.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][24][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: 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).
us, 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.

Calibration of the FDEM
In this section, the FDEM is employed to model the Brazilian tensile strength (BTS) test. e modelled results are then compared with the experiential result to calibrate the proposed method. Figure 6 illustrates the geometrical model and the numerical model for the BTS test. e model is built following the ISRM suggested method [26], i.e., the diameter is 50 m.
e 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. e material properties used for the modelling can be found in Table 1. Figure 7(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)). e force-displacement curves in Figure 8 indicate a typical behaviour of brittle rock under compression. e 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 P Max at point D, the tensile strength of the rock specimen can be calculated as follows:  Shock and Vibration It can be seen that the modelled dynamic tensile strength (34 MPa) is much higher than the input static tensile strength (20 MPa) due to the effects of loading rate.
To better calibrate the hybrid finite-discrete element method, the modelled BTS result (Figure 9(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. us, the modelled results well agree with the typical fracture patterns of BTS tests in the literature. Figure 10 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. e 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. e complete stress distribution along the loading diameter for the BTS test is given by Hondros [30] as follows:  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.

Shock and Vibration
As illustrated in Figure 11, 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  [17]; (c) typical rock failure pattern of BTS after Li and Wong [27] and Hobbs [28].

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. Figure 12(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. e lower two rolls are fixed in both horizontal and vertical directions. During the test, the toploading roll is moving downward at a constant displacement increment. 3PB test is a standard procedure for the determination of the mode-I (tensile) fracture toughness K IC , which can be calculated according to Brown and Srawley [31] as follows: where K IC is the mode-I fracture toughness, P Max 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). Figure 12(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 finitediscrete modelling is simplified as plane strain problems and only the vertical sections are considered.

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. Figure 13 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. e corresponding forceloading displacement curve, force-loading crack mount opening displacement (CMOD) curve, and CMOD time curve are recorded in Figure 15. e alphabets in Figure 15 correspond to those in Figure 14.
As the top-loading roll contacts the beam, the compressive stress at the loading vicinity is produced immediately ( Figure 13 (at 2.5 µs)). en, it propagates downwards ( Figure 13 (at 27.5 µs)). e 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  100 µs)). en, stresses are mainly produced from the three contact points. e 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 Figure 15(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. erefore, there is no opening for a prefabricated crack initially. at is the reason why the force loading-CMOD curve and CMODtime 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)).
e 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 Figures 15(a)  and 15(b), the pure mode-I fracture toughness can be calculated as follows: Liu et al. [32]  e relationship of the mode-I fracture toughness and the energy release rate can be expressed as follows: According to equation (18), the mode-I fracture toughness can be obtained. us, 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. Figure 16 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. e 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.

Under the Loading Rate of 1 m/s and 5 m/s.
(A) in Figures 16(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)). e 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 Shock and Vibration 13 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)). e 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 Figures 17(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)). en, 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. e 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]. e three curves for the 3PB test under the  [33,34]. It should be noted that when force loading increases, the crack initiates before the force reaches its peak. erefore, the hybrid finitediscrete element method shows its ability to model the rock failure process in the 3PB test in terms of crack initiation and propagation. e 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. e 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 Figure 17(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. e fracture toughness for 3PB test under the loading rates of 1 m/s and 5 m/s can be calculated as follows:  Figure 18 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. en, 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. Figures 19(b) and 19(c) depict the prefabricated crack opening process and the corresponding relationship and time, respectively. e 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)). en, 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 (Figure 19(a)), the mode-I fracture toughness under the loading of 10 m/s can be calculated as follows:  Figure 20 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.
(A) in Figure 20 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). en, 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. e 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. Figure 21 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). en, 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). erefore, 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.

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.
e 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. e 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 23 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. is 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.
e 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.

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. us, 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. e time steps can be obtained by the following equation: where l min is the smallest elemental length, λ and μ are Lamé constants, and ρ is the density of geomaterials. us, 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 Figure 24. e top-loading plate moves at the constant displacement of 1 m/s, while the bottom plate is fixed in both horizontal and vertical directions. e 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. us, the mesh orientation significantly affects the rock fracture patterns.

Conclusion
e 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. ree 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. e 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. e FDEM is calibrated by modelling the Brazilian tensile strength test and comparing the modelled results with those well documented in the literature. en, the FDEM is implemented to model the three-point bending test under quasistatic and dynamic loading conditions. e 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) e 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) e 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) e 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.