Discrete Element Modeling of Crack Initiation Stress of Marble Based on Griffith’s Strength Theory

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.


Introduction
With the growing utilization of underground space, the construction of deep-buried tunnels and mining of deep mineral resources has received considerable attention. e 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 (σ ci ) 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]. erefore, 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 σ ci of rocks. Martin et al. conducted uniaxial compression tests on Lac du Bonnet granite and found that σ ci and damage stress (σ c d ) 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 σ ci and σ c d 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 ree Gorges granite to observe that σ ci was consistent with the variation range of confining pressure, and σ ci was generally maintained at 30% to 50% of the peak stress (σ cu ) [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 σ ci 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 σ ci 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. e 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. e 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. e 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 σ ci of the rock is obtained by the volumetric strain versus the axial stress curve and the number of cracks.
e 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 σ ci 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 σ ci is proposed based on the proportional relationship between σ cu and σ ci . Overall, this study provides useful insights for solving the problems of slab cracking and rockburst encountered in underground space engineering.

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. e 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, σ ci 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]. e stress state of the crack edge is where σ θ and σ r are the tangential stress and radial stress at the edge of the crack, respectively. e 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 σ 1 and σ 3 are the maximum and minimum principal stress, respectively; S t represents the tensile strength of rocks.

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. e model size is 100 mm × 100 mm. ere is a prefabricated elliptical crack in the center of the model, which is 10 mm long. e 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. e 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. e PFC characterizes the macroparameters of rock by setting the mesoparameters of the particle and contact model. e 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) e 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]. e uniaxial compression test is used to obtain the uniaxial compressive strength, Poisson's ratio, elastic modulus, and other mechanical parameters of the rock. e Advances in Civil Engineering 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. rough this curve, the uniaxial compressive strength and elastic modulus of the rock can be obtained. e 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. e 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]. e tensile strength of the rock can be expressed as follows: where p t is the load corresponding to damage and L is the sample thickness. e 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. e corresponding mesoparameters are shown in Table 2.

Stress-Strain Curve and Crack Propagation Distribution of
Marble. Figure 5 shows the numerically simulated stressstrain 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. e peak    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 ε v , ε 1 , and ε 2 represent the volumetric, axial, and lateral strains of the rock, respectively. e 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 (σ cc ), during which the rock begins to deform elastically. e 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 (σ ci ). In this stage, new cracks are formed and gradually increase. e 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.
is 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 (σ c d ). 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. e 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 (σ cu ) 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 σ ci . 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 σ c d . Here, the rock experiences a steady growth of cracks. e 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 σ cu . e crack propagates rapidly, and the rock is destroyed in a short period of time during the accelerated growth of the crack.

Measurement of σ ci .
e 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, σ ci of rock is generally 30% to 50% of the rock peak strength. In this paper, two methods are used for the measurement for σ ci of the rock: volumetric strain versus axial stress curve and crack number versus axial stress curve.

Calculation of σ ci 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. erefore, the entire stress-strain curve is divided into curved-straightcurved three segments. A reference line is used to determine σ ci , where the point at which the reference line starts to deviate from the straight corresponds to σ ci of the rock.

Calculation of σ ci 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.  Advances in Civil Engineering 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. e stress is defined as σ ci 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.

eoretical and Calculated Values of σ ci .
To verify Griffith's strength theory, the two methods described in the previous subsection are used to calculate the values of σ ci for the rock under different confining pressures, which are compared with the theoretical value. e theoretical and calculated values of σ ci under each confining pressure are shown in Figure 11 and Table 3. e values of σ ci 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. e average value of σ ci 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. Table 4 shows σ ci and σ cu of the rock under different confining pressures and their ratio R σ . It can be observed that σ ci is 35-51% of σ cu , which is consistent with the research result of Martin et al. [3], in which rock cracks when σ ci � 0.4σ c 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 σ ci and σ cu and confining pressure. It is clear that both the stresses are linearly proportional to the confining pressure, but σ cu changes more significantly than σ ci . is 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 σ ci .

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 σ ci of rock; that is, σ ci − σ 3 � Aσ c , where A is a parameter related to the type of rock, σ 3 is the confining pressure, and σ c is the uniaxial compressive strength of the rock. In this equation, the coefficient of confining pressure is considered as 1; that Advances in Civil Engineering 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, σ ci 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 σ 3 in equation (5) is not 1. e 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 σ ci at the end of the crack and obtained the crack initiation damage criterion at the microscopic level as follows: where μ � tanφ, μ 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: e internal friction angles at the initiation stage (φ 0 ) 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. φ 0 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. e 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. e 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.
ere 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.

Axial Stress-Strain Ratio at σ cu and σ ci Points.
e ratio (R σ ) of σ ci and σ cu of rock can be expressed as follows:     Advances in Civil Engineering e specific data of R σ is shown in Table 4. e ratio of axial strain (ε ci ) corresponding to σ ci point to the axial strain (ε cu ) corresponding to σ cu point, that is, the axial strain ratio R ε , is expressed as follows: Table 5 shows the values of ε ci , ε cu , and R ε 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. e slope of this linear curve is approximately 1.218, which can be used to obtain the following equation: where Here, k is the slope of the line from the origin to σ cu , and k 1 is the slope of the line from the origin to the σ ci point.
As shown in Figure 14, the slope k of the straight line is obtained through the line connecting the origin and σ cu , and a straight line with slope of 1.218k is drawn through the origin. e intersection point with the stress-strain curve is the σ ci point. is method can be used to predict σ ci of the  Advances in Civil Engineering rock just through the stress-strain curve. Table 6 shows the predicted values of σ ci under different confining pressures through this method. e error Δ is defined as follows: where σ ci ′ is the value of σ ci that is predicted by the present method and σ ci is the calculated value of σ ci . It is clear that the predicted values of σ ci 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.

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 (σ ci ) determined by different methods. Further, the calculation results were compared with those based on Griffith's strength theory using a single elliptical crack. e main results of the study can be summarized as follows: (1) Crack initiation stress (σ ci ) of marble was obtained by the volume-strain curve versus the axial stress curve and the crack number versus the axial strain curve. e calculated values were consistent with those based on Griffith's theory, which verified the applicability of Griffith's strength criterion under compressive stress conditions. e calculated value from the crack number versus the axial strain curve was closer to the theoretical value than that from the volume-strain curve.
e calculated values were slightly higher and lower than the theoretical values when the confining pressure was smaller and greater than 30 MPa, respectively.
(2) σ ci of marble was between 35% and 52% of σ cu . is 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 φ 0 � φ/2. e 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) e stress ratio was linearly proportional to the axial strain ratio between σ ci and σ cu points. According to the proportional relationship between these points, the connection method was proposed and used to directly obtain σ ci of the rock from the stress-strain curve. is provided a simple, rapid, and accurate method for predicting σ ci .

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have conflicts of interest regarding the paper.