Experimental and Numerical Analysis of Mode I Fracture Process of Rock by Semi-Circular Bend Specimen

The semi-circular bend (SCB) specimen is widely used to measure fracture toughness of brittle materials such as rock. In this work, the stress field, fracture process zone (FPZ), and crack propagation velocity of SCB specimen are analyzed during the fracture process of rock specimens. The FPZ of specimen is obtained by experimental and numerical methods under a three-point bend test. The stress concentration zones of specimen present a heart shape at peak load points. FPZ forms before macro fracture occurs. The macro fracture form inside FPZ in a post-peak region of a load–displacement curve. The crack propagation process of specimen include two stages, namely the rapid crack initial development stage, and the final crack splitting stage. The maximum crack propagation velocity of specimen is about 267 m/s, and the average crack propagation velocity is about 111 m/s.


Introduction
Fracture toughness is defined to describe the capacity of materials to resist crack propagation, which include three categories under different loading statuses [1]. Crack propagation in brittle materials such as rock and concrete can be easily produced by mode I loading because tensile strength of these materials is lower than compressive strength. In addition, Mode I fracture toughness of rock is proportional to tensile strength under quasistatic and low-rate impact loads [2]. Therefore, most studies of rock fracture toughness focus on mode I [3][4][5][6][7][8][9][10]. Since the semi-circular bend (SCB) specimen is easy to fabricate and is convenient to test, it has been suggested to determine mode I fracture toughness by the International Society for Rock Mechanics [11]. It is widely used in testing fracture toughness of brittle materials such as rock. Kuruppu et al. [12] believed that the SCB specimen can be well used to measure dynamic fracture toughness of rocks. Funatsu et al. [13] obtained effects of temperature and confining pressure on mixed-mode (I-II) toughness of SCB specimens. Gao et al. [14] calculated the fracture initiation toughness and fracture propagation toughness of SCB specimens by using the digital image correlation method. Wei et al. [15] predicted SCB specimen's fracture using the maximum tangential strain criterion.
However, it is not easy to fabricate sharp enough notch tips for the SCB specimen made of fine-grained hard rocks, and the notch tip is often an arc which may lead to the change of crack initiation position. In the numerical model, there is an angle at the notch tip to ensure that the crack initiation position is in the center of tip because the model size can be accurately determined by coordinates. Xu et al. [16] calculated fracture initiation toughness and fracture energy of SCB specimens by the discrete element method. Dai et al. [17] regarded that the fracture of the SCB specimen agreed with the measuring principle by numerical investigation. Wei et al. [18] observed the fracture process zone of SCB specimens by both acoustic emission monitoring and numerical simulation. On The loading diagram of SCB granite specimen is shown in Figure 1. A 5 mm long notch is machined at the bottom of the specimen with a cutting machine and the radius of notch tip is 0.5 mm. The corresponding geometric parameters of specimen are listed in Table 2. Test is carried out on MTS Landmark electro-hydraulic servo control material testing machine at Central South University. Load velocity is 1 mm/min. Loading diagram of the SCB specimen (Note: R is radius of specimen; B is thickness; a is notch length; s is distance between two supporting rollers; P is applied load). π Specimen Lighting Figure 1. Loading diagram of the SCB specimen (Note: R is radius of specimen; B is thickness; a is notch length; s is distance between two supporting rollers; P is applied load). In addition, the DIC method is applied to record the fracture process of granite specimens. White paint is sprayed on the surface of specimen, and then black speckle is brushed. Speckle diameter on the specimen surface is 0.33 mm (0.013 inches). Two CMOS cameras (GS3-U3-123S6M-C) are used to capture the failure process of granite specimens. The photography frequency is set to 10 fps and takes a picture every 0.1 s. The speckle displacement change is analyzed by a software (Vic-2D 6.0) to determine the evolution of the strain field of the specimen surface. The two cameras are arranged on the left and right sides, and the photos taken by the camera may make it feel that the loading pin is not in the notch line. It can be seen from the last figure in Figure 1 that the loading pin is in notch line.
The calculation formula and processing of test data are based on the suggested method by ISRM. Fracture toughness K IC of SCB specimen is calculated by Equation (1) [11]: where P max is the peak load of specimen failure; and Y is dimensionless stress intensity factor. Table 3 lists test results of three SCB specimens, where the peak displacement is displacement corresponding to peak load. The fracture toughness is 1.627, 1.586, and 1.597 MPa · m 1/2 , respectively, indicating that fracture toughness of the granite is relatively homogeneous. It provides a basis for subsequent numerical simulation.  Figure 2 shows the load-displacement curve of three specimens. The load-displacement curve of three specimens are slightly different due to the difference in geometrical dimensions, but their general shape and variation trend are same. It can be seen that the slope of curve is small during initial loading, and increases significantly with the increase of displacement (after about 0.03 mm). When reaching the peak load point, specimens suddenly fail and load capacity suddenly drops to zero. It is difficult to get a post-peak curve of specimens due to the sudden failure of granite, which is supplemented in the next numerical results.

Experimental Results
The post processing of the DIC procedure is as follows: firstly, a study area is selected on the surface of specimen; then, the region is divided into several subsets, and the initial displacements of these subsets are calculated; the corresponding subset is searched based on assumption, and the displacement is calculated; finally, the deformation field of specimen is obtained by processing. The calculation configurations of VIC-2D 6.0 software in DIC post-processing are as follows [32]: subset: 15, step: 4, subset weights: Gaussian weights, error estimation: about 26 µε, correlation criterion: normalized squared differences. ment curve of three specimens are slightly different due to the dimensions, but their general shape and variation trend are sa slope of curve is small during initial loading, and increases sign of displacement (after about 0.03 mm). When reaching the pe suddenly fail and load capacity suddenly drops to zero. It is d curve of specimens due to the sudden failure of granite, whi next numerical results. The post processing of the DIC procedure is as follows: firs on the surface of specimen; then, the region is divided into seve displacements of these subsets are calculated; the correspondin on assumption, and the displacement is calculated; finally, the imen is obtained by processing. The calculation configurations DIC post-processing are as follows [32]: subset: 15, step: 4, weights, error estimation: about 26 με, correlation criterion: n ences. Figure 3 shows the maximum principal strain evolution p during the test. It can be seen from Figure 3b,c that the tensile the tip of the notch (indicated by red in the figure), and then te rapidly upward, resulting in instantaneous failure of the speci about 0.5% (0.005), this area can be regarded as damaged, so Figure 3c can be considered as the fracture process zone. Figure mode of the specimen; therefore, there is no colorbar. Figure 3e, of the specimen.  Figure 3 shows the maximum principal strain evolution process of specimen SCB-1 during the test. It can be seen from Figure 3b,c that the tensile strain zone first occurs at the tip of the notch (indicated by red in the figure), and then tensile strain zone develops rapidly upward, resulting in instantaneous failure of the specimen. When strain reaches about 0.5% (0.005), this area can be regarded as damaged, so the blue-enclosed part in Figure 3c can be considered as the fracture process zone. Figure 3d shows the final failure mode of the specimen; therefore, there is no colorbar. Figure 3e,f show the fracture surface of the specimen.
Due to the rapid failure process of specimen, the static CMOS camera captures a few pictures (10 fps) during the crack propagation process of specimens. The failure process and crack propagation speed are very difficult to obtain accurately, and stress field of the specimen cannot be obtained by the present experiment. Therefore, the FDEM numerical simulation is used in this study to reveal more details during the failure process. The FDEM numerical method has been well applied in rock mechanics. Mahabadi et al. [19] used a coupled finite-discrete element method to simulate behavior of Brazilian disk specimens under dynamic loading. Li et al. [33] studied mechanical properties and the failure process of marble rings under static and dynamic loads by finite-discrete element. The FDEM numerical method can well simulate mechanical properties and the failure process of rocks under external load.  Due to the rapid failure process of specimen, the static CMOS camera captures a few pictures (10 fps) during the crack propagation process of specimens. The failure process and crack propagation speed are very difficult to obtain accurately, and stress field of the specimen cannot be obtained by the present experiment. Therefore, the FDEM numerical simulation is used in this study to reveal more details during the failure process. The FDEM numerical method has been well applied in rock mechanics. Mahabadi et al. [19] used a coupled finite-discrete element method to simulate behavior of Brazilian disk specimens under dynamic loading. Li et al. [33] studied mechanical properties and the failure process of marble rings under static and dynamic loads by finite-discrete element. The FDEM numerical method can well simulate mechanical properties and the failure process of rocks under external load.

Numerical Model and Results
The FDEM numerical modeling is performed by software ELFEN. The modeling geometric size is set according to specimen SCB-1 as shown in Figure 4. The radius, thickness, and notch length are 23.37 mm, 25.28 mm, and 5.85 mm, respectively. The model is subdivided into triangular meshes with a side length of 1 mm, and the end of prefabricated notch is sharpened. The international society of rock mechanics and rock engineer-

Numerical Model and Results
The FDEM numerical modeling is performed by software ELFEN. The modeling geometric size is set according to specimen SCB-1 as shown in Figure 4. The radius, thickness, and notch length are 23.37 mm, 25.28 mm, and 5.85 mm, respectively. The model is subdivided into triangular meshes with a side length of 1 mm, and the end of prefabricated notch is sharpened. The international society of rock mechanics and rock engineering recommends that the specimen diameter be more than 10 times the material particle size [11]. The specimen diameter in the numerical model is about 50 mm, larger than 10 times of the mesh size (10 mm), which basically meets the requirements. In addition, the influence of the mesh size on the loading curves is discussed in Appendix A.
The material parameters in model are the same as the granite's. The material properties and discrete contact parameters are listed in Table 4. The Mohr-Coulomb criterion with tension cut-off is used to judge the material's damage as shown in Figure 5, which can judge both tensile and shear fracture of rock. The fracture energy of hard and brittle materials is usually between 0.01 and 0.3 N/mm [19]. Based on previous research, 0.08 N/mm is suitable to describe the hard and brittle granite in this work.
with tension cut-off is used to judge the material's damage as shown in Figure 5, which can judge both tensile and shear fracture of rock. The fracture energy of hard and brittle materials is usually between 0.01 and 0.3 N/mm [19]. Based on previous research, 0.08 N/mm is suitable to describe the hard and brittle granite in this work.
Displacement control is applied on top loading platen. Constant velocity is 0.3 mm/s and strain rate is about 0.012 s −1 . Bottom platen is a fixed in the model. Upper platen is a non-fixed with an applied constant velocity.      Figure 6 shows a load-displacement curve of specimens in laboratory and numerical models. The two curves agree well with each other. It indicates that the numerical model can simulate the mechanical properties of granite specimens quite well. Therefore, the numerical simulation is reliable, and it can provide us with more details about the mechanical behavior of specimens. Displacement control is applied on top loading platen. Constant velocity is 0.3 mm/s and strain rate is about 0.012 s −1 . Bottom platen is a fixed in the model. Upper platen is a non-fixed with an applied constant velocity. Figure 6 shows a load-displacement curve of specimens in laboratory and numerical models. The two curves agree well with each other. It indicates that the numerical model can simulate the mechanical properties of granite specimens quite well. Therefore, the numerical simulation is reliable, and it can provide us with more details about the mechanical behavior of specimens. Figure 6 shows a load-displacement curve of specimens in laboratory and numer models. The two curves agree well with each other. It indicates that the numerical mo can simulate the mechanical properties of granite specimens quite well. Therefore, numerical simulation is reliable, and it can provide us with more details about the chanical behavior of specimens. The simulation is performed quasi-static. The crack is produced when the ten stress exceeds the tensile strength, and the strain energy reaches the energy required failure. The length of crack propagation is equal to the strain energy released of the no divided by the failure energy per unit area of the material. As shown in Figure 7, in crack is generated from the tip of notch of the specimen, and then grows zigzag tow the top. Eventually, the specimen split into two halves along the middle. Numerical s ulation results reproduce failure process of the SCB specimen. The simulation is performed quasi-static. The crack is produced when the tensile stress exceeds the tensile strength, and the strain energy reaches the energy required for failure. The length of crack propagation is equal to the strain energy released of the notch divided by the failure energy per unit area of the material. As shown in Figure 7, initial crack is generated from the tip of notch of the specimen, and then grows zigzag toward the top. Eventually, the specimen split into two halves along the middle. Numerical simulation results reproduce failure process of the SCB specimen.

Stress Field Evolution
During loading process, the stress field on specimen surface can be used to judge th Macro-cracks Figure 7. The crack evolution process of numerical simulation and actual failure mode.

Stress Field Evolution
During loading process, the stress field on specimen surface can be used to judge the local stress concentration and the overall stress distribution, which is an important information to judge whether the specimen is damaged or not. It is difficult to get the stress field on the surface of specimens in laboratory tests, but the information of the stress field on the surface of specimens can be well obtained in the numerical calculation software.
As shown in Figures 6 and 8, the peak load is 4357 N with the maximum vertical displacement of 0.097 mm under numerical simulation. The specimen is in the tensile state, except the part in contact with platen. The tensile stress concentration is initially created at the tip of pre-notch as shown in Figure 8a. The tensile stress concentration zone of the specimen presents a heart shape at the peak stress point as shown in Figure 8c. After the crack is generated (shown in Figure 8d-f), the maximum tensile stress concentration zone moves up with the crack and expands upward until the crack penetrates the specimen.

Stress Field Evolution
During loading process, the stress field on specimen surface can be used to judge the local stress concentration and the overall stress distribution, which is an important information to judge whether the specimen is damaged or not. It is difficult to get the stress field on the surface of specimens in laboratory tests, but the information of the stress field on the surface of specimens can be well obtained in the numerical calculation software.
As shown in Figures 6 and 8, the peak load is 4357 N with the maximum vertical displacement of 0.097 mm under numerical simulation. The specimen is in the tensile state, except the part in contact with platen. The tensile stress concentration is initially created at the tip of pre-notch as shown in Figure 8a. The tensile stress concentration zone of the specimen presents a heart shape at the peak stress point as shown in Figure 8c. After the crack is generated (shown in Figure 8d-f), the maximum tensile stress concentration zone moves up with the crack and expands upward until the crack penetrates the specimen. For Mode I dominated problems, the failure of rock is dominated by the formation of tensile cracks (Mode I) when the tensile stress exceeds the tensile strength of the material. It can be seen from Figure 9 that the shear stress is concentrated on the top of the specimen and does not cause cracking. The cracks and branches are caused by tensile stress; they are all Mode I crack. For Mode I dominated problems, the failure of rock is dominated by the formation of tensile cracks (Mode I) when the tensile stress exceeds the tensile strength of the material. It can be seen from Figure 9 that the shear stress is concentrated on the top of the specimen and does not cause cracking. The cracks and branches are caused by tensile stress; they are all Mode I crack. For Mode I dominated problems, the failure of rock is dominated by the formation of tensile cracks (Mode I) when the tensile stress exceeds the tensile strength of the material. It can be seen from Figure 9 that the shear stress is concentrated on the top of the specimen and does not cause cracking. The cracks and branches are caused by tensile stress; they are all Mode I crack. Figure 9. The shear stress at the peak load and the unit is MPa.

Evolution of the Fracture Process Zone
The result difference in rock fracture toughness test obtained by different specimen types has attracted the attention of many scholars. Most scholars believe that the fracture process zone (FPZ) near the crack tip is responsible for variation of fracture toughness [35][36][37][38][39]. Therefore, it is necessary to analyze distribution and evolution of the specimen's FPZ during the loading process, in order to provide more understanding of rock fractures.
The damaging behavior and plastic strain law in the numerical model are shown in Figure 10. During the loading process, the increase of stress causes the elastic strain of material. When the stress of a point reaches the tensile strength , the elastic strain reaches the maximum . The strain after the peak elastic strain is plastic strain . The plastic strain will continue to increase when the energy stored is not enough to produce cracks. The is the strain when the material is fractured. The area under the Figure 9. The shear stress at the peak load and the unit is MPa.

Evolution of the Fracture Process Zone
The result difference in rock fracture toughness test obtained by different specimen types has attracted the attention of many scholars. Most scholars believe that the fracture process zone (FPZ) near the crack tip is responsible for variation of fracture toughness [35][36][37][38][39]. Therefore, it is necessary to analyze distribution and evolution of the specimen's FPZ during the loading process, in order to provide more understanding of rock fractures.
The damaging behavior and plastic strain law in the numerical model are shown in Figure 10. During the loading process, the increase of stress causes the elastic strain ε e of material. When the stress of a point reaches the tensile strength σ t , the elastic strain ε e reaches the maximum ε peak e . The strain after the peak elastic strain ε peak e is plastic strain ε p . The plastic strain will continue to increase when the energy stored is not enough to produce cracks. The ε f is the strain when the material is fractured. The area G f under the curve is the energy required for material fracture. When the energy stored is enough to destroy the mesh, the crack will appear in the plastic strain zone. curve is the energy required for material fracture. When the energy stored is enough to destroy the mesh, the crack will appear in the plastic strain zone. The plastic strain in the numerical model is shown in Figure 11. The area where the strain reaches about 0.5% (0.005) is considered damaged and can be defined as FPZ. The FPZ is first produced around the tip of the notch and then develops vertically upwards. When the specimen is subjected to a load of 4357 N (peak load), FPZ is obvious (red part as shown in Figure 11b). The plastic strain is about 5 mm, and no macroscopic fractures can be observed at this time. When FPZ (light blue part) grows to about 12 mm as shown in Figure 11d, the tip of the notch began to produce macro fractures inside FPZ. Then, the FPZ continues to develop with loading, and the macro fracture continues to occur inside FPZ. When FPZ penetrates the specimen, the bearing capacity of specimen has been reduced to zero. The plastic strain in the numerical model is shown in Figure 11. The area where the strain reaches about 0.5% (0.005) is considered damaged and can be defined as FPZ. The FPZ is first produced around the tip of the notch and then develops vertically upwards. When the specimen is subjected to a load of 4357 N (peak load), FPZ is obvious (red part as shown in Figure 11b). The plastic strain is about 5 mm, and no macroscopic fractures can be observed at this time. When FPZ (light blue part) grows to about 12 mm as shown in Figure 11d, the tip of the notch began to produce macro fractures inside FPZ. Then, the FPZ continues to develop with loading, and the macro fracture continues to occur inside FPZ. When FPZ penetrates the specimen, the bearing capacity of specimen has been reduced to zero.
When the specimen is subjected to a load of 4357 N (peak load), FPZ is obvious (red part as shown in Figure 11b). The plastic strain is about 5 mm, and no macroscopic fractures can be observed at this time. When FPZ (light blue part) grows to about 12 mm as shown in Figure 11d, the tip of the notch began to produce macro fractures inside FPZ. Then, the FPZ continues to develop with loading, and the macro fracture continues to occur inside FPZ. When FPZ penetrates the specimen, the bearing capacity of specimen has been reduced to zero. From the above analysis, it can be known that FPZ is ahead of macro fractures, and macro fractures are generated after peak load. When FPZ penetrates the specimen, the bearing capacity of specimen has been reduced to 0, and the macro fracture has not completely penetrated the specimen at this time. This indicates FPZ control fracture behavior of the specimen.

The Crack Propagation Velocity
The crack propagation velocity is an important piece of information in specimen failure. In order to obtain the crack propagation velocity, the crack length was measured every 30 microseconds. Figure 12 shows a relationship between the crack length and the time increment. It can be calculated that the average crack propagation velocity is about 111 m/s. The crack propagation process includes two stages according to crack propagation velocity, namely the rapid crack initiation development stage (A to B) and the final crack splitting stage (B to C). The crack propagates rapidly at the crack initiation development stage and velocity is about 267 m/s, which is also the main stage of the crack formation. The crack propagation velocity is slowed down at the final splitting stage. The crack propagation experienced fast and slow stages in the whole process. From the above analysis, it can be known that FPZ is ahead of macro fractures, and macro fractures are generated after peak load. When FPZ penetrates the specimen, the bearing capacity of specimen has been reduced to 0, and the macro fracture has not completely penetrated the specimen at this time. This indicates FPZ control fracture behavior of the specimen.

The Crack Propagation Velocity
The crack propagation velocity is an important piece of information in specimen failure. In order to obtain the crack propagation velocity, the crack length was measured every 30 microseconds. Figure 12 shows a relationship between the crack length and the time increment. It can be calculated that the average crack propagation velocity is about 111 m/s. The crack propagation process includes two stages according to crack propagation velocity, namely the rapid crack initiation development stage (A to B) and the final crack splitting stage (B to C). The crack propagates rapidly at the crack initiation development stage and velocity is about 267 m/s, which is also the main stage of the crack formation. The crack propagation velocity is slowed down at the final splitting stage. The crack propagation experienced fast and slow stages in the whole process. time increment. It can be calculated that the average crack propagation velocity is about 111 m/s. The crack propagation process includes two stages according to crack propagation velocity, namely the rapid crack initiation development stage (A to B) and the final crack splitting stage (B to C). The crack propagates rapidly at the crack initiation development stage and velocity is about 267 m/s, which is also the main stage of the crack formation. The crack propagation velocity is slowed down at the final splitting stage. The crack propagation experienced fast and slow stages in the whole process.

Conclusions
The mode I fracture process of the semi-circular bend (SCB) rock specimen is analyzed by experimental and numerical methods in this work. The main conclusions were obtained as follows:

Conclusions
The mode I fracture process of the semi-circular bend (SCB) rock specimen is analyzed by experimental and numerical methods in this work. The main conclusions were obtained as follows: (1) The notch tip of is a SCB specimen, often an arc, which may lead to the change of the crack initiation position. In the numerical model, there is an angle at the notch tip to ensure that the crack initiation position is in the center of tip. The FDEM numerical calculation method can well simulate the mode I failure mode and the stress-strain curve of the SCB specimen. It can provide important information of the mode I failure process, such as the stress field and strain field, which is difficult to obtain in the laboratory test, and provides powerful help for the study of mode I fracture.
(2) During the loading process, the tensile stress concentration zone of the SCB specimen generates and grows at the notch tip of the specimen. The tensile stress concentration zone of specimens presents a heart shape at peak load. After macro cracks were generated, the maximum tensile stress concentration zone moves upward with crack propagation.  Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The influence of the mesh size on the loading curves is shown in Figure A1. It can be seen that the loading curves of different mesh sizes are almost identical before the peak, and there is a small difference at the peak point. The mesh size of 1 mm can basically reflect the mechanical properties and failure of SCB specimens.
Generally, the smaller the mesh size can make, the finer the result. However, the crack propagation velocity of SCB specimen is fast. It needs to output a graphic data every 30 µs in the numerical model to capture the crack propagation process, resulting in a total of 10,000 graphic data outputs in the whole calculation process. The output of a large number of graphic data and too small mesh size will lead to numerical calculation and very slow post-processing.
Considering the running speed of the computer, in order to obtain the crack propagation process and reduce the influence of mesh size on the loading curve, the 1 mm mesh size and outputting 10,000 graphic data are selected in the numerical model.