In order to investigate the feasibility and reliability of the three-dimensional particle flow method in simulating the type I fracture toughness test, four types of numerical samples were established by particle flow code PFC3D: straight crack three-point bending (SC3PB), edge cracked flattened semicircular disc (ECFSD), cracked chevron notched Brazilian disc (CCNBD), and edge cracked flattened ring (ECFR). Three models with different strength parameters (group A, group B, and group C) were established for each type, in which group A parameters are obtained from the concrete model, group B parameters are applied for simulating marble, and group C parameters are for granite. The type I fracture toughness and the failure form of each model are obtained by conducting the numerical test, and the curves of load versus displacement of loading point are recorded. The numerical test results show that, with the same strength parameter, the maximum difference in test results of each specimen type is 0.39 MPa·m1/2. The KIC of ECFR specimen is 0.13–0.28 MPa·m1/2 smaller than that of CCNBD specimen, and the KIC of ECFSD specimen is slightly higher than that of CCNBD sample. The KIC of SC3PB specimen is 0.06–0.21 MPa·m1/2 smaller than that of the CCNBD sample. When the loading rate is less than 0.01 m/s, the effect of loading rate on fracture toughness can be reduced to less than 0.1 MPa·m1/2.

1. Introduction

The research of theoretical and laboratory test for type I fracture toughness KIC of rock materials is relatively mature [1]. Typical test methods include three-point bending test, compact tensile test, and Brazilian splitting test, and the specimen types of SC3PB, CCNBD, and SR are commonly used.

Zhang et al. [2] found that the fracture toughness of the sample without notch is higher, and the discreteness in results is obtained. Cui et al. [3] reviewed the testing methods of type I fracture toughness of rocks, compared the results by different methods, and explained the causes of these differences. Wei et al. [4] think that the results of three-point bending test are more stable and reliable, and it is recommended to use the formula proposed by ASTM. Ayatollahi et al. [5] found that the larger the diameter specimen, the higher the type I fracture toughness. The research results of Meng et al. [6] and Yang et al. [7] show that as the center angle of the platform Brazil sample increases, the failure mode becomes complex. The most appropriate platform angle is 20°. Huang et al. [8, 9] study the influence of different particle sizes on load-displacement curves and failure modes by PFC2D. The results show that the generation of secondary cracks is mainly affected by particle sizes. The size effect in particle flow software simulation has been studied by Wong et al. [10]. It was found that the increase of particle size would cause the increment in fracture toughness and crack initiation stress. The deformation characteristics of sandstone and granite are researched under different stress path by Peng et al. [1114], Wang et al. [15], and Shang et al. [16], and the energy characteristics during crack propagation are studied. It is found that the most sensitive parameter for stage identification is the volumetric strain. With the increment of confining pressure, the development and connection direction of cracks are inhibited, and the length of the final fracture surface decreases. The test process and results of SR and CCNBD specimens are studied by Cui et al. [17, 18]. The results suggest significant size effects of CCNBD specimens. As the specimen diameters increase, the variance in fracture toughness KIC of SR and CCNBD specimen becomes smaller. The facture surface of SR specimen is rougher than that of CCNBD. It is also found that fracture toughness test results can be more consistent by the specimen with larger diameter than the ISRM suggested “minimum effective diameter 75 mm.” Based on the boundary effect theory, the fracture toughness of rock is studied by Guan et al. [19]. The real material parameters without size effect of different rock types are constructed, and the results are used to predict the fracture trend.

Experimental study on the fracture toughness of CCNBD specimen and the size effect was conducted by Wu et al. It is proposed that the results can be modified by geometric shape function, and then the real fracture toughness of rock can be obtained. In semicircular bending (SCB) test, the support type influence on rock fracture toughness is researched by Bahrami et al. [20]. The results of finite element method and laboratory tests show that different support types and the friction between the specimen and the bottom supports have great influence on results. The fracture toughness KIC of lapilli-ash tuff is researched by Wong et al. [21] with two typical semicircular bending specimens. The fracture toughness KIC measured by semicircular bending (SCB) method is found to be lower than that using cracked chevron notched semicircular bending (CCNSCB) method. The CCNBD method produces more scattered results.

In laboratory tests of rock fracture toughness KIC, the sample processing is relatively difficult, and the test results are scattered. Therefore, numerical simulation method is widely used to test rock fracture toughness. Most of the numerical tests are simulated by two-dimensional software, and three-dimensional numerical calculation method has not been widely discussed. In this study, four types of sample are selected, and three-dimensional particle flow numerical simulation is used to test the type I fracture toughness of the samples. The results are compared and analyzed, which provides a reference for numerical test of fracture toughness.

2. Establishing a Numerical Model

2.1. Microparameters of the Model

Particle flow software PFC is widely used to simulate the deformation and failure process of elastic-plastic materials, such as rock, soil, and concrete, which can show the mechanical properties and failure mechanism from a micro perspective. By adjusting the parameters of particle and bond model, the mechanical characteristics of numerical model can be similar to actual materials. Failure process can be obtained by monitoring the number and location of microcracks and the stress in model. In order to study the applicability of numerical simulation in rock materials with different strength, three groups of micro parameters (group A, B, and C) are selected for calculation and analysis (see Tables 13), among which the strength parameters are lowest in group A and highest in group C.

The parameters in group A are calibrated according to the results of direct shear test in the laboratory. The parameters in group B are the micro parameters for marble taken by Huang et al. [22]. The parameters of group C are modified based on the Brazil split test curve for granite. The uniaxial compressive strength of numerical models can be used to calibrate parameters, and the test results show that the strength is 78.4 MPa with group A parameters, 106.8 MPa with group B parameters, and 142.1 MPa with group C parameters.

As PFC2D is used in reference [17], the rationality of micro parameters in three-dimensional particle flow simulation needs to be further tested and modified. The details of parameter calibration in group A are as follows:(1)Mix cement, sand, and water evenly by the mass ratio of 1 : 1 : 0.4, put them into a cubic mold sized 15 cm × 15 cm × 15 cm (as shown in Figure 1), and vibrate them tightly. After placing for 1 day, demold them and cure them at room temperature for 28 days; then, a similar model sample can be obtained [23], whose physical and mechanical properties are close to actual rock.(2)Carry out uniaxial compress test and direct shear test on samples under different normal stress (as shown in Figure 2); then, normal stress and shear stress during the test are recorded and plotted in Figure 3. The test results show that the cohesion c of the sample is 5.09 MPa, and the internal friction angle φ is 44.8° (as shown in Figure 3). The uniaxial compressive strength of samples is 82.3 MPa.(3)Establish PFC3D numerical models of the same size, and carry out numerical test with the micro parameters of Huang et al. [13]. Modify these parameters until the test results are closer to those of the laboratory test; then, plot the results in Figure 3. The test results show that the cohesion c of the sample is 4.93 MPa and the internal friction angle φ is 46.7° (as shown in Figure 3); take the micro parameters as group A.

The calibration and modification of parameters in group C are as follows:(1)Cylindrical specimens sized ϕ50 mm × 25 mm are made with intact granite; then, a series of Brazil disk split tests and uniaxial compress tests are conducted with these specimens (as shown in Figure 4). The load-displacement curves of the test are recorded in Figure 5. According to the procedures, during the Brazil tests, the loading rate is controlled as 0.1 MPa/s. The uniaxial compress strength of granite samples is 142.3 MPa.(2)Establish the numerical model of particle flow and carry out the numerical test according to step (1), and then adjust the micro parameters of the model repeatedly, until the load-displacement curves of numerical test agree well with those of laboratory tests (as shown in Figure 5). The parameters are taken as group C.

Based on the above micro parameters, 12 numerical models of SC3PB, ECFSD, CCNBD, and ECFR were established, respectively, among which SC3PB and CCNBD models have been widely used, and ECFSD and ECFR are new type models studied by Zhang [24]. The size and loading form of models are shown in Figures 69. According to Potyondy et al. [25], the particle size can be selected as 3–4% of the model size, and numerical tests are carried out for each model.

2.2. Controlling the Loading Rate

Displacement loading is used in the test. In order to select an optimal loading rate, the SC3PB and CCNBD sample models are established with parameters in group A, and the test is executed under three loading rates of 0.05 m/s [9], 0.01 m/s, and 0.002 m/s [26]. The strength curves of the two samples under different loading rates are obtained, respectively (as shown in Figures 10 and 11).

Elastic stages of the curves, under different loading rates, are almost the same. The slower the loading rate, the lower the peak strength of the curve. For SC3PB sample, KIC is calculated by substituting the first load platform in the curve into formulas (1) and (2) [27]:where P is the load value, kN; s is the nominal span, mm; W is the specimen height, mm; B is the specimen thickness, mm; and a is the prefabricated crack length, mm.

For CCNBD samples, KIC is calculated by substituting the maximum load [5] in test by formulas (3)–(6) [28]:where is critical strength factor; α0 = a0/R, the dimensionless initial crack length; α1 = a1/R, the dimensionless maximum crack length; αB = aB/R, the dimensionless specimen thickness; and Pmax, the local maximum load, kN.

The fracture toughness of the samples under various loading rates is listed in Table 4. The results of SC3PB samples, with the error of 0.16 MPa·m1/2, are greatly affected by loading rate. The results of CCNBD samples are almost close to the error of 0.08 MPa·m1/2. With the decrease of loading rate, the influence of loading rate becomes smaller. The slower the loading rate is, the closer the test is to the quasi-static loading. With the same displacement, there is less damage in the sample, and the results are more accurate. Therefore, theoretically, the loading rate should be controlled relatively low, but considering the calculation efficiency of numerical simulation, the loading rate should be increased properly under the condition that the results are accurate enough. According to the above test, the loading rate of this test is controlled as 0.01 m/s.

3. Results of the Numerical Test

3.1. Three-Point Bending Test Results

For numerical test, first of all, the rationality and validity of the test results should be preliminarily judged. Through monitoring the generation and distribution of microcracks in the model, the failure form of the sample can be observed (as shown in Figures 1214. The black part in the figure represents the microcracks produced by the bond failure between particles. Among them, Figures 12(a), 13(a), and 14(a) show the crack initiation, Figures 12(b), 13(b), and 14(b) show the crack propagation, and Figures 12(c), 13(c), and 14(c) show that the specimen is damaged when the crack propagates to a certain stage. It is found that, in the samples with different strength, the microcracks are first generated from the tip of the preformed groove and gradually expand, indicating that the failure of the model is caused by the growth of the groove, so the numerical test is reasonable and effective.

Then, record the curve of load P-crack opening displacement V in the test (as shown in Figure 15). According to the test results, PA = 1.11 kN, PB = 2.99 kN, PC = 3.27 kN; substituting these values into formula (1), the KIC of model A1 is calculated as 0.92 MPa·m1/2, KIC of model B1 is 2.49 MPa·m1/2, and KIC of model C1 is 2.72 MPa·m1/2.

3.2. ECFSD Test Results

During the test, the distribution of microcracks in the model is shown in Figures 1618. It is found that a small number of microcracks, at the crack initiation stage, appear at the loading point and the prefabricated crack end. Then, the prefabricated crack propagation stops after about 5 mm, and a large number of microcracks appear at the two loading points. Crushing along the loading direction is the main failure mode of the sample.

In B2 model, the crack occurs at the loading points firstly and then occurs at the end of precrack (as shown in Figure 17). The fracture in the direction of diameter caused by stress concentration at both ends of the specimen is the main cause of specimen failure.

Record the curve of load P-displacement V in the test (as shown in Figure 19), and substitute the local minimum value on the curve into the following formula [19]:where Ymax is the dimensionless stress intensity factor [17], Ymax = 1.0756; R is the radius of sample, mm; Pmin is the local minimum load, kN; and the other symbols are the same as before.

According to the test results, PAmin = 6.22 kN, PBmin = 14.9 kN, PCmin = 16.3 kN; then, KIC of models can be calculated. The results show that KIC = 1.11 MPa·m1/2 for A2 model, KIC = 2.66 MPa·m1/2 for B2 model, and KIC = 2.91 MPa·m1/2 for C2 model. The observation of the loading process shows that the local minimum value Pmin of the load curve corresponds to the crack initiation and propagation stage of the prefabricated crack. At this time, the specimen is not completely damaged. As the loading continues, the load curve increases, and the second peak value may exceed the first.

3.3. CCNBD Test Results

The distribution of microcracks in the models is shown in Figures 2022. It is found that microcracks occur at the loading points of the sample and the end of precut groove firstly, then the precut groove begins to propagate, and the microcracks run through the sample along the diameter direction finally. A number of microcracks near the loading platform occur because of the local fracture caused by the stress concentration.

The curve of load P-displacement is recorded in the test (as shown in Figure 23), and the maximum load in the test is substituted into formulas (3)–(6), for calculation.

According to the test results, PAmax = 12.2 kN, PBmax = 18.4 kN, PCmax = 20.4 kN; then, KIC of models are calculated. The results show that KIC = 1.13 MPa·m1/2 for A3 model, KIC = 2.55 MPa·m1/2 for B3 model, and KIC = 2.83 MPa·m1/2 for C3 model.

3.4. ECFR Test Results

The distribution of microcracks in the model is shown in Figures 2426. In the process of loading, the vertical microcracks are produced on the inner wall of the samples firstly; then, the cracks propagate gradually from the inside to the outside along the loading direction and run through the samples. At the same time, the precut cracks also propagate to the inside of the samples. Tensile failure also occurs in the right half of the ring.

The curve of load P-displacement is shown in Figure 27. According to reference [17], the local minimum load Pmin is substituted into formula (7), where Ymax = 1.0468, PAmin = 5.71 kN, PBmin = 13.0 kN, PCmin = 15.2 kN, and the results are as follows: KIC = 1.00 MPa·m1/2 for A4 model, KIC = 2.27 MPa·m1/2 for B4 model, and KIC = 2.65 MPa·m1/2 for C4 model.

4. Analysis of Test Results

In each group of samples, SC3PB and CCNBD samples are damaged by the propagation of precut cracks, while ECFSD and ECFR samples are damaged by not only the propagation of precut cracks, but also the coalescence of vertical cracks. The load-displacement curve of each specimen will decrease slightly and then increase again, which is caused by stress release after crack initiation. By comparing the fracture toughness of each model (see Table 5), it can be found that, in group A, the difference between the four sample types is 0.21 MPa·m1/2. In group B, the difference between the four sample types is 0.39 MPa·m1/2. In group C, the difference between the four sample types is 0.26 MPa·m1/2. On the whole, the fracture toughness of the ECFR samples is smaller than that of the ECFSD and CCNBD samples (Tables 4 and 5).

By comparing the distribution of microcracks in each model, it can be found that, in ECFSD and CCNBD samples, both the crack growth and the fracture of the loading point can be obtained. The short crack propagation distance in the ECFSD specimen is not the main reason for specimen failure. It is because the crack tip is within the width covered by specimen loading platform as precut crack propagating. In the contact diagram (as shown in Figure 28), the black line represents the pressure between particles, and the thickness of the line represents the magnitude of force. It is explicit that the specimen is compressed along the vertical direction, and the precut crack tip will not undergo tensile failure anymore, so the crack will not further expand. Under a lower loading rate, there are fewer microcracks outside the propagation path in the CCNBD specimen, which makes the results relatively accurate.

The ECFR specimen can be regarded as a platform Brazilian disk specimen with edge crack and central circular hole. Compared with the disk or half-disk specimen, the ring-like specimen is more prone to compression deformation.

The calculated results are closely related to the ratio of internal and external radius r/R, the ratio of crack length to external diameter a/R, and other geometric parameters. There are many microcracks outside the crack propagation path, which has a certain influence on the test results.

5. Conclusions

The fracture toughness test of rock with 4 types of specimens is researched by the numerical simulation method of particle flow. The failure forms and load-displacement curves are analyzed and compared. The conclusions in this research are listed as follows:(1)The results of the numerical test are reasonable and effective. The maximum difference between the test results of different samples with the same strength parameter is 0.39 MPa·m1/2.(2)When the loading rate is reduced to 0.01 m/s, the effect of loading rate on fracture toughness can be reduced to less than 0.1 MPa·m1/2. So, the loading rate of 0.01 m/s is reasonable.(3)During the loading process, the microcracks occur in multiple areas in ECFR specimens. The test results of ECFR specimens are 6%–11% smaller than those of CCNBD specimens. For ECFSD specimens, there are many microcracks generated along the loading direction as the propagation of prefabricated crack. The test load is larger than the other specimens, so the fracture toughness of ECFSD specimens is 0.08–0.11 MPa·m1/2 larger than that of CCNBD specimens.(4)The SC3PB and CCNBD specimens are spilt along the loading direction because of the propagation of prefabricated crack, and the test results of SC3PB are 0.06–0.21 MPa·m1/2 less than those of CCNBD specimens.

Data Availability

The article data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This work was supported by the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant nos. 19KJD410001 and 18KJB440002), the Science and Technology Project of Housing and Construction in Jiangsu Province (Grant No. 2018ZD199), and the Science and Technology Project of Changzhou (Grant No. CJ20190020).