Abstract
Investigating the crack initiation stress of rocks is vital for understanding the gradual damage process of rocks and the evolution law of internal cracks. In this paper, the particle flow code method is used to conduct biaxial compression tests on a marble model with an elliptical crack under different confining pressures. According to the evolution status of microcracks in the rock during compression, four characteristic stresses are defined to reflect the gradual damage process of the marble. Two different methods are used to obtain crack initiation stress of rocks, and the calculation results are compared with those based on Griffith’s strength theory to verify the accuracy of this theory under compressive stress. Based on the numerical simulation results, the evolution law for the strength parameters of marble with the degree of damage is described. According to the proportional relationship between the peak stress and crack initiation stress, a new method for predicting the initiation stress is proposed, whose effectiveness is verified. Overall, the results of this study can serve as a useful guide for solving the important problems of slab cracking and rockburst encountered in underground space engineering.
1. Introduction
With the growing utilization of underground space, the construction of deep-buried tunnels and mining of deep mineral resources has received considerable attention. The excavation of deep rocks can cause the brittle instability of the surrounding rock and even rockburst because the deep rock mass is in a complex environment such as high in situ stress [1]. Crack initiation stress () of rock is defined as the stress corresponding to the initiation of microcracks in the rock. Stress concentration occurs in the rock when its stress is larger than the crack initiation strength, leading to cracks and failure of the rock, and the strength of rock also decreases rapidly [2]. Therefore, it is of great significance to study the gradual damage process of rocks, understand the evolution process of internal cracks, and determine the initiation stress of rocks.
Several theoretical and experimental studies have focused on an in-depth investigation of of rocks. Martin et al. conducted uniaxial compression tests on Lac du Bonnet granite and found that and damage stress () were 0.4 and 0.8 times of the peak strength of the rock, respectively [3, 4]. Everitt and Lajtai experimentally studied the influence of the internal structure of rock on and and proved that the strength of the rock is significantly affected by its internal structure [5]. Zhu et al. performed conventional triaxial tests on the Three Gorges granite to observe that was consistent with the variation range of confining pressure, and was generally maintained at 30% to 50% of the peak stress () [6]. Liu et al. established crack initiation stress criteria for Jinping marble through acoustic emission tests and compared the results with those based on field test [7]. Cai et al. correlated the disturbance zone of surrounding rock with by using on-site acoustic emission test and microseismic monitoring [8]. Nicksiar and Martin used a database of 336 samples and examined the effects of geology and loading conditions on crack initiation stress [9]. Wu et al. introduced the investigation of AE cumulative count to identify the crack stress thresholds in the rock damage process [10]. Gong et al. carried out the USLUC test of ten rocks based on the principle of LURR theory and proposed a method for determining the crack compaction and crack damage stress thresholds of rocks [11, 12]. Until now, most researchers have primarily conducted indoor experiments to examine the gradual damage process of rocks. However, the crack propagation of rocks is an extremely complicated process, and rock mass is heterogeneous with natural cracks. Indoor experiments are random and require huge manpower as well as material resources. Numerical simulation methods are faster and more flexible than indoor experiments, and the corresponding results are not affected by environmental factors. With the advancement of computer technologies, several numerical simulation methods have been proposed to investigate the gradual damage process and to obtain of rocks. Liu et al. proposed a constitutive model by considering the gradual damage of granite and implemented the model into the large-scale finite element analysis software ABAQUS to simulate indoor compression tests of granite under different confining pressures. The simulated stress-strain curve was in good agreement with the test curve, and the model could reveal the prepeak nonlinear mechanical behavior of brittle rocks such as granite [13]. In the traditional finite element method (FEM), the direction of crack propagation must be preset, which causes the calculation results to deviate from the actual situation, so it needs to be remeshed each time the crack propagates. On the contrary, the extended finite element method (XFEM) employs an integral method to obtain the stress intensity factor at the crack tip and adopts the maximum energy release rate criterion to determine the propagation and direction of the crack. During the calculation process, the description of the discontinuous field is completely independent of the grid boundary, which is highly beneficial for dealing with fracture problems. Rabczuk et al. proposed a new method for treating crack growth by particle methods which is based on the cracking-particle method where the crack is treated as a collection of cracked particles. And this method split particles where the cracking criteria are met and introduce crack segments instead of enriching the particle and adding additional degrees of freedom [14, 15]. Ren et al. presented a dual-horizon peridynamics formulation which allows for simulations with dual-horizon with minimal spurious wave reflection and completely solved the “ghost force” issue [16, 17]. Boundary element method (BEM) has been applied in engineering calculations since 1960. It is considered to be more accurate and efficient than the finite element method for the problems with a large gradient of boundary variables, such as stress concentration problems or crack problems with singular boundary variables. In addition, the grid needs to be set only on the calculation boundary and the crack surface in BEM, which facilitates the advantages of small data processing workload and low memory consumption. However, the BEM is not suitable to solve the heterogeneous or elastoplasticity problems. The above numerical methods are typical representatives of continuum mechanics method. In fact, several discontinuities such as joints and cracks exist in the rock. Most continuum methods are based on some assumptions and simplifications, and they cannot be used to accurately describe the mechanical properties of rock from the viewpoint of continuum mechanics. In contrast, the discontinuous method is beneficial compared to continuum mechanics methods for the simulation of discontinuities in rocks. The discrete element method (DEM) proposed by Cundall [18] is a typical discontinuous method, which has been developed to a mature commercial numerical calculation software called particle flow code (PFC) [19]. PFC can be used to investigate the fracture damage of rock and soil from a microperspective, thereby providing unique advantages in the simulation of the rock fracture and crack propagation.
Although numerical simulation methods are frequently used to examine the gradual damage and crack propagation processes of rocks, only a few studies based on DEM have been reported. In this paper, Griffith’s strength theory is first introduced, and then PFC2D is used to construct a simulation model for marble with prefabricated elliptical crack. In addition, the biaxial compression test is conducted under different confining pressures, and of the rock is obtained by the volumetric strain versus the axial stress curve and the number of cracks. The calculation results are compared with those based on the crack initiation strength criterion deduced from Griffith’s strength theory, which verifies the accuracy and reliability of Griffith’s strength theory for predicting of rock under compressive stress conditions. Subsequently, based on the numerical simulation results, the evolution law for the strength parameters of marble as a function of the degree of damage is presented, and a new method for predicting is proposed based on the proportional relationship between and . Overall, this study provides useful insights for solving the problems of slab cracking and rockburst encountered in underground space engineering.
2. Griffith’s Strength Theory
Rock strength theory is a vital part of rock mechanics, which explains the effect of rock parameters, external loads, and environmental factors on the strength of rock when it loses stability under complex stress conditions. The traditional Mohr-Coulomb strength theory regards rock as a continuous homogeneous medium without fractures. Although this theory can describe the mechanism of plastic failure of rocks, the rocks experience many violent and complex tectonic movements during their long geological ages, and many discontinuous interfaces such as faults, joints, and fissures destroy their original integrity, so the rocks cannot be treated as an ideal continuous medium. Besides, most rocks exhibit brittle failure under load. Hence, these rock strength theories have certain limitations, and they can explain neither the failure mechanism of rocks with microcracks nor the brittle failure mechanism of rocks [20]. Griffith’s strength theory solves this problem well, which regards materials as discontinuous media and fully considers the difference between rocks and other mechanical media. Further, the discontinuities, that is, internal cracks of the rocks, are considered in this theory. Consequently, it provides an effective way to study the failure mechanics of rock masses.
Considering that the actual tensile strength of materials is usually much lower than that predicted theoretically, Griffith [21] suggested that the fracture damage of brittle materials is a result of the propagation of existing cracks. In the surrounding area of these cracks, especially at the crack tip, stress concentration occurs under the action of an external force, and once the tangential tensile stress becomes highly concentrated near the crack tip and reaches the molecular cohesive strength value of the material, brittle damage occurs in a certain direction. In addition, of the material has been determined under tensile stress by taking a flat elliptical crack as an example, which was then extended to the biaxial compression loading test [22]. Hoek verified Griffith’s crack initiation strength criterion of rocks derived from an elliptical single crack. However, he suggested that Griffith’s strength theory is only applicable when the minimum principal stress is zero or equal to tensile stress [23].
In Griffith’s strength theory, to determine the stress around the crack, it is assumed that there are microcracks in a flat elliptical shape within the material with a small axial ratio, as shown in Figure 1. Adjacent cracks are supposed to not affect each other, and the elliptical crack is treated as a single hole in a semi-infinite elastic medium. Further, the ellipse and the stress system acting on its surrounding materials are described using a two-dimensional problem [21]. The stress state of the crack edge iswhere and are the tangential stress and radial stress at the edge of the crack, respectively.
The extremum principle is used to calculate the value and position of the maximum dangerous stress around the elliptical crack, and the direction and stress of the long axis of the most dangerous crack are determined subsequently. Finally, the obtained ultimate stress and the uniaxial tensile strength are compared, and Griffith’s strength criterion is established as follows:where and are the maximum and minimum principal stress, respectively; represents the tensile strength of rocks.
3. Calculation Model and Results
3.1. Establishment of the Calculation Model and Calibration of Mesoscopic Parameters
Figure 2 shows a schematic of the prefabricated elliptical crack model established by servo control in the PFC. The model size is 100 mm × 100 mm. There is a prefabricated elliptical crack in the center of the model, which is 10 mm long. The ratio of the short half axis to the long half axis of an ellipse is the axial ratio. In materials such as rock mass, it can be considered that the ellipse has a very small axial ratio; that is, the shape of the ellipse is extremely flat. The model contains 24747 particles, and the particle radius R is uniformly distributed from 0.5 to 0.83 mm.
Choosing appropriate mechanical parameters is the key aspect of numerical simulation. The PFC characterizes the macroparameters of rock by setting the mesoparameters of the particle and contact model. The mesoparameters are generally determined by the “trial and error” method, and many researchers use the results of several numerical simulations and laboratory experiments to obtain the relationship between mesoparameters and macroparameters: (1) the material’s macroelastic modulus has a linear relationship with the meso-Young’s modulus at each particle-particle contact and parallel-bond contact (assuming that the two modulus values are the same) and has a logarithmic relationship with the normal-to-shear stiffness ratio of the particle-particle and parallel-bond contacts (assuming that the two ratios are the same). (2) Poisson’s ratio is mainly determined by the particle stiffness, and there was a logarithmic relationship between them. (3) The compressive strength is determined by the particle bond strength ratio. We calibrate according to the law between microscopic parameters and macroparameters [24]. Here, the macroparameters of the marble of Jinping Hydropower Station are calibrated through the numerical experiment of uniaxial compression and Brazilian splitting [25].
The uniaxial compression test is used to obtain the uniaxial compressive strength, Poisson’s ratio, elastic modulus, and other mechanical parameters of the rock. The test model used in the simulation has a height of 100 mm and a width of 50 mm. Figure 3 shows a schematic of the uniaxial compression test, where the red part in the model denotes cracks produced in the rock during compression process, and the curve is the stress-strain curve. Through this curve, the uniaxial compressive strength and elastic modulus of the rock can be obtained.
The tensile strength of the rock can be obtained through the Brazilian splitting test. Figure 4 shows a schematic of this test, where the diameter of the test model is 50 mm. The red part in the model shows the distribution of cracks in the rock, and the curve is the load-strain curve in the process of rock compression. In the simulation of rock test of PFC, the loading rate and time step of loading plate will affect the simulation results. We follow the proper loading rate of 0.05 m/s and the time step of 1.1 × 10−8 step−1 for uniaxial compression test and Brazilian splitting test [26]. The tensile strength of the rock can be expressed as follows:where is the load corresponding to damage and L is the sample thickness.
The macromechanical parameters of Jinping marble obtained through numerical experiments and actual indoor tests are shown in Table 1. It is clear that the macromechanical parameters obtained using the two sets of tests are basically consistent with the indoor test parameters. The corresponding mesoparameters are shown in Table 2.
3.2. Stress-Strain Curve and Crack Propagation Distribution of Marble
Figure 5 shows the numerically simulated stress-strain curve obtained using the biaxial compression test of the marble with an elliptical crack under confining pressure of 0–50 MPa. It is evident that the confining pressure significantly affects the peak strength of marble. The peak strength of the rock increases as the confining pressure increases, while the stiffness of the rock is hardly affected by the confining pressure. The rock undergoes brittle damage when the confining pressure is zero or lower, and the damage gradually changes from brittle to ductile form with the increase in the confining pressure. The marble exhibits a transition from brittleness to ductility when the confining pressure increases to 40 MPa and shows plastic flow when it increases to 50 MPa.
Figure 6 shows the relationship between the three strains of the rock and the axial stress when the confining pressure is 40 MPa. Among them, the volumetric strain can be hardly measured, and hence it is obtained by the following approximate formula [27]:where , , and represent the volumetric, axial, and lateral strains of the rock, respectively.
The curves in Figure 6 can be divided into four stages from the beginning of compression to the peak strength within the rock. (1) Stage I indicates the closing stage of the original cracks, during which the original cracks in the rock begin to close under compression. (2) Stage II is the elastic compression stage and starts from the closing stress (), during which the rock begins to deform elastically. The axial stress is linearly proportional to the three strains, and no new cracks are formed in this stage. (3) Stage III indicates stable crack growth and starts from crack initiation stress (). In this stage, new cracks are formed and gradually increase. The axial stress basically remains linearly proportional to the axial strain, but the axial stress-lateral strain relationship begins to change from a straight line to a curve. This is attributed to the dilatancy effect caused by the lateral expansion of the main axial crack. At the end of this stage, a typical conjugate shear plane begins to form in the rock sample. (4) Stage IV represents the process of accelerated crack growth and starts from damage stress (). During this stage, the lateral strain increases faster than the axial strain, and the increased volume of the crack connecting and expanding exceeds the elastic deformation formed by compression, so the volume-strain curve begins to deflect to the left, and the axial stress-axial strain relationship also changes to nonlinear. The typical conjugate shear surface of the rock sample develops into a macroscopic shear surface; the number of cracks increases gradually. Further, the volume of the rock increases until the peak strength () of the rock is reached when macroscopic damage occurs.
Figure 7 shows the crack propagation diagram corresponding to each characteristic stress point of the rock during compression. Figure 7(a) shows the crack propagation diagram at the point . It is clear that the tensile stress is concentrated at the joint end of the prefabricated elliptical crack, so tensile stress damage occurs at this end, and a tensile stress crack is initiated. Figure 7(b) shows the crack propagation diagram at the point . Here, the rock experiences a steady growth of cracks. The tensile cracks initiated from the crack tips propagate along the direction of the maximum principal stress and develop into tensile wing cracks. Figure 7(c) shows the crack propagation diagram at the point . The crack propagates rapidly, and the rock is destroyed in a short period of time during the accelerated growth of the crack.
(a)
(b)
(c)
4. Results and Discussion
4.1. Measurement of
The biaxial compression test of rock is performed by moving the upper and lower walls for applying axial pressure to the rock with a loading rate of 0.05 mm/s. Based on a large number of tests, of rock is generally 30% to 50% of the rock peak strength. In this paper, two methods are used for the measurement for of the rock: volumetric strain versus axial stress curve and crack number versus axial stress curve.
4.1.1. Calculation of by Volumetric Strain versus Axial Stress Curve
Figure 8 shows the volumetric strain versus the axial stress curve of the rock in the biaxial compression test when the confining pressure is 40 MPa. In the elastic compression stage, the volume of the rock steadily increases, and the curve is linear. However, the relationship between stress and strain is nonlinear in other stages. Therefore, the entire stress-strain curve is divided into curved-straight-curved three segments. A reference line is used to determine , where the point at which the reference line starts to deviate from the straight corresponds to of the rock.
4.1.2. Calculation of Using Crack Number versus Axial Strain Curve
Potyondy and Cundall [28] established a parallel-bond model (PBM) to simulate the behavior of rock in which parallel bonding (used to characterize bonds) is used to cohere particles that are in contact with each other to form an aggregate and characterize the whole rock. A schematic of this model is shown in Figure 9. Parallel bonds can be considered as a series of springs with a certain stiffness, which act together on the contact of particles, which can withstand the action of force and moment and have a specific tensile and shear strength. The parallel bonding contact breaks when the force and moment generated at the contact exceed its strength. Contact fracture is used to represent the generation of new cracks. Figure 10 shows the number of cracks in the rock during biaxial compression. The stress is defined as when the number of microcracks is 1% of the microcracks at the peak strength [29]. In the stage of stable crack growth, the number of cracks increases slowly with the increase in axial strain, and the relationship between the number of cracks and the axial strain is approximately linear. After reaching the stage of rapid crack growth, the slope of this curve increases significantly.
4.2. Theoretical and Calculated Values of
To verify Griffith’s strength theory, the two methods described in the previous subsection are used to calculate the values of for the rock under different confining pressures, which are compared with the theoretical value. The theoretical and calculated values of under each confining pressure are shown in Figure 11 and Table 3. The values of calculated by the two methods are close to those obtained by Griffith’s strength criterion, which verifies the feasibility of this criterion. Further, it can be seen that the calculated value from the crack number curve is closer to the theoretical value than that from the volume-strain curve. The average value of is slightly higher and lower than the theoretical value when the confining pressure is less and greater than 30 MPa, respectively, and it is closest to the theoretical value when the confining pressure is 30 MPa.
4.3. Relationship between and and Confining Pressure
Table 4 shows and of the rock under different confining pressures and their ratio . It can be observed that is 35–51% of , which is consistent with the research result of Martin et al. [3], in which rock cracks when for Lac du Bonnet granite in Canada. When the confining pressure is less than 30 MPa, the ratio increases with the increase in the confining pressure. When it is greater than 30 MPa, the ratio remains constant at 0.51.
Figure 12 shows the relationship between and and confining pressure. It is clear that both the stresses are linearly proportional to the confining pressure, but changes more significantly than. This is mainly because the closure degree of primary microcracks increases with the increase in confining pressure, and the contact area of these microcracks under high confining pressure is larger than that under low confining pressure, which increases the friction strength and .
4.4. Evolution of Strength Parameters of Marble
Using the measured characteristic strength values of marble under different confining pressures, the strength parameters under different stages of crack evolution can be analyzed quantitatively. To characterize the effect of confining pressure, a simple linear criterion is generally used to predict of rock; that is, , where is a parameter related to the type of rock, is the confining pressure, and is the uniaxial compressive strength of the rock. In this equation, the coefficient of confining pressure is considered as 1; that is, it is assumed that the internal friction angle of the rock at the crack initiation stage is 0, which is obviously inconsistent with the actual situation. According to the previous subsection, can be expressed as
Under the action of confining pressure, the original cracks in the rock are in a closed state under compression, and the new cracks need to overcome the frictional force on the contact surface during the crack propagation process, so the coefficient of in equation (5) is not 1. The coefficient of confining pressure is related to the friction properties of the material, and hence the influence of friction should be considered while analyzing the crack evolution process under confining pressure. Rabczuk and Belytschko [15] considered the influence of friction for deriving the local at the end of the crack and obtained the crack initiation damage criterion at the microscopic level as follows:where , is the friction coefficient of rocks, and is the internal friction angle.
After simplifying equation (6), the coefficient of confining pressure can be obtained as follows:
The internal friction angles at the initiation stage () and peak stage () of the marble can be obtained as 12.7° and 28.5° by substituting the fitting results into equation (7) and using the Mohr-Coulomb criterion, respectively. is nearly half of . It can be concluded from the above results that the internal friction angle of the rock does not remain constant during the damage process but gradually increases as the degree of damage increases. The internal structure of the rock is dense, and it is generally difficult to generate relative sliding between mineral crystals, so the rock strength is primarily governed by the cohesion of the rock. The friction strength of the rock is not fully utilized, and the internal friction angle of the rock is small when the rock does not crack. When tensile damage of the rock occurs, it starts to crack and expand, and its crystal structure is destroyed. There is a tendency of squeezing and relative sliding between the mineral particles of the rock, and the frictional resistance between the crystals needs to be overcome when the cracks further expand, which further increases the effective friction resistance. At this time, both the internal friction angle and the friction strength of the rock gradually increase.
4.5. Axial Stress-Strain Ratio at and Points
The ratio () of and of rock can be expressed as follows:
The specific data of is shown in Table 4. The ratio of axial strain () corresponding to point to the axial strain () corresponding to point, that is, the axial strain ratio , is expressed as follows:
Table 5 shows the values of , , and under different confining pressures. It is clear that the variation trend of axial strain with the confining pressure is the same as that of the stress. When the confining pressure is less than 30 MPa, the ratio of the two increases as the confining pressure increases. When the confining pressure is greater than 30 MPa, the ratio of the two remains stable.
Figure 13 shows that the stress ratio is linearly proportional to the axial strain ratio of the rock under different confining pressures, indicating that the linear relationship is not affected by the confining pressure or the nature of the rock. The slope of this linear curve is approximately 1.218, which can be used to obtain the following equation:where
Here, is the slope of the line from the origin to , and is the slope of the line from the origin to the point.
As shown in Figure 14, the slope k of the straight line is obtained through the line connecting the origin and , and a straight line with slope of 1.218k is drawn through the origin. The intersection point with the stress-strain curve is the point. This method can be used to predict of the rock just through the stress-strain curve. Table 6 shows the predicted values of under different confining pressures through this method. The error Δ is defined as follows:where is the value of that is predicted by the present method and is the calculated value of . It is clear that the predicted values of are close to the actual calculated values, and the error ranges between 3.1% and 10.7%, which verifies the reliability of this method. It should be noted that this method used here is not suitable for prepeak cyclic loading, unloading, and double-humped conditions.
5. Conclusions
In this paper, PFC software was used to establish a calculation model for marble with prefabricated elliptical crack. We conducted biaxial compression tests under different confining pressures and compared the values of crack initiation stress () determined by different methods. Further, the calculation results were compared with those based on Griffith’s strength theory using a single elliptical crack. The main results of the study can be summarized as follows:(1)Crack initiation stress () of marble was obtained by the volume-strain curve versus the axial stress curve and the crack number versus the axial strain curve. The calculated values were consistent with those based on Griffith’s theory, which verified the applicability of Griffith’s strength criterion under compressive stress conditions. The calculated value from the crack number versus the axial strain curve was closer to the theoretical value than that from the volume-strain curve. The calculated values were slightly higher and lower than the theoretical values when the confining pressure was smaller and greater than 30 MPa, respectively.(2) of marble was between 35% and 52% of . This is consistent with the results of many studies. When the confining pressure was less than 30 MPa, the ratio increased with the increase in the confining pressure. When it was greater than 30 MPa, the ratio remained constant at 0.51.(3)Using the linear fitting equation, it was obtained that . The internal friction angle of the rock gradually increased as the degree of damage increased, and the evolution process of the internal friction angle affected the brittle failure tendency of the rock.(4)The stress ratio was linearly proportional to the axial strain ratio between and points. According to the proportional relationship between these points, the connection method was proposed and used to directly obtain of the rock from the stress-strain curve. This provided a simple, rapid, and accurate method for predicting .
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have conflicts of interest regarding the paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 51879245, 41920104007, and 11672360) and the Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan) (nos. CUGCJ1821 and CUG170645).