# Effect of the Number and Distribution of Circular Holes on the Damage and Fracture of Limestone Using PFC2D

Authors: Handong Liu; Shuai Liu; Haining Liu; Yunyou Fu; Minglei Kang; Zhe Bai; Jianglong Xu; Chao Zheng
Journal: Geofluids
 (2022)

---

## Abstract

The common occurrence of various number and distribution of hole flaws complicates the mechanical behavior and fracture mode of the rock masses that contain them. This study develops seven numerical models of limestone samples with different numbers and distributions of circular hole flaws using 2-D particle flow discrete element code (PFC2D) to investigate their impact on the mechanical properties of limestone while maintaining the same flaws area. In addition, it analyzes the effect of these factors on the mesomechanical features of each model, including the characteristic stress values (peak stress, crack initiation, and damage stress), crack evolution, stress, and displacement field. The results showed that the peak stress, crack initiation, and damage stress of the single-hole model are between those of multihole models. As the arrangement of dip angle increases, the peak stress, crack initiation, and damage stress of models with the same number of multihole flaws exhibit a V-shaped change. The characteristic stress values are the largest when the holes are vertically aligned. Model differences in crack development path, shape, and number, as well as stress concentration area and failure mode, are primarily due to the number and distribution of holes. The circular holes are arranged at approximately 45°, and the greater the number of defects, the more likely the model is to fail. The study’s findings can provide support and reference for the research system of deformation and failure of rock mass with hole flaws.


---

## Body

## 1. Introduction

Under different geological processes, naturally occurring rock masses often contain various holes, fissures, and other flaws [1, 2]. For example, in limestone areas, small karst caves with diverse distribution patterns are common [3, 4]. These flaws are often numerous and randomly distributed, making the physical properties, deformation, and failure behavior of the rock mass that contain them complex and dynamic. This seriously affects the stability and security of rock masses during construction [5–7]. Consequently, it is necessary to investigate the effect of the number and distribution of hole flaws in limestones on their damage-deformation-failure behavior.Researchers have obtained many positive results in studying rock masses containing holes. The variable factors related to hole flaws are hole shape [8], hole size [9], hole number [10, 11], and their combination with other flaw types [12, 13]. Zeng et al. [14], Gui et al. [15], and Xia et al. [16] applied the discrete element numerical simulation method to analyze the influence of a single-hole flaw with multiple regular forms (circular, square, and trapezoidal) on the mesomechanical properties and failure mode of a rock mass. Li et al. [17] utilized digital image correlation to record and discuss the deformation and fracture of marble containing circular and elliptical holes. Zhao et al. [18] proposed a failure strength model for brittle materials with single circular hole flaws by improving the Sammis-Ashby model [19] and using a numerical model to test and verify the modified model. Similarly, Lin et al. [20] analyzed the failure mode of multiple-holed granite subjected to uniaxial loading and developed the Sammis-Ashby model. Yang et al. [21] and Tian et al. [22] examined the different distribution patterns of two elliptical holes based on laboratory and PFC2D simulation tests; they investigated the effect of hole spacing and hole arrangement angle on the damage evolution and failure mechanism of sandstone. Du et al. [23], Zhao et al. [24], and Ma et al. [25] conducted uniaxial compressive strength (UCS) tests on rock samples with two holes and studied the strength and crack propagation features of samples under various hole spacing configurations and different dip angles.For comparative analysis of the different number of hole flaws, it is common for most researchers to list one to three holes. For example, He et al. [26] compared red sandstone’s energy evolution and mechanical behavior with one to three circular holes under loading and unloading states. Tang et al. [27] conducted experiments on rock material with single and three-round holes and simulated the splitting failure observed in the test samples. Liu et al. [28] used PFC2D to reveal the damage evolution and distinguishing strength features of a composite rock mass when the weak interlayer contained one to three square holes and conducted research on the filling and reinforcement of the holes. Wong et al. and Wong and Lin [29, 30] conducted an array of indoor tests on rock samples with one, three, and multiple circular holes and investigated the failure mode and distinguishing strength features of rock samples with various distributions of circular holes. Gui et al. [31] developed single circular and multicircular hole numerical models with staggered distributions of various circular hole aperture sizes and numbers. They systematically analyzed the effect of the number and distribution of hole flaws on the mechanical behavior of the rock containing them. Huang et al. [32] drilled three noncollinear holes in granite samples and analyzed the rock samples’ crack development and failure regularity with a specific arrangement and combination of multiple small circular holes. Jespersen et al. [33] conducted laboratory tests and numerical simulations on gypsum specimens with seven uniformly distributed circular holes and discussed the effect of large hole spacing and arrangement on the strength and failure mode of the samples. Lin et al. [20] and Huang et al. [34] performed UCS tests on granite rock samples, including multiple holes, and focused on the impact of the holes on crack initiation, propagation, and coalescence. Using AE technology combined with PFC2D numerical simulation, the failure mechanism of the rock samples was revealed.The studies above have achieved fruitful results regarding the mechanical properties of lithosome containing hole flaws and introduced scientific and reliable research approaches such as the PFC2D method. Existing research on the quantity and distribution of hole flaws focuses primarily on holes of the same size. For example, in previous studies, the size of the single versus multiple-hole flaws in the rock sample was the same, or the influence of a fixed number of holes on the rock’s strength properties was investigated. However, increasing the number of holes will inevitably raise the total defect area in the rock mass, and rock samples will undoubtedly deteriorate. With the increase in the number of holes, the overall distribution of holes in the rock mass will inevitably change, and the distribution and number of holes are organic and unified research factors. The above literature indicates that few studies have been conducted on the influence of the number and distribution of holes on the fracture and distinguishing failure features of rock samples with the same flaw area. In addition, there is little research on the distinguishing mesomechanical features of rock, such as damage stress, based on the distribution and number of holes. Under the same total area of hole flaws, the design and analysis of different numbers and distributions of holes can more accurately represent the scientific nature of the research on their influence on the mechanical properties of a rock mass. Under the assumption of equal hole flaw area, this paper numerically simulates the mesomechanical properties of limestone samples with various numbers and distributions of circular holes using PFC2D code to reveal the damage and deformation mechanism of a rock mass.


## 2. Mesoparameter Calibration and Numerical Models

### 2.1. Determining of the Macro-Mesoparameters of Limestone

The parallel bond model in the PFC2D code was applied to simulate the rock material numerically, because it can simultaneously transfer force and moment and has an extensive contact range [35–37]. Using the parallel bond model will result in issues such as low compressive-tensile strength ratio and internal friction angle. However, modifying the parallel bond failure criteria to reduce the moment contribution of bond stress has solved the above problems better. The particle flow discrete element numerical simulation method has become exceedingly reliable and scientific in international rock mechanics [38, 39]. In order to determine limestone’s macro and mesoscopic parameters, its macroscopic physical and mechanical parameters must first be determined through laboratory tests. The laboratory tests in this study were conducted according to the ISRM specification [40]. In the UCS test, cylindrical rock samples with standard dimensions of 50 mm in diameter and 100 mm in height were used, and the axial deformation rate was set for loading to 0.001 mm/s. After obtaining the limestone’s basic physical and mechanical parameters, an isologous model was developed for the numerical simulation test and comparison. The mesoscopic parameters in the PFC2D code were repeatedly adjusted until the results of the experimental and numerical tests were perfectly consistent [41, 42]. Figure 1 depicts the limestone rock sample and its typical uniaxial compression numerical model and parallel bond contact in the PFC2D code.Figure 1
Limestone rock sample and corresponding intact numerical model.The UCS test was performed on the established numerical model of the intact rock sample using the identical method and loading conditions as the physical experiment. Figure2 illustrates the results of comparing the failure forms and uniaxial stress-strain curves of the rock models in the experimental and numerical tests. The peak strength and deformation modulus of the stress-strain curves of the experimental and numerical samples are identical, and the failure modes of the two are similar. This demonstrates that the microscopic parameters in the PFC model can accurately simulate the mechanical properties of macro rock samples [43]. The microscopic parameters are obtained as shown in Table 1.Figure 2
Verification of the experimental and numerical results.Table 1
Microscopic parameters of the rock.
ParametersRock (value)Minimum particle diameter (mm)0.2Particle diameter ratio1.5Grain density (g/cm3)3.05Contact modulus of the particle (GPa)36.21Contact bond gap (mm)0.05Porosity0.1Parallel bond friction angle (°)42.5Parallel bond tensile strength (MPa)10.95Normal critical damping ratio0.5Parallel bond cohesive force (MPa)52.23


### 2.2. Establish Numerical Models

The area of the holes in the established models remained unchanged, while the holes’ number and location were altered to illustrate better the effect of the number and distribution of holes on the mechanical features of the rock samples. Circular holes are a representative flaw shape in rock masses, so the hole shape in the models in this study was set as circular. The aperture diameter for a single circular hole is 15 mm (the hole area is about 1.77 cm2), so the aperture diameters for two and three circular holes are approximately 10.61 mm and 8.66 mm, respectively. In the models, the arrangement and distribution of two and three circular holes in the horizontal, 45°, and vertical directions are considered, along with the setting of equal distances between the hole boundaries and the circular holes and the model boundary. Thus, seven numerical models, designated NM-1 through NM-7, were established for the rock samples with circular hole flaws. Figure 3 depicts the numerical models. The NM-1 model contains a single-hole flaw with a diameter of 15 mm in the model’s center. NM-2 includes two circular holes with a horizontal distribution, and the distance between the two holes and the left and right boundaries of the model is equal (about 9.60 mm). NM-3 contains two circular holes arranged at 45° in the middle, and the distances between the two circular holes and the left and right boundaries of the model are about 16.50 mm. Two circular holes are arranged vertically in the center of the NM-4 model, and the distances between the two circular holes and the upper and lower boundaries of the model are about 26.26 mm. Similar to NM-2 to NM-4, NM-5 to NM-7 contains three circular holes with a diameter of 8.66 mm, and the distances in the horizontal, 45°, and vertical directions were 6.01, 11.18, and 18.51 mm, respectively.Figure 3
Numerical simulation models: (a–g) represent the numerical models NM-1 to NM-7, respectively.
(a)(b)(c)(d)(e)(f)(g)


### 2.3. Method for Monitoring the Cracks and Simulating the Acoustic Emission

The cohesion of particles is determined in the PFC2D code by the normal bond strength (n_bond), tangential bond strength (s_bond), and friction coefficient between particles [44–46]. When the transmitted stress exceeds the bond strength between rock particles, the contact bond between particles is destroyed under the load, and microcracks appear in the numerical model [47–49]. As the microcracks continue to expand in the model, the deformation energy will be released, forming acoustic waves. In this study, the AE events were simulated based on the change in the number of cracks generated in the PFC2D code, i.e., the AE hits were generated during the breaking of the contact bonds [50, 51].During the damage evolution of rock samples, crack initiation and damage stress are significant boundary values representing the rock’s strength. This study utilized the method of monitoring crack development and evolution to represent the damage stress and crack initiation stress, and the two define the boundary points of different periods in the crack growth process. The crack initiation stress refers to the original cracks in the compacted rock, after which the latent cracks begin to sprout and grow. Meanwhile, a crack was just initiated in the model, and the stress value corresponding to the first crack’s initiation was monitored [28, 52]. The damage stress refers to the convergence and connection of the initial and expanded microcracks in the previous stage, resulting in an unstable propagation of cracks. From this time until the peak stress in the model, the number of cracks shows a continuously increasing acceleration trend, and the corresponding AE events also begin to have a continuous and greater stable activity [53]. Previous studies [54–56] have revealed that the damage stress and crack initiation stress determined by the method of monitoring microcracks are identical to the results of laboratory tests. Although the simulated AE events cannot achieve the mechanical response state of a real rock mass, the reflected change aids in understanding the AE phenomenon of the rock.


## 2.1. Determining of the Macro-Mesoparameters of Limestone

The parallel bond model in the PFC2D code was applied to simulate the rock material numerically, because it can simultaneously transfer force and moment and has an extensive contact range [35–37]. Using the parallel bond model will result in issues such as low compressive-tensile strength ratio and internal friction angle. However, modifying the parallel bond failure criteria to reduce the moment contribution of bond stress has solved the above problems better. The particle flow discrete element numerical simulation method has become exceedingly reliable and scientific in international rock mechanics [38, 39]. In order to determine limestone’s macro and mesoscopic parameters, its macroscopic physical and mechanical parameters must first be determined through laboratory tests. The laboratory tests in this study were conducted according to the ISRM specification [40]. In the UCS test, cylindrical rock samples with standard dimensions of 50 mm in diameter and 100 mm in height were used, and the axial deformation rate was set for loading to 0.001 mm/s. After obtaining the limestone’s basic physical and mechanical parameters, an isologous model was developed for the numerical simulation test and comparison. The mesoscopic parameters in the PFC2D code were repeatedly adjusted until the results of the experimental and numerical tests were perfectly consistent [41, 42]. Figure 1 depicts the limestone rock sample and its typical uniaxial compression numerical model and parallel bond contact in the PFC2D code.Figure 1
Limestone rock sample and corresponding intact numerical model.The UCS test was performed on the established numerical model of the intact rock sample using the identical method and loading conditions as the physical experiment. Figure2 illustrates the results of comparing the failure forms and uniaxial stress-strain curves of the rock models in the experimental and numerical tests. The peak strength and deformation modulus of the stress-strain curves of the experimental and numerical samples are identical, and the failure modes of the two are similar. This demonstrates that the microscopic parameters in the PFC model can accurately simulate the mechanical properties of macro rock samples [43]. The microscopic parameters are obtained as shown in Table 1.Figure 2
Verification of the experimental and numerical results.Table 1
Microscopic parameters of the rock.
ParametersRock (value)Minimum particle diameter (mm)0.2Particle diameter ratio1.5Grain density (g/cm3)3.05Contact modulus of the particle (GPa)36.21Contact bond gap (mm)0.05Porosity0.1Parallel bond friction angle (°)42.5Parallel bond tensile strength (MPa)10.95Normal critical damping ratio0.5Parallel bond cohesive force (MPa)52.23


## 2.2. Establish Numerical Models

The area of the holes in the established models remained unchanged, while the holes’ number and location were altered to illustrate better the effect of the number and distribution of holes on the mechanical features of the rock samples. Circular holes are a representative flaw shape in rock masses, so the hole shape in the models in this study was set as circular. The aperture diameter for a single circular hole is 15 mm (the hole area is about 1.77 cm2), so the aperture diameters for two and three circular holes are approximately 10.61 mm and 8.66 mm, respectively. In the models, the arrangement and distribution of two and three circular holes in the horizontal, 45°, and vertical directions are considered, along with the setting of equal distances between the hole boundaries and the circular holes and the model boundary. Thus, seven numerical models, designated NM-1 through NM-7, were established for the rock samples with circular hole flaws. Figure 3 depicts the numerical models. The NM-1 model contains a single-hole flaw with a diameter of 15 mm in the model’s center. NM-2 includes two circular holes with a horizontal distribution, and the distance between the two holes and the left and right boundaries of the model is equal (about 9.60 mm). NM-3 contains two circular holes arranged at 45° in the middle, and the distances between the two circular holes and the left and right boundaries of the model are about 16.50 mm. Two circular holes are arranged vertically in the center of the NM-4 model, and the distances between the two circular holes and the upper and lower boundaries of the model are about 26.26 mm. Similar to NM-2 to NM-4, NM-5 to NM-7 contains three circular holes with a diameter of 8.66 mm, and the distances in the horizontal, 45°, and vertical directions were 6.01, 11.18, and 18.51 mm, respectively.Figure 3
Numerical simulation models: (a–g) represent the numerical models NM-1 to NM-7, respectively.
(a)(b)(c)(d)(e)(f)(g)


## 2.3. Method for Monitoring the Cracks and Simulating the Acoustic Emission

The cohesion of particles is determined in the PFC2D code by the normal bond strength (n_bond), tangential bond strength (s_bond), and friction coefficient between particles [44–46]. When the transmitted stress exceeds the bond strength between rock particles, the contact bond between particles is destroyed under the load, and microcracks appear in the numerical model [47–49]. As the microcracks continue to expand in the model, the deformation energy will be released, forming acoustic waves. In this study, the AE events were simulated based on the change in the number of cracks generated in the PFC2D code, i.e., the AE hits were generated during the breaking of the contact bonds [50, 51].During the damage evolution of rock samples, crack initiation and damage stress are significant boundary values representing the rock’s strength. This study utilized the method of monitoring crack development and evolution to represent the damage stress and crack initiation stress, and the two define the boundary points of different periods in the crack growth process. The crack initiation stress refers to the original cracks in the compacted rock, after which the latent cracks begin to sprout and grow. Meanwhile, a crack was just initiated in the model, and the stress value corresponding to the first crack’s initiation was monitored [28, 52]. The damage stress refers to the convergence and connection of the initial and expanded microcracks in the previous stage, resulting in an unstable propagation of cracks. From this time until the peak stress in the model, the number of cracks shows a continuously increasing acceleration trend, and the corresponding AE events also begin to have a continuous and greater stable activity [53]. Previous studies [54–56] have revealed that the damage stress and crack initiation stress determined by the method of monitoring microcracks are identical to the results of laboratory tests. Although the simulated AE events cannot achieve the mechanical response state of a real rock mass, the reflected change aids in understanding the AE phenomenon of the rock.


## 3. Analysis of Simulation Test Results

### 3.1. Analysis of AE Events and Characteristic Stress Values of each Model

The stress-strain, AE, and the totality of cracks for each model are summarized together, and the three characteristic stress values of crack initiation stress, damage stress, and peak stress of each model are marked on the stress-strain curve. The stress-strain-acoustic emission-crack number curves of each model are obtained, as shown in Figure4. The crack initiation, damage, and peak stress are denoted by σci, σcd, and σc, respectively. Table 2 lists the values of the crack initiation, damage, and peak strain-stress of the models.Figure 4
Stress-strain-acoustic emission-crack number curves of each model: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Table 2
Crack initiation, damage, and peak stress and strain of each model.
Model nameCrack initiation strain (%)Crack initiation stress (MPa)Damage strain (%)Damage stress (MPa)Peak strain (%)Peak stress (MPa)Intact model0.062548.160.1402106.520.1727125.05NM-10.019213.500.093461.430.149884.68NM-20.022114.990.084155.590.118570.15NM-30.015610.540.081852.870.137670.51NM-40.022715.670.103468.440.158392.28NM-50.025717.230.080252.750.113169.52NM-60.017011.350.080549.960.115263.86NM-70.029220.270.109573.300.150194.35The curves in Figure4 reveal that the UCS of the intact model is about 125.05 MPa, while the uniaxial compressive strengths of NM-1 to NM-7 are about 84.68, 70.15, 70.51, 92.28, 69.52, 63.86, and 94.35 MPa, respectively. The strength of the intact model is the largest, and the UCS of NM-1 to NM-7 is significantly reduced. Given the numbers and distribution of hole flaws, the UCS of the models containing the circular hole flaws is different, and the UCS of the two models with the vertical arrangement of circular holes is high. The strength of NM-2 and NM-3 is equivalent, while the strength of NM-5 is greater than that of NM-6. This indicates that when the two circular holes are arranged at 45°, the influence of the related flaw types, such as fissures and soft structural planes, on the rock mass is different [57, 58]. The main reason for this is that although the two circular holes are arranged at a 45° dip angle, the spacing between the circular holes is still relatively large, and penetration under uniaxial compression is insufficient. When the three circular holes are arranged at 45°, the distance between them decreases, and the effect of promoting the failure of the rock samples increases, resulting in a reduction in its bearing capacity; therefore, the NM-6 model has the smallest UCS of all the models. In general, the strength of the two-hole model with the horizontal and 45° distribution of holes is stronger than that of the corresponding three-hole model. However, when the circular hole flaws are arranged in the vertical direction, the UCS of the model with the three vertical holes is the largest, and the intact model is reduced by 24.55%. The strength of the model with the two holes is the second smallest, and the UCS of the intact model decreases by 26.21%. The three-hole model with the 45° inclination angle has the smallest UCS, which is 48.93% lower than the intact model.In the AE characteristic histograms of each model, the AE events do not occur when the initial loading stress is minor. As the loading stress increases, cracks appear within the rock sample, and acoustic emission phenomena appear in each model. The crack initiation stress of each model is about 48.16, 13.50, 14.99, 10.54, 15.67, 17.23, 11.35, and 20.27 MPa, respectively. The hole flaws significantly reduce the crack initiation stress of the models and cause the AE events to occur earlier. Two-hole and three-hole models with the 45° arrangement begin to exhibit AE activity at strains of 0.0156% and 0.0170%. Internal cracks in the rock samples gradually develop and expand with continuous loading, and the number of cracks maintains a stable growth rate. The number of AE counts of each model also continues to increase slowly. When loading to 70%-80% of the peak stress, the totality of cracks in each model gradually begins to increase rapidly, and the number of AE also counts suddenly increases. At this time, the model is considered to have reached the damage state, and the corresponding stress level is the damage stress. The damage stress of each model is about 106.52, 61.43, 55.59, 52.87, 68.44, 52.75, 49.96, and 73.30 MPa, respectively. The damage stress of the models with circular hole flaws is significantly reduced compared to the intact model, while the damage strain is also significantly increased. In terms of the damage changes and the crack initiation stress, the variations under various distributions of the two and three circular holes are similar, with the smallest changes occurring when the holes are arranged at 45°. Accordingly, the damage degree is greater under the same loading condition, and the values are the largest when the holes are distributed vertically, resulting in a lower damage degree. After the damage stress of each model is reached, the total number of cracks in each model increases rapidly, and the number of acoustic counts rises sharply at the peak stress, where it reaches the maximum. The acoustic emission is “active” around the peak stress until the postpeak failure stage, and the AE events are concentrated around the peak stress. The number of cracks in each model increases rapidly again, but the total number of AE counts decreases in the postpeak stage. Overall, the intact model has the highest number of AE counts and the most concentrated AE event distribution. Each model’s maximum number of AE counts is small, and the AE events are scattered around the peak stress. The distribution and number of the circular hole flaws also influence the distinguishing AE features. They affect the early and late occurrence of AE events at the crack initiation stress level and the number of acoustic emission shots at the peak stress. The changing trend is equivalent to each model’s peak stress, and the model’s AE events with the circular holes arranged vertically are active.


### 3.2. Analysis of Crack Evolution

Based on the description of the total number of cracks in Figure4, to better show the crack evolution in each model, the crack growth in each model is analyzed from the perspective of the crack initiation stress, damage stress, peak stress, and 60% of the UCS stress level after model failure. The distinguishing crack evolution features of each model are shown in Figure 5. The tensile cracks are depicted in red, while the shear cracks are depicted in black.Figure 5
Crack evolution characteristics in each model: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)The crack initiation in the models with hole flaws primarily occurs at the top and bottom points of the circular holes, and each model has tensile cracks. The existence of the circular hole flaws is the main factor that induces the early initiation of the cracks. When the damage stress of each model is reached, the number of cracks in the intact model increases, the crack distribution becomes more dispersed, and a small number of shear cracks appear. The cracks in the NM-1 to NM-7 models are gathered around the circular holes, and the crack distribution patterns are distinct due to the influence of the varying number and distribution of circular holes. For example, in the NM-1 model, the cracks exhibit an “Ж”-shaped distribution centered on the single circular hole, while in the NM-2 model, the cracks are mainly distributed above and below the two circular holes and concentrated in the middle part of the vertical direction of the model, resulting in greater crack propagation in the middle part. A small number of shear cracks appear in each flawed model under damage stress, and they are concentrated near or between the circular holes. At the peak stress, cracks continue to develop and propagate along the distribution state of the damage stress in all the models. In the NM-1 through NM-7 models, the cracks at the upper and lower points of the circular holes spread slowly and are near the circular holes. These cracks propagate faster through the circular holes and penetrate the entire model. However, the models are in the limit equilibrium state, and the cracks are about to penetrate and destroy the model. The number of shear cracks shows a particular increase, and the distribution range is still mainly gathered around or between the holes at the peak stress. At the postpeak failure stage, the extensive cracks in each model develop rapidly and form a concentrated area of crack penetration, affecting each model’s failure path.


### 3.3. Variation of the Vertical Stress Field

In order to investigate the distinguishing stress field evolution features of each model, the vertical stress field is analyzed under various stress levels. 10%, 50%, 100%, and 60% of the UCS in the postpeak failure stage are selected to illustrate the distinguishing distribution features of the stress field in each model, as shown in Figure6.Figure 6
Distinguishing distribution features of the vertical stress field: (a) represents the intact model; (b–h) represent the NM-1 to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Figure6(a) indicates that the stress concentration of the intact model before the peak stress is uniformly distributed in the entire model and only slightly dispersed until peak stress. The stress concentration area in the postpeak stage shows a roughly “conical” distribution, and the main distribution area of the low stress is located in the lower right corner of the model. Before the peak stress, in the NM-1 to NM-7 models, regardless of the size and distribution of the circular hole, there is a compressive stress concentration phenomenon on both sides of the holes. In contrast, the area above and below the circular hole is a low-stress or tensile stress area. With the increase in the loading stress, the compressive stress concentration region on the two sides of the circular hole shows a decreasing trend, but the degree of stress concentration is greater. The low-stress or tensile stress area above and below the circular hole shows a gradually expanding trend. At the peak stress, the compressive stress concentration region of the NM-1 model is distributed in the form of a “butterfly” centered on the circular hole. The maximum concentration stress is still close to the two sides of the circular hole, and the maximum concentration stress is about 191.68 MPa. The maximum compressive stress concentration area between the two horizontal circular holes in the NM-2 model is gradually fused, and the maximum concentration stress is about 187.67 MPa. The maximum stress concentration area in the NM-3 model is located between the two 45° arranged circular holes and they are connected, while the stress concentration area outside the two circular holes diffuses outward to the model. In the NM-4 model, the large stress concentration region spreads away from the positions on both sides of the circular hole, focusing primarily on both sides of the model. The low-stress area above and below the hole also extends to an area between the top and bottom of the model and two holes. At peak stress, the stress concentration distribution of the NM-5 to NM-7 models is comparable to that of the NM-2 to NM-4 models. In the postpeak failure stage, the stress concentration areas in each model with hole flaws are significantly reduced, and the low-stress areas above and below the circular hole penetrate the entire model. It is worth noting that in the NM-4 and NM-7 models, the compressive stress has a wide range of stress concentration on both sides of the model, whether at the peak stress or after the postpeak stage. This demonstrates that the models have good bearing capacity, which is also the reason why both models have high uniaxial compressive strengths.


### 3.4. Displacement Field Evolution and Failure Mode

To analyze and compare the change of the displacement field and the distinguishing features of each model’s deformation and failure process, 50% of the UCS, peak stress, and 60% of the UCS in the postpeak failure stage of each model are selected. Variations in the displacement fields of each model are shown in Figure7.Figure 7
Variations in the displacement fields: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Figure7 indicates that under uniaxial compression, the displacement field of the intact model before failure is the smallest in the middle of the horizontal direction, and the displacement gradually increases towards the top and bottom of the model. The top and bottom are the maximum displacement area of the model. The displacement field is distributed symmetrically along the horizontal centerline of the model. During the failure period, the displacement at both ends of the bottom of the model is large, and the model’s small and medium displacement areas exhibit a “V” type distribution. The overall failure mode of the model is conical, which corresponds to the crack concentration area during the failure period of the intact model. When the NM-1 model is at 50% of the UCS, the displacement field is also symmetrically distributed along the model’s horizontal centerline, and areas with small displacements appear on either side of the circular hole. Additionally, the small-displacement change areas on both sides are distributed symmetrically along the circular hole in an “M” shape. When the load reaches its peak stress, the displacement of the small displacement area in the model gradually increases and spreads, forming an “X”-shaped small-displacement area. At the failure stage, the model displays inclined splitting failure along the vicinity of the circular hole. In the NM-2 to NM-4 models, the maximum displacement before peak stress is primarily distributed at the top and bottom of the model, and the minimum displacement is concentrated between the horizontal connections of the two circular holes, i.e., the distribution of the minimum displacement in the NM-3 model is at 45° between the two circular holes. At the peak stress, the small-displacement area is still concentrated between the two circular holes, whereas the displacement outside the two circular holes increases rapidly. For example, the displacement outside the two circular holes in the NM-3 model changes into a concentrated area of the maximum displacement. At this time, the displacement field distribution in the middle of the model is chaotic, indicating that large critical failure cracks occur in some areas, affecting the displacement field distribution. At the postpeak failure stage, NM-2 shows apparent fracture failure on both sides of the model (outside of the two holes) and in the middle of the circular hole. The NM-3 model is split near the two circular holes in the model’s lower left and upper right areas. In the small-displacement area of the NM-4 model at the failure stage, the main circular holes below show a small “V” distribution with vertices, while the entire model reveals fracture failure at the outer and lower circular holes on both sides of the “V” distribution to the bottom of the model.Before the peak stress, the displacement field variation pattern of the NM-5 to NM-7 models is similar to that of the corresponding two-hole models. The regions with small displacement are also mainly distributed between the holes. However, the displacement field distribution at the peak stress is more disorderly than that of the corresponding two-hole models due to the influence of the three circular holes. In the failure period, the NM-5 model undergoes fracture failure along the outer side of the middle hole and the circular holes on both sides, and the maximum displacement area is distributed outside the circular holes on both sides of the model. In the NM-6 model, there is a fracture connection between the three circular holes, and splitting failure occurs along the upward and downward circular holes. The fracture failure of the NM-7 model mainly occurs on both sides of the three circular holes, demonstrating a relatively broken splitting failure. The analysis of the displacement field evolution shows that the distribution and number of circular holes play a vital role in the evolution of the displacement field and the failure mode of the models.


## 3.1. Analysis of AE Events and Characteristic Stress Values of each Model

The stress-strain, AE, and the totality of cracks for each model are summarized together, and the three characteristic stress values of crack initiation stress, damage stress, and peak stress of each model are marked on the stress-strain curve. The stress-strain-acoustic emission-crack number curves of each model are obtained, as shown in Figure4. The crack initiation, damage, and peak stress are denoted by σci, σcd, and σc, respectively. Table 2 lists the values of the crack initiation, damage, and peak strain-stress of the models.Figure 4
Stress-strain-acoustic emission-crack number curves of each model: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Table 2
Crack initiation, damage, and peak stress and strain of each model.
Model nameCrack initiation strain (%)Crack initiation stress (MPa)Damage strain (%)Damage stress (MPa)Peak strain (%)Peak stress (MPa)Intact model0.062548.160.1402106.520.1727125.05NM-10.019213.500.093461.430.149884.68NM-20.022114.990.084155.590.118570.15NM-30.015610.540.081852.870.137670.51NM-40.022715.670.103468.440.158392.28NM-50.025717.230.080252.750.113169.52NM-60.017011.350.080549.960.115263.86NM-70.029220.270.109573.300.150194.35The curves in Figure4 reveal that the UCS of the intact model is about 125.05 MPa, while the uniaxial compressive strengths of NM-1 to NM-7 are about 84.68, 70.15, 70.51, 92.28, 69.52, 63.86, and 94.35 MPa, respectively. The strength of the intact model is the largest, and the UCS of NM-1 to NM-7 is significantly reduced. Given the numbers and distribution of hole flaws, the UCS of the models containing the circular hole flaws is different, and the UCS of the two models with the vertical arrangement of circular holes is high. The strength of NM-2 and NM-3 is equivalent, while the strength of NM-5 is greater than that of NM-6. This indicates that when the two circular holes are arranged at 45°, the influence of the related flaw types, such as fissures and soft structural planes, on the rock mass is different [57, 58]. The main reason for this is that although the two circular holes are arranged at a 45° dip angle, the spacing between the circular holes is still relatively large, and penetration under uniaxial compression is insufficient. When the three circular holes are arranged at 45°, the distance between them decreases, and the effect of promoting the failure of the rock samples increases, resulting in a reduction in its bearing capacity; therefore, the NM-6 model has the smallest UCS of all the models. In general, the strength of the two-hole model with the horizontal and 45° distribution of holes is stronger than that of the corresponding three-hole model. However, when the circular hole flaws are arranged in the vertical direction, the UCS of the model with the three vertical holes is the largest, and the intact model is reduced by 24.55%. The strength of the model with the two holes is the second smallest, and the UCS of the intact model decreases by 26.21%. The three-hole model with the 45° inclination angle has the smallest UCS, which is 48.93% lower than the intact model.In the AE characteristic histograms of each model, the AE events do not occur when the initial loading stress is minor. As the loading stress increases, cracks appear within the rock sample, and acoustic emission phenomena appear in each model. The crack initiation stress of each model is about 48.16, 13.50, 14.99, 10.54, 15.67, 17.23, 11.35, and 20.27 MPa, respectively. The hole flaws significantly reduce the crack initiation stress of the models and cause the AE events to occur earlier. Two-hole and three-hole models with the 45° arrangement begin to exhibit AE activity at strains of 0.0156% and 0.0170%. Internal cracks in the rock samples gradually develop and expand with continuous loading, and the number of cracks maintains a stable growth rate. The number of AE counts of each model also continues to increase slowly. When loading to 70%-80% of the peak stress, the totality of cracks in each model gradually begins to increase rapidly, and the number of AE also counts suddenly increases. At this time, the model is considered to have reached the damage state, and the corresponding stress level is the damage stress. The damage stress of each model is about 106.52, 61.43, 55.59, 52.87, 68.44, 52.75, 49.96, and 73.30 MPa, respectively. The damage stress of the models with circular hole flaws is significantly reduced compared to the intact model, while the damage strain is also significantly increased. In terms of the damage changes and the crack initiation stress, the variations under various distributions of the two and three circular holes are similar, with the smallest changes occurring when the holes are arranged at 45°. Accordingly, the damage degree is greater under the same loading condition, and the values are the largest when the holes are distributed vertically, resulting in a lower damage degree. After the damage stress of each model is reached, the total number of cracks in each model increases rapidly, and the number of acoustic counts rises sharply at the peak stress, where it reaches the maximum. The acoustic emission is “active” around the peak stress until the postpeak failure stage, and the AE events are concentrated around the peak stress. The number of cracks in each model increases rapidly again, but the total number of AE counts decreases in the postpeak stage. Overall, the intact model has the highest number of AE counts and the most concentrated AE event distribution. Each model’s maximum number of AE counts is small, and the AE events are scattered around the peak stress. The distribution and number of the circular hole flaws also influence the distinguishing AE features. They affect the early and late occurrence of AE events at the crack initiation stress level and the number of acoustic emission shots at the peak stress. The changing trend is equivalent to each model’s peak stress, and the model’s AE events with the circular holes arranged vertically are active.


## 3.2. Analysis of Crack Evolution

Based on the description of the total number of cracks in Figure4, to better show the crack evolution in each model, the crack growth in each model is analyzed from the perspective of the crack initiation stress, damage stress, peak stress, and 60% of the UCS stress level after model failure. The distinguishing crack evolution features of each model are shown in Figure 5. The tensile cracks are depicted in red, while the shear cracks are depicted in black.Figure 5
Crack evolution characteristics in each model: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)The crack initiation in the models with hole flaws primarily occurs at the top and bottom points of the circular holes, and each model has tensile cracks. The existence of the circular hole flaws is the main factor that induces the early initiation of the cracks. When the damage stress of each model is reached, the number of cracks in the intact model increases, the crack distribution becomes more dispersed, and a small number of shear cracks appear. The cracks in the NM-1 to NM-7 models are gathered around the circular holes, and the crack distribution patterns are distinct due to the influence of the varying number and distribution of circular holes. For example, in the NM-1 model, the cracks exhibit an “Ж”-shaped distribution centered on the single circular hole, while in the NM-2 model, the cracks are mainly distributed above and below the two circular holes and concentrated in the middle part of the vertical direction of the model, resulting in greater crack propagation in the middle part. A small number of shear cracks appear in each flawed model under damage stress, and they are concentrated near or between the circular holes. At the peak stress, cracks continue to develop and propagate along the distribution state of the damage stress in all the models. In the NM-1 through NM-7 models, the cracks at the upper and lower points of the circular holes spread slowly and are near the circular holes. These cracks propagate faster through the circular holes and penetrate the entire model. However, the models are in the limit equilibrium state, and the cracks are about to penetrate and destroy the model. The number of shear cracks shows a particular increase, and the distribution range is still mainly gathered around or between the holes at the peak stress. At the postpeak failure stage, the extensive cracks in each model develop rapidly and form a concentrated area of crack penetration, affecting each model’s failure path.


## 3.3. Variation of the Vertical Stress Field

In order to investigate the distinguishing stress field evolution features of each model, the vertical stress field is analyzed under various stress levels. 10%, 50%, 100%, and 60% of the UCS in the postpeak failure stage are selected to illustrate the distinguishing distribution features of the stress field in each model, as shown in Figure6.Figure 6
Distinguishing distribution features of the vertical stress field: (a) represents the intact model; (b–h) represent the NM-1 to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Figure6(a) indicates that the stress concentration of the intact model before the peak stress is uniformly distributed in the entire model and only slightly dispersed until peak stress. The stress concentration area in the postpeak stage shows a roughly “conical” distribution, and the main distribution area of the low stress is located in the lower right corner of the model. Before the peak stress, in the NM-1 to NM-7 models, regardless of the size and distribution of the circular hole, there is a compressive stress concentration phenomenon on both sides of the holes. In contrast, the area above and below the circular hole is a low-stress or tensile stress area. With the increase in the loading stress, the compressive stress concentration region on the two sides of the circular hole shows a decreasing trend, but the degree of stress concentration is greater. The low-stress or tensile stress area above and below the circular hole shows a gradually expanding trend. At the peak stress, the compressive stress concentration region of the NM-1 model is distributed in the form of a “butterfly” centered on the circular hole. The maximum concentration stress is still close to the two sides of the circular hole, and the maximum concentration stress is about 191.68 MPa. The maximum compressive stress concentration area between the two horizontal circular holes in the NM-2 model is gradually fused, and the maximum concentration stress is about 187.67 MPa. The maximum stress concentration area in the NM-3 model is located between the two 45° arranged circular holes and they are connected, while the stress concentration area outside the two circular holes diffuses outward to the model. In the NM-4 model, the large stress concentration region spreads away from the positions on both sides of the circular hole, focusing primarily on both sides of the model. The low-stress area above and below the hole also extends to an area between the top and bottom of the model and two holes. At peak stress, the stress concentration distribution of the NM-5 to NM-7 models is comparable to that of the NM-2 to NM-4 models. In the postpeak failure stage, the stress concentration areas in each model with hole flaws are significantly reduced, and the low-stress areas above and below the circular hole penetrate the entire model. It is worth noting that in the NM-4 and NM-7 models, the compressive stress has a wide range of stress concentration on both sides of the model, whether at the peak stress or after the postpeak stage. This demonstrates that the models have good bearing capacity, which is also the reason why both models have high uniaxial compressive strengths.


## 3.4. Displacement Field Evolution and Failure Mode

To analyze and compare the change of the displacement field and the distinguishing features of each model’s deformation and failure process, 50% of the UCS, peak stress, and 60% of the UCS in the postpeak failure stage of each model are selected. Variations in the displacement fields of each model are shown in Figure7.Figure 7
Variations in the displacement fields: (a) represents the intact model; (b–h) represent the NM-1 model to the NM-7 model, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Figure7 indicates that under uniaxial compression, the displacement field of the intact model before failure is the smallest in the middle of the horizontal direction, and the displacement gradually increases towards the top and bottom of the model. The top and bottom are the maximum displacement area of the model. The displacement field is distributed symmetrically along the horizontal centerline of the model. During the failure period, the displacement at both ends of the bottom of the model is large, and the model’s small and medium displacement areas exhibit a “V” type distribution. The overall failure mode of the model is conical, which corresponds to the crack concentration area during the failure period of the intact model. When the NM-1 model is at 50% of the UCS, the displacement field is also symmetrically distributed along the model’s horizontal centerline, and areas with small displacements appear on either side of the circular hole. Additionally, the small-displacement change areas on both sides are distributed symmetrically along the circular hole in an “M” shape. When the load reaches its peak stress, the displacement of the small displacement area in the model gradually increases and spreads, forming an “X”-shaped small-displacement area. At the failure stage, the model displays inclined splitting failure along the vicinity of the circular hole. In the NM-2 to NM-4 models, the maximum displacement before peak stress is primarily distributed at the top and bottom of the model, and the minimum displacement is concentrated between the horizontal connections of the two circular holes, i.e., the distribution of the minimum displacement in the NM-3 model is at 45° between the two circular holes. At the peak stress, the small-displacement area is still concentrated between the two circular holes, whereas the displacement outside the two circular holes increases rapidly. For example, the displacement outside the two circular holes in the NM-3 model changes into a concentrated area of the maximum displacement. At this time, the displacement field distribution in the middle of the model is chaotic, indicating that large critical failure cracks occur in some areas, affecting the displacement field distribution. At the postpeak failure stage, NM-2 shows apparent fracture failure on both sides of the model (outside of the two holes) and in the middle of the circular hole. The NM-3 model is split near the two circular holes in the model’s lower left and upper right areas. In the small-displacement area of the NM-4 model at the failure stage, the main circular holes below show a small “V” distribution with vertices, while the entire model reveals fracture failure at the outer and lower circular holes on both sides of the “V” distribution to the bottom of the model.Before the peak stress, the displacement field variation pattern of the NM-5 to NM-7 models is similar to that of the corresponding two-hole models. The regions with small displacement are also mainly distributed between the holes. However, the displacement field distribution at the peak stress is more disorderly than that of the corresponding two-hole models due to the influence of the three circular holes. In the failure period, the NM-5 model undergoes fracture failure along the outer side of the middle hole and the circular holes on both sides, and the maximum displacement area is distributed outside the circular holes on both sides of the model. In the NM-6 model, there is a fracture connection between the three circular holes, and splitting failure occurs along the upward and downward circular holes. The fracture failure of the NM-7 model mainly occurs on both sides of the three circular holes, demonstrating a relatively broken splitting failure. The analysis of the displacement field evolution shows that the distribution and number of circular holes play a vital role in the evolution of the displacement field and the failure mode of the models.


## 4. Discussion

This study used the PFC2D code to establish limestone models with different numbers and distributions of circular hole flaws under seven operating conditions. It considered the arrangement and number of holes as the main variables. If the number of holes increases, the distribution of holes will be involved. To better represent the effect of the number of holes on the damage and distinguishing failure features of the rock samples, the total flaw area in the model was kept consistent, and the distance between holes and the model boundary was the same; most researchers did not consider these conditions. However, no specific analysis of other multihole models with different inclination arrangements, such as 30° or 60°, was performed to facilitate this study.The models with various numbers and distributions of circular hole flaws in this research and the models in related studies [23–25, 27, 30, 31, 34] have some similar distribution arrangements. Although the aperture sizes of the circular holes set in previous models are frequently identical, the derived mechanical properties and deformation and failure behaviors of rocks with the same number and arrangement of hole defects are comparable to the results obtained, indicating the reliability of this study. Zhang et al. [59] also utilized the PFC2D code to investigate the changes in the mechanical parameters of the sandstone model with two circular holes arranged at various inclination angles and analyzed the damage variables of the models using the AE characteristics. However, the established models only contain two circular holes, large and small, without analyzing the evolution law of the entire process of the cracks and using the AE characteristics to determine the damage stress of the models. Based on the difference in the number and distribution of holes, previous studies rarely addressed the effect of crack initiation stress and damage stress of defective rock samples. The study’s results demonstrated that under the same number of circular holes, the crack initiation stress and damage stress of the model change with the distribution of holes in the shape of a “V”, providing a valuable reference for analyzing the damage and failure of rock masses.In addition, since the hole shape in this study was only circular, it cannot represent the complexity of natural rock masses with irregular defect shapes and uneven distribution of other defect shapes or types. However, based on the premise of the same flaw area, the variation of the mesomechanical behavior of a rock mass, such as the damage evolution caused by different numbers and distributions of circular hole flaws, was comprehensively analyzed. To a certain extent, this study provides insights into the interaction between the number and distribution of holes and their influence on the complex mechanical responses of rock masses.


## 5. Conclusions

(1)
The number and distribution of circular hole defects greatly influence the peak stress, crack initiation stress, damage stress, and acoustic emission characteristics of rock samples. Due to the varying number and distribution of circular hole flaws, the peak stress, crack initiation, and damage stress of the single-hole model are between those of multihole models. For the models with two or three circular holes, the peak stress, damage stress, and crack initiation stress exhibit a “V”-shaped change as the angle of the hole arrangement rises, and the characteristic stress values are the highest when the hole flaws are arranged vertically(2)
The presence of hole flaws causes crack initiation at the top and bottom of the circular holes in the models, and the two models with circular holes arranged at 45° begin to crack earlier. In the process of each model reaching its peak stress, the crack development of the model with circular hole defects concentrates primarily near or between circular holes and extends outward from near or between circular holes until the entire model is essentially coalesced. Variations in the path, shape, and number of crack development are primarily due to differences in the number and distribution of circular holes(3)
Before the peak stress, stress concentration areas are formed near the circular holes in the models. With the different distribution of the circular hole defects, the stress release areas or the tensile stress areas remain concentrated at the upper and lower vertices of the circular hole, whereas the compressive stress concentration areas of the multihole arranged at 45° rotate with the connection area of the circular hole’s center. The stress concentration areas in the two models with vertical circular holes are still distributed on both sides of the models after their failure(4)
The number and distribution of the circular holes play a vital role in the evolution of the displacement field and the failure mode of the models. Before the peak stress, the small displacements of circular-holed models are concentrated near and between the circular holes and vary with the arrangement angle of the circular holes. In the failure stage, the fracture failure path of the entire model is affected by the different distributions of the circular holes, resulting in various failure patterns of the rock samples with distinct numbers and distributions of hole flaws



---

*Source: 9382577-2022-12-05.xml*