Investigation of Microcrack Propagation and Energy Evolution in Brittle Rocks Based on the Voronoi Model

The cracking of rock mass under compression is the main factor causing structural failure. Therefore, it is very crucial to establish a rock damage evolution model to investigate the crack development process and reveal the failure and instability mechanism of rock under load. In this study, four different strength types of rock samples from hard to weak were selected, and the Voronoi method was used to perform and analyze uniaxial compression tests and the fracture process. The change characteristics of the number, angle, and length of cracks in the process of rock failure and instability were obtained. Three laws of crack development, damage evolution, and energy evolution were analyzed. The main conclusions are as follows. (1) The rock’s initial damage is mainly caused by tensile cracks, and the rapid growth of shear cracks after exceeding the damage threshold indicates that the rock is about to be a failure. The development of micro-cracks is mainly concentrated on the diagonal of the rock sample and gradually expands to the middle along the two ends of the diagonal. (2) The identification point of failure precursor information in Acoustic Emission (AE) can effectively provide a safety warning for the development of rock fracture. (3) The uniaxial compression damage constitutive equation of the rock sample with the crack length as the parameter is established, which can better reflect the damage evolution characteristics of the rock sample. (4) Tensile crack requires low energy consumption and energy dispersion is not concentrated. The damage is not apparent. Shear cracks are concentrated and consume a large amount of energy, resulting in strong damage and making it easy to form macro-cracks.


Introduction
Rock fracture has always been one of the most concerning issues in the field of rock engineering practice, such as underground excavation [1,2], rock slope stability [3], seismic influence [4], and enhanced geothermal system [5], etc. The existing mechanical characteristic experiments can only explore the macroscopic failure and cannot directly explore the rock's compression fracture process. The initiation, accumulation, and development of internal microscopic cracks in the rock are essential for the irreversible damage and brittle fracture of the rock during compression [6,7]. Therefore, the same physical experiment condition is created by numerical simulation to reappear the stress-strain path of the rock material and to study the crack development process, which cannot be observed in the experiment.
There are many ways to explore the evolution of rock cracks using numerical simulation methods, typically classified as continuum methods and discrete methods [8]. The problem of the ratio of unconfined compressive strength to tensile strength [28,34]. The Voronoi model has the advantage of using simple logical combinations to achieve close to real rock materials.
In this study, the UDEC-Voronoi approach was used to build the numerical models for uniaxial/triaxial compressive and Brazilian split tests, and the microscopic parameters of contacts were obtained for four typical rock samples by calibration and inversion. Then, the uniaxial compressive tests of rock with four different compressive strengths were conducted. The characteristics of changes in the number, angle, and length of cracks in the process of rock failure and instability were monitored. The laws of crack development, damage evolution, and energy evolution were analyzed using the Voronoi model. The failure and instability mechanism of rock material under the discrete block model was discussed. It reveals the failure laws of rock materials at the micro-level and reproduces the whole process of crack evolution under natural conditions.

Voronoi Model Mechanical Behavior
At the microstructure level, the rock material is generally considered to be a combination of many discrete polygonal grains (Figure 1a). The Voronoi tessellation framework has been verified that naturally agrees with the granular meso-or micro-structure of rock materials [29,40,41]. It is mainly composed of the lattice element method and spring network methods [40][41][42]. UDEC has an in-built, automatic generator of the Voronoi tessellation pattern, where a particular region in a model can be subdivided into randomly sized polygons. This assemblage of distinct deformable polygons model is called the Voronoi model. The Voronoi model is made up of a set of lattice nodes that can be regularly and irregularly distributed along the domain [40,41]. Lattice nodes are connected with lattice elements, usually forming a Delaunay triangulation, and each triangle is circumscribed within a circle containing its three vertexes. Finally, the Voronoi polygons are created by constructing perpendicular bisections of all the triangles that share a common side [40,43] ( Figure 1a). the tensile strength of rock materials by adjusting the tensile strength of the block contact surface. It solves the problem of the ratio of unconfined compressive strength to tensile strength [28,34]. The Voronoi model has the advantage of using simple logical combinations to achieve close to real rock materials. In this study, the UDEC-Voronoi approach was used to build the numerical models for uniaxial/triaxial compressive and Brazilian split tests, and the microscopic parameters of contacts were obtained for four typical rock samples by calibration and inversion. Then, the uniaxial compressive tests of rock with four different compressive strengths were conducted. The characteristics of changes in the number, angle, and length of cracks in the process of rock failure and instability were monitored. The laws of crack development, damage evolution, and energy evolution were analyzed using the Voronoi model. The failure and instability mechanism of rock material under the discrete block model was discussed. It reveals the failure laws of rock materials at the micro-level and reproduces the whole process of crack evolution under natural conditions.

Voronoi Model Mechanical Behavior
At the microstructure level, the rock material is generally considered to be a combination of many discrete polygonal grains (Figure 1a). The Voronoi tessellation framework has been verified that naturally agrees with the granular meso-or micro-structure of rock materials [29,40,41]. It is mainly composed of the lattice element method and spring network methods [40][41][42]. UDEC has an in-built, automatic generator of the Voronoi tessellation pattern, where a particular region in a model can be subdivided into randomly sized polygons. This assemblage of distinct deformable polygons model is called the Voronoi model. The Voronoi model is made up of a set of lattice nodes that can be regularly and irregularly distributed along the domain [40,41]. Lattice nodes are connected with lattice elements, usually forming a Delaunay triangulation, and each triangle is circumscribed within a circle containing its three vertexes. Finally, the Voronoi polygons are created by constructing perpendicular bisections of all the triangles that share a common side [40,43]  The Voronoi polygon blocks are assigned an isotropic elastic deformable material model, which means that failure on the micro-scale only takes place at the block boundaries, not inside them. The mechanical properties of the connection between Voronoi polygon blocks obey the Coulomb friction law. Normal stiffness, shear stiffness, cohesion, friction, and tensile strength can be assigned to the contacts. These properties are referred to as the micro-properties. The Voronoi micro-mechanical behaviors and properties are shown in Figure 1 and Table 1.  [29,43]); (b) constitutive behavior of Voronoi model (modified from [44]) Note: The normal force (F n ) is controlled by the normal stiffness (k n ), and the shear force (F S ) is controlled by the shear stiffness (k s ); (c) numerical models of uniaxial/triaxial compressive and Brazilian split tests.
The Voronoi polygon blocks are assigned an isotropic elastic deformable material model, which means that failure on the micro-scale only takes place at the block boundaries, not inside them. The mechanical properties of the connection between Voronoi polygon blocks obey the Coulomb friction law. Normal stiffness, shear stiffness, cohesion, friction, and tensile strength can be assigned to the contacts. These properties are referred to as the micro-properties. The Voronoi micro-mechanical behaviors and properties are shown in Figure 1 and Table 1. In the normal direction, the contact force-displacement relation is assumed to be linear and governed by its normal stiffness (k n ) such that where ∆σ n is effective normal force increment and ∆u n is the normal displacement increment. A limiting tensile strength, σ max n , is assumed for any contact. If the tensile strength is exceeded (σ n ≤ −σ max n ), then σ n = 0, and the contact is marked as a tensile crack. In the opposite direction, blocks should be compressed and overlap at these contact points and the overlap is controlled by k n .
In the shear direction, the response is controlled by contact shear stiffness (k s ). The shear stress (τ s ) is limited by a combination of the cohesion C cont and friction ϕ cont .
or else, if |τ s |≥ τ max s (4) where a c is contact areas, ∆u e s is the elastic component of the incremental shear displacement and ∆u s is the total incremental shear displacement. The contact is marked as a shear crack if Equation (4) is achieved.
The Voronoi model's time step calculation is based on stiffness, and the connected nodes are generally given higher normal and shear stiffness values to prevent free movement along the joints. The normal stiffness and shear stiffness for the contacts in the numerical model can be derived from the equations in [25]: After one contact breaks, forces are redistributed, and it might cause adjacent contacts to break. Tension spring (k n ) and shear spring (k s ) work independently of each other. So, when the tension spring fails, the shear spring is still working, and vice versa. It can be further explained as follows: shear and tensile cracks can occur at the same contact at the same time, or only tensile cracks can occur. During the process, there is no need for a complex constitutive model to control cracking behavior. This implies that only contact properties are controlling material response on the macro-scale [34]. Figure 1b illustrates the contact normal and shear behavior.

Modelling Characteristics
According to different compressive strengths, Loc du Bonnet (Ldb) Granite, Augig Granite, Transjuane Sandstone (TS), and Coal, the four typical rock samples from hard rock to weak, were selected as the numerical simulation test samples of this paper. The detailed experimental parameters are shown in Table 2. According to the experimental sample standard recommended by ISRM [45], the numerical models of uniaxial/triaxial compression and Brazilian tension tests were established, as shown in Figure 1c. The specimen's size for the uniaxial/triaxial compression model is 50 mm × 100 mm (diameter × height), and the diameter of the sample for the Brazilian tension model is 50 mm. The Voronoi block is assigned as elasticity. The connection between the blocks follows the Coulomb slip model with residual strength properties [46]. The setting of the Voronoi block size also influences the failure mode [28,34,43]. Only when the block size is less than 1/10 of the minimum model sample size can the influence of the size be ignored [43,46]. Considering the calculation performance and time consumption, a too-small grain block size is also impractical, making the model calculation time-consuming. Referring to the actual grain size distribution of the selected samples, the average block size of LdB granite and Augig granite samples in the model is 4 mm [28,29,34,47], and the TS and Coal rock samples are 2 mm [28,34,48], which meets the size requirements. The loading rate is set to maintain the displacement rate at 0.01 m/s in both compression and tension tests [34,43]. Triaxial compression can be obtained by applying a constant lateral force to both sides of the uniaxial compression model. The stress and strain of the model during compression can be obtained through the displacement and force monitoring points as shown in Figure 1b. The numerical model for this paper took approximately 28 min to calculate (PC, Inter Core i9-9900k 3.6 GHz, Win10 64 bit). The model is in equilibrium and the calculation converges when the real-time maximum unbalance force is reduced to 1% of the initial maximum unbalance force.

Microscopic Mechanical Properties of Voronoi Model
In the Voronoi model, using the UDEC built-in FISH program to track and monitor the changes of the connection properties between blocks, the fracture evolution process of micro-cracks can be observed effectively. However, it is necessary to verify the validity of the model before numerical simulation. Therefore, we need to conduct a sensitivity analysis on the microscopic parameters of the model to prove that the microscopic mechanical properties of the numerical model can match the macroscopic properties of the rock. The micro-properties to be calibrated in the model are shown in Table 2. Previous studies [34,43,44] have analyzed the parameters of Voronoi to some extent. In this section, the classical LdB granite sample was used to conduct supplementary studies on the sensitivity of microscopic parameters by means of controlling variables, and feasible correction suggestions were proposed. Limited to space, only the representative simulation results are listed below. From Figure 2a, Young's modulus of the micro-block E b plays a leading role in Young's modulus E of the entire rock and the relationship between the two is increasing linearly, while E b will cause a slight decrease in Poisson's ratio, and has little effect on the strength. When the value of E b gradually increases, the polygons behave more stiffly. The cumulative deformation mainly focuses on contacts and ignores the elastic deformation of the block. Thus, the material's Young's modulus increases and Poisson's ratio decreases. When E b is assigned a lower value, the polygons are softer and the cumulative deformation mainly focuses on blocks, and thus the material's Young's modulus decreases and Poisson's ratio decreases. It is also shown that E b has a negligible influence on strength.
• Contact normal stiffness ( ) and shear stiffness ( ) It can be seen from Figure 2 that the Young's modulus E of the rock is also affected by the and between the microscopic joints. The increase in and ks in the early stage will cause E to increase exponentially, which far exceeds the influence of on E. However, for the Poisson's ratio v of the material, becomes an exponentially positive correlation and becomes an exponentially negative correlation. The Young's modulus (E) and Poisson's ratio (v) of the rock are mainly affected by and in the stiffness ratio between 0.1 and 1.1 (kn/ks is close to 0.9-10). When this ratio is exceeded, both E and v of the rock will tend to a constant value.
• Contact stiffness ratio (ks/kn) Some scholars have verified that the ratio of shear stiffness to normal stiffness determines Poisson's ratio of macro samples [28,34]. However, different stiffness ratios have different effects on the material's Poisson ratio. The smaller stiffness ratio will increase the elasticity of the rock, resulting in a larger Poisson's ratio. If the stiffness ratio taken is too small, the axial stress is much greater than the vertical stress. It will result in an internal force imbalance that cannot be loaded statically.

Micro Strength Properties
• Contact-Cohesion ( ) and Contact-Friction angle ( ) In Figure 2, the contact-cohesion and contact-friction angles do not affect the elastic • Contact normal stiffness (k n ) and shear stiffness (k s ) It can be seen from Figure 2 that the Young's modulus E of the rock is also affected by the k n and k s between the microscopic joints. The increase in k n and k s in the early stage will cause E to increase exponentially, which far exceeds the influence of E b on E. However, for the Poisson's ratio v of the material, k n becomes an exponentially positive correlation and k s becomes an exponentially negative correlation. The Young's modulus (E) and Poisson's ratio (v) of the rock are mainly affected by k n and k s in the stiffness ratio between 0.1 and 1.1 (k n /k s is close to 0.9-10). When this ratio is exceeded, both E and v of the rock will tend to a constant value.

•
Contact stiffness ratio (k s /k n ) Some scholars have verified that the ratio of shear stiffness to normal stiffness determines Poisson's ratio of macro samples [28,34]. However, different stiffness ratios have different effects on the material's Poisson ratio. The smaller stiffness ratio will increase the elasticity of the rock, resulting in a larger Poisson's ratio. If the stiffness ratio taken is too small, the axial stress is much greater than the vertical stress. It will result in an internal force imbalance that cannot be loaded statically.

•
Contact-Cohesion (c j ) and Contact-Friction angle (ϕ j ) In Figure 2, the contact-cohesion and contact-friction angles do not affect the elastic modes, Poisson's ratio, and tensile strength, which only have the linear incremental relationship of UCS. It is shown that contact-cohesion and contact-friction angle are the main factors affecting strength. To further explore the influence of contact-cohesion and contact-friction angle on the strength of intact rock, triaxial compression simulations are established. The confining pressure range is suggested in 0 ≤ σ 3 ≤ UCS/10 [50]. Therefore, the confining pressure value of this experiment is 0, 4, 8, 12 Mpa, respectively. The contactcohesion is assumed to be from 35 Mpa to 60 Mpa and increments at 5 Mpa intervals while keeping the other parameters constant. Variation of the contact-friction angle are between 10 • and 60 • with increments of 10 • . The friction angle (ϕ) and cohesion (c) of intact rock were obtained [27,44] . N ϕ was defined [27] by the peak strength (σ f ) and for confinements in the range P 0 (0 Mpa) − P 1 (12 Mpa) via Figure 3, the contact-cohesion is only linearly positively related to the cohesion of intact rock. The contact-friction angle not only plays a major role in the internal friction angle of the material but also has a greater influence on the cohesion of intact rock, especially when the contact-friction angle is larger. The explanation may be that the Voronoi model is a convex polygon block, and the friction angle between the blocks is too large, which will cause the interlocking effect to be formed, and also make the strength of the sample unrealistically large. Therefore, it is not desirable to take too large a value of the internal friction angle.
In the Brazilian tension test, the indirect tensile strength of a cylindrical sample is given by σ t = P max πRt . where P max is the maximum load at failure, R is the radius of the sample, and t is the thickness of the sample. In Figure 2g, the tensile strength of intact rock is only controlled by the micro tensile strength. Therefore, the Brazilian tension numerical model can be established to adjust the micro tensile parameters so that the tensile strength of the model reaches the experimental tensile strength. where is the maximum load at failure, R is the radius of the sample, and t is the thickness of the sample. In Figure 2g, the tensile strength of intact rock is only controlled by the micro tensile strength. Therefore, the Brazilian tension numerical model can be established to adjust the micro tensile parameters so that the tensile strength of the model reaches the experimental tensile strength.

Post-Peak Stress Parameters
When the peak strength exceeds, the connection bond between blocks will be broken to generate cracks. Currently, there is no residual cohesion ( ) and residual tensile strength ( ) in cracks, so the value is 0. The residual friction angle ( ) is an inherent property of the polygon block and still exists. In Figure 2f, the residual friction angle does not affect the deformation and strength before peak strength. The residual friction angle mainly plays a certain carrying role in the post-peak damage [51]. The carrying capacity is mainly provided by generating greater frictional resistance for the micro-block's dislocation motion. When the residual friction angle value is too large, it will also cause the rock sample to appear interlocked, which is a phenomenon and result of the post-peak enhancement effect. To meet the actual rock sample damage pressure, by trial and error, the value should not be greater than half of the friction angle.

Calibration Procedure and Results Analysis
The elastic parameters ( , , , ) and the strength parameters ( , , ) of microscopic blocks control the deformation ( , ) and strength ( , , ) of intact rocks. The initial parameters of the model were set according to the experimental data (i.e., = , = ). The block size was chosen according to mineralogy. The contact stiffness ratio (ks/kn) is determined to macro-Poisson's ratio. Once the contact stiffness ratio was set, both the

Post-Peak Stress Parameters
When the peak strength exceeds, the connection bond between blocks will be broken to generate cracks. Currently, there is no residual cohesion (c jr ) and residual tensile strength (σ tjr ) in cracks, so the value is 0. The residual friction angle (ϕ jr ) is an inherent property of the polygon block and still exists. In Figure 2f, the residual friction angle does not affect the deformation and strength before peak strength. The residual friction angle mainly plays a certain carrying role in the post-peak damage [51]. The carrying capacity is mainly provided by generating greater frictional resistance for the micro-block's dislocation motion. When the residual friction angle value is too large, it will also cause the rock sample to appear interlocked, which is a phenomenon and result of the post-peak enhancement effect. To meet the actual rock sample damage pressure, by trial and error, the value should not be greater than half of the friction angle.

Calibration Procedure and Results Analysis
The elastic parameters (E b ,v b ,k n ,k s ) and the strength parameters (c j ,ϕ j ,σ tj ) of microscopic blocks control the deformation (E,v) and strength (c,ϕ,σ t ) of intact rocks. The initial parameters of the model were set according to the experimental data (i.e., The block size was chosen according to mineralogy. The contact stiffness ratio (k s /k n ) is determined to macro-Poisson's ratio. Once the contact stiffness ratio was set, both the normal stiffness (k n ) and block deformability (E b ) were altered to fit the macro-Young's Modulus (E). The calibration procedure of the numerical model followed the procedure outlines by Christianson et al. [52], Kazerani and Zhao [34], and Gao and Stead [28]. By trial and error, a unique set of contact parameters, which satisfies the material properties, including Young's model, Poisson's ratio, tensile strength, friction angle, cohesion, and uniaxial compressive strength, were established (see Table 3).  Original cracks and pores exist in real rock samples, so the stress-strain curve in the experimental data has a crack closure phase in the early stages, which is not considered in the Voronoi model. When the rock sample is compressed, the stress-strain curve jumps over the crack closure and directly goes through the elastic, elastic-plastic, and plastic yield. Therefore, it is necessary to remove the nonlinear phase (crack fracturing stage) in the early stage of the experiment when correcting the model. The stress-strain curve of the modified model was moved to the right to the elastic part of the experimental stress-strain curve (the thumbnail in the upper left corner of the model), and the stress-strain curve of the corrected model is consistent with the experimental results, as shown in Figure 4. Due to the high degree of denseness of LdB rock samples, there are very few original cracks and pores, and the pre-pressing stage is not obvious. Therefore, the LdB rock numerical model can ignore the pre-nonlinear stage and correct the model directly. The error values of macro parameters obtained by the four simulation samples are all less than ±7%, which is close to the actual results (see Table 4). The tensile compression ratio is mainly between 0.05 and 0.08, which is consistent with the experimental results. Based on the above conclusions, the microscopic calibration parameters can accurately reflect the macroscopic mechanical properties of rock samples.

Micro-Crack Evolution in the Rock Fracture Process
Uniaxial compression experiments are extensively used to study the failure process of intact rock. Meanwhile, rock is a brittle material, and its failure process cannot be obtained by direct observation. To understand the evolution of fracture development, it is necessary to reconstruct the fracture development process by numerical simulation. According to the failure criterion of the Voronoi model principle, the contact surface between the two blocks could be monitored and judged. If the tensile strength/shear stress is exceeded, the contact will be recorded as a tensile crack/shear crack.

Strain-Stress Crack Analysis
Eberhardt et al. [53] conducted a series of unconfined compressive tests on brittle rock and found that there are two main thresholds before the peak strength: the crack initiation threshold ( ) and the crack damage threshold ( ). Both were identified through acoustic emission and strain monitoring. The crack initiation threshold is determined by the first appearance of the crack and the beginning of stable growth. The crack damage threshold is mainly determined by the unstable development of cracks, which is mainly manifested as the cracks begin to increase sharply. Both the crack initiation and damage thresholds were successfully captured in the Voronoi models. Cai and Hoek et al. [54,55] obtained the distribution range of the crack initiation threshold and crack damage threshold in the range of 30-60% and 70-90% of uniaxial compressive strength. Xue et al. [56], through the data statistics of the uniaxial compression experiment, found that the crack damage threshold of most sedimentary rocks is between 60-90%. In Figure 5, the simulation result shows that the crack initiation thresholds (Point A) of the four rock samples under the uniaxial compression model were 35.97%, 31.88%, 23.25%, and 23.83% of UCS, respectively. These values are slightly lower than 30-60% of the experimental statistical value of UCS, which is mainly due to the fact that the compaction stage in the early stage is not considered, which makes the initial threshold of cracks move forward. The

Micro-Crack Evolution in the Rock Fracture Process
Uniaxial compression experiments are extensively used to study the failure process of intact rock. Meanwhile, rock is a brittle material, and its failure process cannot be obtained by direct observation. To understand the evolution of fracture development, it is necessary to reconstruct the fracture development process by numerical simulation. According to the failure criterion of the Voronoi model principle, the contact surface between the two blocks could be monitored and judged. If the tensile strength/shear stress is exceeded, the contact will be recorded as a tensile crack/shear crack.

Strain-Stress Crack Analysis
Eberhardt et al. [53] conducted a series of unconfined compressive tests on brittle rock and found that there are two main thresholds before the peak strength: the crack initiation threshold (σ ci ) and the crack damage threshold (σ cd ). Both were identified through acoustic emission and strain monitoring. The crack initiation threshold is determined by the first appearance of the crack and the beginning of stable growth. The crack damage threshold is mainly determined by the unstable development of cracks, which is mainly manifested as the cracks begin to increase sharply. Both the crack initiation and damage thresholds were successfully captured in the Voronoi models. Cai and Hoek et al. [54,55] obtained the distribution range of the crack initiation threshold and crack damage threshold in the range of 30-60% and 70-90% of uniaxial compressive strength. Xue et al. [56], through the data statistics of the uniaxial compression experiment, found that the crack damage threshold of most sedimentary rocks is between 60-90%. In Figure 5, the simulation result shows that the crack initiation thresholds (Point A) of the four rock samples under the uniaxial compression model were 35.97%, 31.88%, 23.25%, and 23.83% of UCS, respectively. These values are slightly lower than 30-60% of the experimental statistical value of UCS, which is mainly due to the fact that the compaction stage in the early stage is not considered, which makes the initial threshold of cracks move forward. The crack damage threshold values (Point B) of four rock samples were 88.88%, 91.2%, 90.75% and 93.2% of UCS, respectively and the values are in line with the experimental statistics range. capacity. After exceeding the peak strength point (C-D stage), the shear crack and tensile crack increase sharply and lose the bearing capacity in a short time.
Comparing the numerical simulation results of the four rock samples, the number and increment of shear cracks in hard rock are significantly higher than those of tensile cracks, while the number and increment of tensile cracks in the weak rock are higher than those in the shear crack. It shows that shear cracks are the main form of hard rock failure, while tensile cracks are the main form of failure in weak rock. No matter the kinds of rock fractures, the generation, and evolution of shear cracks all occur near the peak strength. When the peak strength is exceeded, the number of shear cracks decreases rapidly and tends to a lower value. The development of the tensile crack runs through the whole failure process. Hoek and Bienawski [57] divided crack development into five stages: consolidation, elastic deformation, stable crack development, unstable crack development, and post-peak stage. Because the Voronoi model did not consider the original cracks or pores, there will be no crack closure in the simulated stress-strain curve. In this model, the development of cracks mainly includes the elastic deformation stage (0-A), stable crack development (A-B), unstable crack development (B-C), and post-peak stage (C-D); see Figure 5.
At the elastic deformation stage, the four rock samples are mainly deformed by pressure and do not produce cracks. From the A-B stage, cracks gradually began to germinate. This period is the initial development stage of crack and mainly focuses on tension cracks. A small number of shear cracks begin to initiate at this stage. When the B-C stage is reached, the damage of the rock sample develops rapidly. The growth rates of tension crack and shear crack increase. When the crack development reaches the peak strength (point C), the rock sample still keeps stability and has an effective bearing capacity. After exceeding the peak strength point (C-D stage), the shear crack and tensile crack increase sharply and lose the bearing capacity in a short time.
Comparing the numerical simulation results of the four rock samples, the number and increment of shear cracks in hard rock are significantly higher than those of tensile cracks, while the number and increment of tensile cracks in the weak rock are higher than those in the shear crack. It shows that shear cracks are the main form of hard rock failure, while tensile cracks are the main form of failure in weak rock. No matter the kinds of rock fractures, the generation, and evolution of shear cracks all occur near the peak strength. When the peak strength is exceeded, the number of shear cracks decreases rapidly and tends to a lower value. The development of the tensile crack runs through the whole failure process.

Crack Evolution Process and Crack Angle Analysis
To better display the fracture trend, the whole process of crack evolution is quantitatively described by selecting rock samples at the initial crack threshold point (σ ci ), crack damage threshold point (σ cd ), peak failure point (σ c )and post-peak residual point (σ r ), as shown in Figure 6. Due to limited space, only LdB granite and coal samples with the highest uniaxial compression strength (184.5 Mpa) and the lowest (13.1 Mpa) of the four samples were selected to characterize hard rock and weak rock. It can be seen from the crack evolution process of LdB granite and coal shown in Figure 6. and 70° (135° and 150°), respectively, and the tendency of crack accumulation began to appear in some areas of the numerical model.
When the peak failure point ( ) is reached, the distribution densities of tensile cracks and shear cracks increase significantly, and the density of the tensile cracks' distribution is higher than that of the shear cracks. Tension cracks and shear cracks inside the model began to accumulate based on the original cracks. The microcracks gradually formed fissures and further expanded into multiple fissures. The development of the fissures was mainly concentrated on the diagonal line of the rock sample. After the fissures formed from the diagonal corners, they gradually expanded to the center of the rock sample along both ends of the diagonal line.
When loaded to the post-peak residual point ( ), the tensile cracks mainly extend along the loading direction of 90°, and the shear cracks mainly extend along the direction of 45-70° (135-150°). The distribution density and distribution area of the cracks increase significantly and finally form many macroscopic penetrating fractures.
Whether it is hard rock or weak rock in uniaxial compression failure, the tensile failure caused by the accumulation of tensile cracks mainly occurs along the direction of loading (failure occurs in the direction of 90°), and the shear failure formed by the shear cracks mainly occurs in the directions of 45° and 70° (135° and 150°), which is consistent with the macroscopic fracture formed in the experiment.

AE Event Count
The cracks inside the rock material release elastic waves when they are generated and extended. The way to count the number of elastic waves sent out by cracks is called acoustic emission (AE) count event. For the Voronoi model, the AE event can be simulated by counting the number of fracture joints. Although the numerical simulation method When the stress reaches the initial crack threshold (σ ci ), the damage distribution of the cracks is discrete, and mainly tensile cracks. The cracking angle is mostly concentrated at 90 • , the cracks are small, and the distribution density is low, in which no shear cracks are generated.
As the load increases and reaches the crack damage threshold (σ cd ), the density of the tensile crack distribution gradually increases, and the developed cracks mainly occur at the tip of the polygonal block and grow steadily in the form of random distribution. At the same time, the shear cracks began to appear, and the cracking angles were mainly 45 • and 70 • (135 • and 150 • ), respectively, and the tendency of crack accumulation began to appear in some areas of the numerical model.
When the peak failure point (σ c ) is reached, the distribution densities of tensile cracks and shear cracks increase significantly, and the density of the tensile cracks' distribution is higher than that of the shear cracks. Tension cracks and shear cracks inside the model began to accumulate based on the original cracks. The microcracks gradually formed fissures and further expanded into multiple fissures. The development of the fissures was mainly concentrated on the diagonal line of the rock sample. After the fissures formed from the diagonal corners, they gradually expanded to the center of the rock sample along both ends of the diagonal line.
When loaded to the post-peak residual point (σ r ), the tensile cracks mainly extend along the loading direction of 90 • , and the shear cracks mainly extend along the direction of 45-70 • (135-150 • ). The distribution density and distribution area of the cracks increase significantly and finally form many macroscopic penetrating fractures.
Whether it is hard rock or weak rock in uniaxial compression failure, the tensile failure caused by the accumulation of tensile cracks mainly occurs along the direction of loading (failure occurs in the direction of 90 • ), and the shear failure formed by the shear cracks mainly occurs in the directions of 45 • and 70 • (135 • and 150 • ), which is consistent with the macroscopic fracture formed in the experiment.

AE Event Count
The cracks inside the rock material release elastic waves when they are generated and extended. The way to count the number of elastic waves sent out by cracks is called acoustic emission (AE) count event. For the Voronoi model, the AE event can be simulated by counting the number of fracture joints. Although the numerical simulation method cannot completely reproduce the same situation as the experiment, the statistical contact faults can also better reproduce the number of cracks changes in the rock to characterize the acoustic emission. This method of AE simulation has been well applied in previous studies [58][59][60].
As can be seen from Figure 7, the four rock samples' acoustic emission characteristics are roughly divided into three stages: the initial silence period, the intermediate stabilization period, and the last peak period.
Initial silence period (segment I): This period corresponds to the line elasticity phase of stress-strain, which is mainly the elastic deformation behavior that occurs in the pressure between blocks, so the acoustic emission signal is very weak or not even in this segment. The small amount of AE event that appears in TS sand and coal in this period is caused by uneven force due to the uneven block shape. [53,61], which further proves that the Voronoi model is effective in simulating the rock fracture process. Combined with Figures 5 and 7, it can be seen that the process of hard rock and weak rock crack generation, propagation and penetration into the macro-crack is very short. The acoustic emission signal is strong and concentrated, showing the brittle failure characteristics of the rock. Using the obvious characteristics of acoustic emission count segmentation, recording the failure precursor information recognition point can also effectively provide a safety warning for rock fracture.

Crack Damage Analysis
On the microscale, the formation of fractures is the result of microscopic crack generation and accumulation. Based on the theory of damage mechanics, the damage degree of intact rock can be defined as = × 100%; where is the total length of cracks that have broken at different times, is the total length of contacts without damage, and is the damage degree of intact rock in crack length. In Figure 8, when the crack damage value reaches 6-10%, the rock is close to the peak strength. When the damage degree exceeds 10%, the shear crack damage value increases greatly, all the rock samples Intermediate stabilization period (segment II): After the accumulation of energy during the line elastic phase, micro-cracks develop gradually inside the rock, and the AE event begins. As the loading increases, the number of AE events begins to increase steadily; this phase takes a long time.
Last peak period (segment III): When loading into the line elasticity to yield the weakening phase is excessive, the number of cracks increase and the AE event count begins to increase dramatically. At this time, before the peak will appear a convex peak value. When the pressure load exceeds the peak strength, the micro-crack begins to expand, evolve and gradually develop into a penetrating macro-crack. The AE event count drops sharply after reaching its maximum peak, at which point the rock sample has broken. The beginning of the last peak stage indicates that the rock has entered the yield failure stage, so the higher peak value of the bulge can be used as the information recognition point of the rock failure precursor (in Figure 7, Point E).
These three stages of AE counts have good consistency with the experimental test [53,61], which further proves that the Voronoi model is effective in simulating the rock fracture process. Combined with Figures 5 and 7, it can be seen that the process of hard rock and weak rock crack generation, propagation and penetration into the macro-crack is very short. The acoustic emission signal is strong and concentrated, showing the brittle failure characteristics of the rock. Using the obvious characteristics of acoustic emission count segmentation, recording the failure precursor information recognition point can also effectively provide a safety warning for rock fracture.

Crack Damage Analysis
On the microscale, the formation of fractures is the result of microscopic crack generation and accumulation. Based on the theory of damage mechanics, the damage degree of intact rock can be defined as L crack = L D L O × 100%; where L D is the total length of cracks that have broken at different times, L O is the total length of contacts without damage, and L crack is the damage degree of intact rock in crack length. In Figure 8, when the crack damage value reaches 6-10%, the rock is close to the peak strength. When the damage degree exceeds 10%, the shear crack damage value increases greatly, all the rock samples reach the peak strength and damage occurs. The rock gradually loses its effective bearing capacity. The large increase of shear crack length is an important signal of rock failure. When the crack damage rate of rock is more than 45%, most of the rock has lost its bearing capacity and all of them are destroyed.
When the crack damage rate of rock is more than 45%, most of the rock has lost its bearing capacity and all of them are destroyed. Based on the above analysis of crack damage rate, we can further calculate the damage model of the rock rupture process. The scholars Kachanov et al. [62] defined the damage variable as = . where D is the cumulative damage variable, is the area of the material damage cross-sectional at a certain period, and A is the cross-sectional material area without initial damage.
Suppose that is the total cumulative crack length of the rock sample from initial no damage to complete loss of bearing capacity. The cumulative crack length per unit area damage is , = . where the cross-sectional area of damage reaches is , and the cumulative crack length is , = = .
Thus, the damage variable can be converted to = .
(a) (b) (c) (d) It is difficult for rock samples to achieve an absolute complete failure mode during compression. Therefore, the total crack length of rock samples under axial residual stress is regarded as the cumulative total crack length of rock samples at the complete failure time. From the previous analysis, it can be seen that the rock will completely lose its bearing capacity when the crack damage rate exceeds 45%. Therefore, the crack length with the crack damage rate of 50% is used as the total cumulative crack length.
Based on the crack length and the strain equivalence principle [63], a damage constitutive model of hard to weak rock specimens under uniaxial compression is established: Based on the above analysis of crack damage rate, we can further calculate the damage model of the rock rupture process. The scholars Kachanov et al. [62] defined the damage variable as where D is the cumulative damage variable, A d is the area of the material damage cross-sectional at a certain period, and A is the cross-sectional material area without initial damage.
Suppose that C 0 is the total cumulative crack length of the rock sample from initial no damage to complete loss of bearing capacity. The cumulative crack length per unit area damage is C w , C w = C 0 A . where the cross-sectional area of damage reaches is A d , and the cumulative crack length is C d , Thus, the damage variable can be converted to D = C d C 0 . It is difficult for rock samples to achieve an absolute complete failure mode during compression. Therefore, the total crack length of rock samples under axial residual stress is regarded as the cumulative total crack length of rock samples at the complete failure time. From the previous analysis, it can be seen that the rock will completely lose its bearing capacity when the crack damage rate exceeds 45%. Therefore, the crack length with the crack damage rate of 50% is used as the total cumulative crack length.
Based on the crack length and the strain equivalence principle [63], a damage constitutive model of hard to weak rock specimens under uniaxial compression is established: Figure 9 shows the stress-strain curve fitted by the damage constitutive equation of rock samples based on crack length. The fitting curve is in agreement with the numerical simulation curve. In summary, it is feasible to use the crack length parameter to reflect the damage evolution characteristics of rock samples.  Figure 9 shows the stress-strain curve fitted by the damage constitutive equation of rock samples based on crack length. The fitting curve is in agreement with the numerical simulation curve. In summary, it is feasible to use the crack length parameter to reflect the damage evolution characteristics of rock samples.

Characteristics of Macroscopic Failure
With the increasing load, the microscopic tension and shear cracks in the rock accumulate continuously and eventually form macro-fractures. According to the trend of the fracture and the relative movement directions of the blocks on both sides, it can be divided into the tensile fracture and shear fracture. If the fracture angle was vertical or sub-vertical (i.e., parallel to the loading direction) and the blocks on both sides of the crack move in the opposite direction (i.e., opening), it is defined as a tensile fracture. If the fracture angle was oblique and the blocks on both sides of the crack are sliding, it is defined as a shear fracture. Combined with the horizontal displacement contour plot, the macroscopic failure trend curves of tensile fracture and shear fracture can be obtained, as shown in Figure 10.

Characteristics of Macroscopic Failure
With the increasing load, the microscopic tension and shear cracks in the rock accumulate continuously and eventually form macro-fractures. According to the trend of the fracture and the relative movement directions of the blocks on both sides, it can be divided into the tensile fracture and shear fracture. If the fracture angle was vertical or sub-vertical (i.e., parallel to the loading direction) and the blocks on both sides of the crack move in the opposite direction (i.e., opening), it is defined as a tensile fracture. If the fracture angle was oblique and the blocks on both sides of the crack are sliding, it is defined as a shear fracture. Combined with the horizontal displacement contour plot, the macroscopic failure trend curves of tensile fracture and shear fracture can be obtained, as shown in Figure 10.

Energy Evolution Pattern
The principle of an AE event is mainly to record the energy information changes in a certain frequency band, but it cannot accurately detect the fracture information of weak joints. In the Voronoi model, the fracture between the contact surfaces of blocks will release the strain energy. The energy evolution characteristics of each crack can be obtained by recording the dissipation information of strain energy when cracks are generated. For each contact, the energy accumulation is governed using Equations (8) and (9): where / is the energy accumulation in tension/shear crack, / and / are the normal/shear stress of contact for current and previous time steps, respectively, and / is the normal/shear displacement of contact for the current time step. At the macro-scale, Figure 11a shows the whole process of energy accumulation and dissipation of four rock samples from hard rock to weak rock. Hard rock gathers more energy before failure than weak rock does, but the energy dissipates faster after failure than weak rock.

Energy Evolution Pattern
The principle of an AE event is mainly to record the energy information changes in a certain frequency band, but it cannot accurately detect the fracture information of weak joints. In the Voronoi model, the fracture between the contact surfaces of blocks will release the strain energy. The energy evolution characteristics of each crack can be obtained by recording the dissipation information of strain energy when cracks are generated. For each contact, the energy accumulation is governed using Equations (8) and (9): where U jct /U jcs is the energy accumulation in tension/shear crack, f n / f s and f n / f s are the normal/shear stress of contact for current and previous time steps, respectively, and u n /u s is the normal/shear displacement of contact for the current time step. At the macro-scale, Figure 11a shows the whole process of energy accumulation and dissipation of four rock samples from hard rock to weak rock. Hard rock gathers more energy before failure than weak rock does, but the energy dissipates faster after failure than weak rock. At the micro-scale, see Figures 11b and 12. Although the energy of tensile cracks keeps gathering, it is far less than that of shear cracks. The changing trend of total energy dissipation is consistent with that of shear cracks. The energy change caused by shear cracks is the main control factor of total energy dissipation. When the tensile crack is broken, the required energy is small and the energy is not concentrated. The damage caused by tensile cracks is not obvious. However, in the energy evolution process of shear cracks, energy dissipation is concentrated, and the macroscopic fractures formed are apparent. This indicates that the shear cracks cause strong damage to the rock. At the micro-scale, see Figure 11b or Figure 12. Although the energy of tensile cracks keeps gathering, it is far less than that of shear cracks. The changing trend of total energy dissipation is consistent with that of shear cracks. The energy change caused by shear cracks is the main control factor of total energy dissipation. When the tensile crack is broken, the required energy is small and the energy is not concentrated. The damage caused by tensile cracks is not obvious. However, in the energy evolution process of shear cracks, energy dissipation is concentrated, and the macroscopic fractures formed are apparent. This indicates that the shear cracks cause strong damage to the rock. At the micro-scale, see Figures 11b and 12. Although the energy of tensile cracks keeps gathering, it is far less than that of shear cracks. The changing trend of total energy dissipation is consistent with that of shear cracks. The energy change caused by shear cracks is the main control factor of total energy dissipation. When the tensile crack is broken, the required energy is small and the energy is not concentrated. The damage caused by tensile cracks is not obvious. However, in the energy evolution process of shear cracks, energy dissipation is concentrated, and the macroscopic fractures formed are apparent. This indicates that the shear cracks cause strong damage to the rock.

Discussion and Conclusions
This paper uses the Voronoi discrete block model to reproduce the whole process of the crack rupture effectively. The microscopic connection parameters between the blocks are the basis for the establishment of the Voronoi model. Many numerical simulations were used to analyze the sensitivity of the microscopic parameters. Therefore, the relationship between the microscopic block's mechanical properties and the macroscopic sample was obtained through sensitivity analysis. Then, the micro-parameters of the material were inverted through the stress-strain curve. Finally, the process of calibrating the model was quantified, and a reliable numerical model was obtained. If the rock sample has a high degree of compactness, the numerical model can ignore the previous non-linear stage and directly correct it.
(1) The crack development process of the four simulated samples is consistent with the experimental results, and all have experienced elastic deformation, stable crack development, unstable crack development, and post-peak failure. The period between the damage threshold and the peak strength is the most critical period for crack development. Before reaching the damage threshold, it is mainly the initiation and stable development of tensile cracks and it is accompanied by a few shear cracks. The appearance of tensile cracks earlier than shear cracks proves that the tensile resistance of rocks is much weaker than the shear resistance, and that the initiation of internal defects is mainly caused by tensile failure. When the damage threshold is exceeded and close to the peak strength, the internal damage rises sharply and develops into tensile-shear mixed damage. Although the macroscopic failure displacement is not obvious, a large amount of damage occurs inside the rock, reaching the limit of bearing capacity.
(2) Tensile cracks first appeared and mainly developed in the tip part of rock blocks, with a cracking angle of 90 • . Tensile cracks formed along the loading direction and developed throughout the whole fracture process. The development of shear cracks is mainly concentrated near the peak strength, with 45 • and 70 • (135 • and 150 • ) as the main cracking angle, which is the main indicator of the rock damage threshold. The damage caused by shear cracks is the main form in hard rock, while tensile cracks are the main form of damage in weak rock. However, no matter what kind of rock is damaged, the process of crack generation and expansion until it becomes a macro fracture is very short, showing the characteristics of brittle fracture of the rock. The development of micro-cracks is mainly concentrated on the diagonal of the rock sample and gradually expands to the middle along the two ends of the diagonal. Finally, the penetrating fractures along the vertical and horizontal dip angles are formed.
(3) The acoustic emission characteristics of the numerical model can be roughly divided into three stages: the initial silent period, the intermediate stable period, and the last peak period, which is consistent with the acoustic emission characteristics of the experimental statistics and can effectively capture the crack changes. The maximum number of acoustic emission (AE) events appears slightly earlier than the peak strength, and the identification point of failure precursor information can effectively provide a safety warning for the development of rock fracture.
(4) When the crack length damage rate reaches 10%, the rock reaches its peak strength, and when it exceeds 45%, the rock is destroyed. The uniaxial compression damage constitutive equation of the rock sample with the crack length as the parameter is established, which can better reflect the damage evolution characteristics of the rock sample.
(5) In the process of rock rupture, tensile failure requires low energy consumption, the energy dispersion is not concentrated, and the damage is not obvious; while the derived energy of shear cracks is concentrated and destructive, macroscopic cracks are easily formed.  Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.