Abstract
The existence of flaws in brittle rocks or rock-like materials has an obvious influence on the material mechanical properties and cracking behavior of civil engineering projects. In this work, the two-dimensional particle flow code PFC2D was used to study the deformation and strength properties, failure processes, and acoustic emission (AE) characteristics of mudstone with a single preexisting flaw. First, the procedure to construct a parallel bond model is introduced. The Weibull distribution is used to reflect the mechanical heterogeneity of rocks. Then, the microscopic parameters used in PFC2D are calibrated to the macroproperties of mudstone obtained from laboratory tests under the uniaxial compression. The results indicate that the increases of the flaw inclination lead to the increasing uniaxial compressive strength and elastic modulus. In terms of microcrack evolution, the initiation, propagation, and coalescence of microcracks are closely related to the force chain. Specifically, an “X” shaped tension force chain concentrated area around the preexisting flaw is founded, which is the most prone area for microcracks to initiate. With an increase in flaw inclination, the b value of AE also shows an increasing trend. By incorporating the AE event numbers into a damage variable, this paper derives a constitutive model, which is verified by numerical results on brittle rocks with a single preexisting flaw under uniaxial compression.
1. Introduction
With the development of productivity and technology, increasingly open-pit coal mines have formed a series of high and steep slopes. The slope contains a rich network of fractures, and the stability of these slopes greatly affects the production and life of open-pit coal mines and the safety of personnel. At present, for low-strength fractured rock masses, such as mudstone, there are relatively few studies on the mechanical properties and microcrack propagation [1]. In an engineering project, such as open-pit coal mines and tunnels, rocks often contain more primary joints and flaws, which are more likely to cause damage to the rock masses. Therefore, it is necessary to conduct a more detailed study on the mechanical properties and microcrack propagation of mudstone containing a single preexisting flaw.
Because it is difficult to obtain and measure the internal fracture network of the rock masses, it has become a feasible research method to study the preexisting flaws with a small number. In terms of experiments, Yang and Jing conducted uniaxial compression tests of single-flaw sandstone and found that its mechanical properties have a certain correlation with the inclination of the fracture [2]. Wong and Einstein used gypsum prefabricated single-flaw specimens for uniaxial compression tests and obtained the propagation and evolution characteristics of wing cracks and secondary cracks [3]. Aliabadian et al. used 3D printing technology and digital image correlation technology, and the mechanical properties and strain field laws of rock-like materials with single and double flaws in Brazilian testing were studied [4]. Zhao et al. found that the evolution of the strain field was consistent with the initiation, propagation, and coalescence of microcracks [5]. Li et al. conducted uniaxial compression tests on a coal sample with a single preexisting flaw and conducted acoustic emission (AE) signal monitoring during the loading processes [6].
With the rapid development of computer technology, numerical simulation technology has become an essential technical method to study the mechanical properties of fractured rocks. As an effective supplement to laboratory tests, it plays an important role in the field of rock mechanics. Zhao et al. used an expanded distinct element method (EDEM) to simulate the crack propagation of a single fractured rock and found that with the change of the inclination of the fracture, the failure mode of the rock also changed accordingly [5]. Guo et al. used FLAC software to simulate the changes in the displacement field in tunnels and foundation pits [7, 8]. Bi et al. used the general particle dynamics (GPD) method to study the crack growth of rock-like materials containing multiple preexisting flaws under different confining pressures under three-dimensional conditions [9]. The discrete element method (DEM) proposed by Cundall and Strack [10] is a very widely used method in numerical simulation. particle flow code (PFC) [11] is commercial software developed based on DEM. To date, many researchers have used PFC to study the mechanical properties and the growth of microcracks in rock materials or rock-like materials with preexisting flaws. Zhang et al. used PFC to study the difference in AE characteristics of single-flaw rock-like materials under different loading rates [12]. Chen et al. used PFC to study the microcrack growth and mechanical properties of the rock containing nonpersistent joints [13].
The main objective of this work is to propose a parallel bond model with the Weibull distribution to model crack propagation in heterogeneous brittle rocks, such as mudstone. In the proposed framework, the micromechanical properties of brittle rocks are specifically defined through the Weibull distribution. This particular setting ensures that material heterogeneities at the grain scale and allows for the crack propagation more random at the microscopic level. In this paper, PFC2D 5.0 is adopted to investigate the mechanical properties, cracking behavior, and AE characteristics of mudstone under uniaxial compression. The main structure of this paper is as follows: in Section 2, the parallel bond model, Weibull distribution, microparameter calibration, and acoustic emission based on moment tensor are introduced. In Section 3, mechanical properties, cracking behavior, and damage constitutive model of mudstone under uniaxial compression are studied. And conclusions are given in Section 4.
2. Numerical Model
2.1. Establishment of the Numerical Model
There are two basic bonding models in PFC, the contact-bond model and the parallel bond model (PBM). In the simulation of rocks or rock-like materials, PBM is widely used. Previous-related numerical studies [12–14] showed that PBM is a suitable contact model in PFC for the simulation of the mechanical properties and failure modes of rocks. In our study, the particle deletion method was adopted to generate an open single preexisting flaw by deleting the particles in a specified region, as shown in Figure 1. The flaw inclination α increases from 0°, 15°, 30°, 45°, 60°, 75° to 90°.
(a)
(b)
2.2. Rock Heterogeneity
As an aggregate composed of multiple groups of mineral grains, the mechanical properties of the rock are affected by the type, size, shape, and primary defects of minerals, making the rocks become a typical heterogeneous material [15]. Therefore, in the processes of numerical simulation, it is necessary to consider the heterogeneity of the rocks.
Many researchers have adopted the Weibull distribution to characterize the heterogeneity of mechanical parameters of rock in numerical simulations [14, 16–18] and theoretical analysis [19, 20]. Therefore, in our study, two parameters, i.e., parallel bond tensile strength and parallel bond cohesion, were assigned following the Weibull distribution which was defined in terms of the characteristic strength [14, 16–18]. For example, the density distribution function for the parallel bond tensile strength (σc) iswhere is the scale parameter giving the characteristic value of the parallel bond tensile strength distribution σc and m is the shape parameter. With increasing m, the generated data are more concentrated.
2.3. Microparameter Calibration
The microparameters used in the PBM (parallel bond modulus, parallel bond tensile strength, parallel bond cohesion, etc.) are generally difficult to obtain directly through laboratory tests. Instead, the macromechanical properties such as uniaxial compressive strength (UCS), elastic modulus (E), and Poisson’s ratio () obtained from laboratory uniaxial compression tests are constantly compared, and the “trial-and-error method“ is used to calibrate the microparameters. Finally, the established numerical model of PFC2D of the intact sample is consistent with the macromechanical properties of the intact sample obtained in laboratory tests. This parameter calibration method has been widely used and has achieved very good results [14, 21–23]. Many researchers have studied the relationship between the microparameters of the PBM and the macromechanical properties of the sample in laboratory tests, which provides a lot of important information about parameter calibration [18, 24].
The microparameters of the intact model established in this paper are shown in Table 1. Figure 2 shows the density distribution of the parallel bond tensile strength with the values of and m being 11 MPa and 5 in Table 1, respectively. The comparison between the numerical model and the laboratory tests [25] are shown in Table 2 and Figure 3.
Then, the calibrated intact model is used to create numerical models with a single preexisting flaw by deleting the particles in a certain region (as shown in Figure 1). The flaw inclination α increases from 0°, 15°, 30°, 45°, 60°, 75°, to 90°.
2.4. Acoustic Emission Simulation by PFC2D
The moment tensor, a well-known method, was used for the quantification of the mechanics of a seismic source [26]. The principle is briefly described as follows (for a more detailed description, refer to the literature [26]). In our study, the AE magnitudes are calculated by the moment tensor:where mj is the jth eigenvalue of the moment tensor matrix.
The magnitude M of the AE event based on the maximum value of the scalar moment is calculated as follows:
In PFC2D, the AE simulation mainly sets two parameters: the acoustic quality factor Q and the shear wave velocity Vs, the values of which are 200 and 3500 m/s, respectively [26].
3. Numerical Results and Discussion
3.1. Mechanical and Microcrack Properties of Model Containing a Single Preexisting Flaw
3.1.1. Stress-Strain Curve Characteristics
The mechanical properties of models containing a single preexisting flaw are shown in Figure 4. It can be seen from Figure 4 that the mechanical response of different models under axial loading is different. As the angle increases, the UCS of each model increases, as shown in Figure 4(b). And the axial strain corresponding to UCS also increases, which is consistent with the numerical research results of Jin et al. [27] and the experiment results of Li et al. [6]. It shows that the increase of the preexisting flaw angle makes the model’s ability to resist external deformation also increased. The UCS and E of the intact model are 21.2 MPa and 4.16 GPa, respectively. Compared with the intact model, the UCS and E of the model with a single preexisting flaw are significantly reduced. For example, when α = 0°, UCS, and E are 13.8 MPa and 3.82 GPa, respectively, with a decrease of 34.9% and 8.17%; when α = 60°, UCS and E are 18.55 MPa and 4.05 GPa, respectively, with a decrease of 12.5% and 2.64%. As shown in Figure 4(c), with the increase of the flaw inclination, the elastic modulus and Poisson’s ratio of the numerical model also tend to increase roughly. In the seven models established, the minimum and maximum values of the elastic modulus are 3.82 GPa and 4.14 GPa, respectively, with an increase of 8.37%.
(a)
(b)
(c)
3.1.2. Evolution of Microcracks
The evolution of axial stress and total microcrack number is recorded during loading, as presented in Figure 5. At the beginning of loading, no microcracks were initiated inside the model; as the loading progressed, a small amount of microcracks were initiated, when the axial strain was 2 × 10−3 approximately. Subsequently, the number of microcracks increased slowly. When the axial stress reached about ∼95% of the peak strength, the number of microcracks began to increase sharply, which proved that the microcracks rapidly expanded and overlapped each other, forming a macroscopic rupture zone and finally leading to the failure of the model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
In the PFC2D model, the force is reflected by the contact force chain between the balls. The width of the force chain line is related to the magnitude of the force; the wider the force chain, the greater the force. The microcrack propagations and the “force chain” in the PFC2D of the seven models are shown in Table 3, corresponding to points A ∼ E in Figure 5.
It can be seen from Table 3 that the initiation, propagation, and coalescence of microcracks are consistent with the evolution of force chains. In PFC, when the contact stress between particles (the value of which is the contact force divided by the contact area) exceeds the parallel bond contact tensile strength (or shear strength) of the contact, a tensile microcrack (or shear microcrack) will be initiated [11]. Because the model established in this work contains the heterogeneity of mechanical microparameters, its parallel bond contact tensile strength and parallel bond cohesion obey the Weibull distribution (as shown in Figure 2). Therefore, in the initial stage of loading, some contacts are due to low strength, and there will be some microcracks randomly distributed inside the model (as shown in Table 3, α = 0° model, point A). At this time (points A ∼ B), there is a “nucleus” with a smaller compression force chain around the preexisting flaw, and the area of this “nucleus” gradually decreases with the increase of the preexisting flaw inclination. When α = 90°, it almost disappears. Therefore, the ability of its model to resist external deformation has gradually increased, and the UCS of the model has also increased, which is consistent with the conclusion in Figure 4. At the same time, it can be seen in the tension force chain diagram that there is an “X” shaped force chain concentration area, and the center point is near the center of the preexisting flaw. Therefore, the two ends of the preexisting flaw are the most prone to the initiation of tension cracks, and the tension cracks are prone to continue to expand toward the tension force chain concentration area, as shown in point B in Table 3. The model of α = 0°∼45° has a certain length of macroscopic cracks at point B. The model of α = 60°∼90° does not have this phenomenon. When the axial stress continues to increase to the peak strength of the model (at point C in Figure 5), the model of α = 0° also has a macroscopic crack at the right end of the preexisting flaw. The crack extends to the upper boundary of the model and finally makes the model occur split failure.
It can be seen from Table 3 that, as the flaw inclination of the preexisting flaw increases, its failure mode also changes. When the flaw inclination is small, the failure mode presents a typical split failure, such as the model α = 0°∼15°; when the flaw inclination is large, the model has a typical shear failure, such as the model α = 60°∼75°; the failure surface is along the preexisting flaw surface for the model with α = 90°, and the cracks are also randomly distributed.
3.1.3. Fractal Dimension of the Microcrack Distribution
The fracture surfaces in rocks can provide valuable information about fracture mechanisms [28]. The fractal dimension offers a value for quantitatively describing the morphology of the failure fracture surfaces. The ImageJ image processing program (https://imagej.net/) is used for digital image processing [29]. The fractal dimension Df is calculated as follows:where δ is the scale length of a mesh box and N(δ) is the number of boxes. The fractal dimension Df is defined as the gradient of the log-log plots of N(δ) and δ.
The microcrack image at point E in Table 3 is imported into the ImageJ image processing program to calculate the fractal dimension, and the results are shown in Figure 6 and Table 4. It can be seen from Figure 6 that the fractal dimension of the failure surface of the models with different preexisting flaw inclination angles is different, which proves that the propensity of microcrack propagation is different, and the stress distribution and magnitude are also different [28]. As the flaw inclination increases, its fractal dimension roughly shows an increasing trend. The fracture surface fractal dimension of the α = 90°model is the highest, reaching 1.475, which proves that the α = 90° model has the roughest fracture surface.
(a)
(b)
3.2. AE Characteristics
3.2.1. b Value
The AE characteristics of seven numerical models calculated using the moment tensor are shown in Figures 7 and 8. It can be seen from Figures 7 and 8 that the number of AE events and the axial stress have good consistency. As the axial stress increases, the AE number gradually increases, and after reaching the peak strength, it increases rapidly. And the distribution of AE numbers with magnitude accords with the normal distribution. The b value of AE can be used to indicate the scale of fracture development when the rock ruptures. The larger the value of b value, the greater the proportion of small-scale fractures in the rock; the smaller the value of b value, the greater the proportion of large-scale fractures in the rock. It can be seen from Figure 8 that the b value of the model with α = 75° is the largest, reaching 2.29, which indicates that the proportion of small-scale fractures in the model with α = 75° is the highest. That is, the final macroscopic fracture surface of the model is composed of many small microfractures. The b value of the model with α = 0° is 1.99, which indicates that its large-scale fracture ratio is relatively high.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
3.2.2. Constitutive Model Based on AE Events
To study the relationship between deformation and strength of rock materials, constructing a damage constitutive model is one of the more basic and effective methods. Therefore, based on the results of laboratory tests or numerical simulations, the introduction of damage variables is to reflect the difference in the degree of internal damage of the rock during the loading processes. There are many methods to define damage variables, such as based on elastic modulus [30], number of cracks [31], and AE characteristics [32]. In our research, the damage variable D is defined by the number of AE events simulated by the theory of moment tensor:where Cn is the total number of AE events of the rock from the original loading point to a specific time and C is the total number of AE events of the rock.
Due to the different test control methods or different rock failure conditions, it is possible that the damage variable D cannot reach 1. Therefore, it is necessary to introduce the damage critical value variable Dn, so the damage variable is corrected towhere σc and σp, are the residual strength and peak strength, respectively.
Therefore, the constitutive equation of shale damage established based on the number of AE events iswhere σ is the axial stress; E is the elastic modulus; and ε is the axial strain.
The results calculated by equations (5)∼(7) are shown in Figures 9 and10. As shown in Figure 9, overall, the damage evolution law of each model is similar. In short, it can be divided into two stages: the slow increase in the early stage and the rapid increase in the later stage. As the loading progresses, the damage variable D increases continuously from 0 to 1. And the growth rate is slow before about 0.13, and the damage variable increases rapidly after 0.13. At the same time, the damage variable of the model with α = 0° has the fastest growth rate, indicating that the damage accumulation inside the model is fast, so the UCS of the model will also be lower.
(a)
(b)
After exceeding 0.13, all models are destroyed quickly. Due to the axial stress, a large number of microcracks in the model are initiated (such as points D ∼ E in Table 3). The microcracks rapidly propagate toward the stress concentration area or weak points of mechanical properties and overlap each other to form macrocracks, resulting in the number of AE events increasing rapidly, so that the value of the damage variable D also increases rapidly. When the damage variable D enters the rapid increase stage, it indicates that the model is about to lose its carrying capacity and will be in failure.
According to equation (5), the theoretical curves of different models can be obtained. The comparison between the theoretical curves and the numerical simulation is shown in Figure 10. It can be seen from Figure 10 that the damage constitutive model of mudstone with a single preexisting flaw established in this work according to the AE event numbers is consistent with the numerical simulation values. In the initial loading stage, the elastic stage of the stress-strain curve fits well. The postpeak stress drop is in good agreement, but the fitting effect to the peak strength is average. Therefore, within the allowable error range, it can be considered that the damage constitutive model constructed by the AE event numbers obtained from the AE phenomenon of rock based on the moment tensor simulation can better describe the complex relationship between the deformation and strength of mudstone containing a single preexisting flaw with different flaw inclination under uniaxial compression.
In the mining of various open-pit coal mines, low-strength rocks such as mudstone are more common and usually contain a large number of joints and cracks, which pose a certain threat to production and people’s life. Therefore, the research motivation of this work is to obtain more credible data through the PFC simulation under the condition of limited experimental data obtained from a small amount of rock samples on site and to build large-scale coal mines in the later stage. The stability analysis of the slope has laid a good foundation. At the same time, it is shown that certain AE monitoring facilities can be deployed at the open-pit coal mine site to monitor the microseismic signal inside the rock slope, and the damage constitutive model can be constructed through the AE characteristic parameters to guide the on-site coal mining activities. The efficient and safe mining of coal mines provides certain technical support and guarantees.
4. Conclusions
In this paper, the PBM with the Weibull distribution calibrated by laboratory tests is used to construct models containing a single preexisting flaw with different flaw inclinations to analyze the deformation and strength properties, the propagation of microcracks, the fractal dimension of the failure surface, and the acoustic emission characteristics. The main conclusions of this work can be summarized as follows:(1)The PBM coupled with the Weibull distribution can well simulate the heterogeneity of mechanical properties inside the rock. The UCS, elastic modulus, and Poisson’s ratio of a rock increase with the increase of flaw inclination.(2)Good agreement is maintained between microcrack development and axial stress. The initiation, propagation, and coalescence of microcracks dominate the distribution of force chains in rocks. In the compression force chain, there is a weak nucleus around the preexisting flaw. Correspondingly, in the tension force chain, there exists an “X” shaped force chain concentration area, which is dominant in the initiation and propagation of microcracks. And with the change of preexisting flaw inclination, the force chain concentration area will also rotate. The fractal dimension of the fracture surface also increases with the increase of flaw inclination.(3)The AE phenomenon of rock is closely related to the development of microcracks. The AE event numbers are small at the early loading stage, but they start to increase sharply near the peak strength, indicating the imminent failure of rocks. The magnitude and frequency of AE events were normally distributed. With the increase of flaw inclination, the AE b value increases substantially. The model of α = 75° has the largest b value (b = 2.29).(4)The damage constitutive model of rock based on AE event numbers can fit the test results well, and it also shows that the AE phenomenon of the rock simulated by the principle of moment tensor is correct. In the mining process of the open-pit coal mine, AE monitoring equipment can be arranged to monitor the microseismic signals inside the rock, to build the damage constitutive model, study the damage evolution law of the high slope formed by open-pit mining, and ensure the safe production.(5)Based on the impacts of preexisting flaws on mechanical properties, cracking behavior, and AE characteristics of rocks, engineering geologists could evaluate the rock slope stability by observing the property of preexisting flaws. For instance, the rock slope with cracks parallel to the direction of the maximum principal stress is relatively stable.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (No. 51778575) and the Special Project on Key Areas for General Universities in Guangdong Province (New Generation of Information Technology) (No. 2021ZDZX1116).