Abstract
Rock has the characteristics of natural heterogeneity and discontinuity. Its failure phenomenon induced by external force involves complex processes, including the microcrack initiation, propagation, coalescence, and the macrocrack formation. In this study, the Weibull random distribution based on the rock microstructure characteristics is introduced into the combined finite-discrete element method (FDEM) to establish the heterogeneous rock model, and the mechanical response and damage evolution of rock samples in uniaxial compression test are simulated. The results show that FDEM simulation with loaded heterogeneous rock model can reflect the progressive development of rock damage, fracture, and acoustic emission (AE) activity in real rock well. Meanwhile, the statistical analysis indicates that the number and energy evolution of AE events with different fracture modes in the model are consistent with the macroscopic failure mode of rock. The change of b-value also agrees with the increasing trend of high-energy events in the loading process. This method provides a new tool for the analysis of rock damage and fracture evolution.
1. Introduction
Rock is a kind of natural heterogeneous and discontinuous medium. Under the influence of external force, microcracks in rock gradually initiate, propagate, nucleate, coalesce, and finally form the macrocracks. The complex process corresponds to a trans-scale mechanical behavior from the microcosmic to the macroscopic, which involves the rock deformation, fracture, interaction of rock blocks, and further fragmentation in interaction process. It is difficult to quantitatively describe the process with theoretical analysis. The combined finite-discrete element method (FDEM) proposed by Munjiza et al. [1, 2] combines the advantages of finite element and discrete element, in which the continuous mechanical behaviors, such as elastic-plastic deformation of materials, can be simulated by the finite element method, while the discontinuous deformation behaviors such as damage and fracture can be analyzed by the softening failure of crack element and the contact calculation of blocks used in the discrete element method. The method can simulate the failure process of rock from continuity to discontinuity.
In this study, a new modelling method of rock heterogeneity in FDEM is proposed referring to Lisjak et al. [3], and rock damage evolution and acoustic emission (AE) phenomenon in uniaxial compression test are simulated.
2. Fundamental Principles of FDEM
In FDEM, the simulated material is discretized as a mesh consisting of nodes and triangular finite elements, and four-node cohesive quadrilateral crack elements are inserted between adjacent triangular elements. Before fracture formation, the stress and strain of the triangular element are calculated with finite element method, the node displacement is calculated with the forward difference method based on the time-explicit integration, and the fracture initiation and growth are reflected with the softening and failure of the crack element. After fracture formation, the interaction of fracture surfaces and the displacement of discrete blocks can be calculated by the discrete element contact method. The core algorithms of FDEM method contain the governing equations, deformation and failure of crack element, contact detection, and interaction of discrete element. The following is a brief introduction to the fundamental principles of FDEM [4, 5].
2.1. Governing Equations
In FDEM, the solid is discretized as triangular finite elements. At the beginning of each simulation time step, the velocity and displacement of nodes are obtained by solving the motion equation with the explicit second-order forward difference time integration method, and the nodal coordinates are updated. The governing equations can be expressed as follows:wherein M and C are mass matrix and damping diagonal matrix of the system, respectively; x is the vector of nodal displacements; F is the vector of nodal unbalanced forces which consist of the internal resisting forces, applied external loads, and contact forces.
2.2. Deformation and Failure Modelling
FDEM can simulate the damage and fracture process of materials by embedding crack elements without initial thickness between the edges of all adjacent triangular elements pairs and realizing the nonlinear mechanical behavior and state transformation from continuity to discontinuity through the softening and failure of crack elements. On the other hand, as shown in Figure 1 [3], the virtual crack model is adopted in FDEM to reflect the mechanical property changes in region near the crack tip. It is considered that there exists a fracture process zone (FPZ) which forms ahead of the crack tip due to interlocking and microcracking. Before the appearance of macrocracks, the typical stress-strain curve can be divided into two parts: the hardening branch before the peak stress and the strain softening partthat represents decreasing stress with increasing strain. The crack element breaks when its strain reaches or exceeds the limit value. According to differences of ultimate strain states when crack elements break, three fracture modes are classified, i.e., the tensile (mode I), the shear (mode II), and the mixed (mode I-II). Constitutive behavior and failure criterion of crack elements in different fracture modes are shown in Figure 2.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
2.3. Contact Interaction of Discrete Element
The NBS (No Binary Search) algorithm proposed by Munjiza [1] was used in contact detection of this work. Meanwhile, parallel algorithm can be also applied to improve detection efficiency greatly.
According to the contact theory, it is generally considered that the contact force between two discrete elements is proportional to their invasion area between each other and acts as the distributed force. To obtain computable contact interaction, the distributed contact force between interactive solids is simplified into the nodal forces of triangular elements in FDEM (as shown in Figure 3). Detailed calculation methods can be found in the monograph written by Munjiza [1].
(a)
(b)
3. AE Simulation and Introduction of Material Heterogeneity
3.1. Definition of AE Events in Numerical Method
Acoustic emission is a physical phenomenon of internal energy release in the form of high-frequency transient elastic waves when rock is subjected to the sufficiently high stress field [6]. For brittle rocks, AE events are generally induced by microcrack initiation, pore fracture, or grain boundary sliding. With the development of numerical simulation technologies, AE activity of the loaded rock can be simulated well and AE parameters can be also analyzed quantitatively.
3.1.1. AE Simulation in Continuous Model and Discontinuous Models
At present, some scholars have conducted AE simulation of brittle materials in failure process with continuous or discontinuous medium model. In the continuous medium model, solid failure is usually simulated by weakening mechanical properties of broken elements, like strength and stiffness parameters. The element fragmentation plays an important role in AE sources and the released elastic energy corresponds to the energy of AE event. This method was proposed by Tang and implemented by his finite element code, RFPA (rock failure process analysis) [7]. In the discrete element model of discontinuous medium, PFC (particle flow code) method simulates the elastic dynamic effects such as stress wave propagation and crack induced AE activity by monitoring the movement of particles. Hazzard and Young [8] defined the radiation energy of AE events with the change of particle kinetic energy in relevant region when the bond breaks, and the b-value change of AE evolution was also analyzed to reflect the feature of rock damage and fracture development. Zhang et al. [9, 10] also studied the damage evolution and AE activity of rock under various loading conditions using the abovementioned cohesive particle model.
3.1.2. AE Event Definition of Crack Element Fragmentation in FDEM Model
In FDEM, the strain energy of materials is stored in triangular elements due to their elastic deformation. When stress states of some local crack elements reach their yield conditions, then new fractures are initiated and the stored strain energy in rock begins to release gradually. Among the released energy, a part is transformed into the energy input needed in crack element fragmentation process from the yield stage, a part is dissipated by the possible friction between adjacent triangular elements after fracture occurrence, and the rest is transmitted into the radiation of kinetic energy which can be regarded as AE event energy.
An AE event can correspond to a crack element breakage. Therefore, the mass center coordinates of crack element when entering the yield state can be considered as the AE source location, and its fracture mode is determined according to the relative displacement of fracture element boundaries. The energy of AE event is equivalent to the maximum kinetic energy change of the broken crack element from entering yield state to complete failure.
Taking the mode I fracture of crack element as an example, its evolution of normal bonding stress and nodal kinetic energy are shown in Figure 4. In the process of crack element failure from the moment when it enters the yield stage to the moment when it reaches the ultimate deformation, its normal bonding stress decreases gradually, while the kinetic energy of element increases at its early stage and then decreases after the peak at the moment . The phenomenon that the kinetic energy of crack element decreases after it reaches the maximum value can be attributed to the continuous radiation of kinetic energy and the influence of damping. The calculation of AE event energy in FDEM also refers to the algorithm developed by Hazzard and Young for discrete element based on particle model [8]. The overall calculation process is as follows:(1)When the crack element stress reaches the peak strength (the moment ), the kinetic energy of its four nodes is calculated and stored as : wherein is the nodal mass and is the nodal velocity at the yielding moment.(2)The kinetic energy of crack element is monitored continuously until it fails. The change of element kinetic energy in each time step is calculated as follows: wherein is the kinetic energy of the crack element at the moment t.(3)The maximum change of in the whole process of crack element failure from the yielding moment to the failing moment can be regarded as the seismic energy of AE event induced by each crack breakage. It is calculated as follows and the initiation time of AE event corresponds to the moment :(4)Finally, the event magnitude, , can be calculated according to Gutenberg [11]:
3.2. Introduction of Material Heterogeneity
The heterogeneity of natural rock includes the nonuniform existence and distribution of microstructure, such as natural defects, holes, and joints, as well as the nonuniform mechanical properties of these microstructures. It causes the discrete fracture response of rock in stress field and the progressive damage development before the critical instability [7]. However, the equivalent uniform distribution of material properties usually leads to no local element damage and AE activity before some macrocracks occur [3]. It is inconsistent with the AE monitoring results of rock failure tests. Therefore, it is necessary to define the heterogeneity of materials for simulation of gradual AE development in rock loading process.
3.2.1. Heterogeneous Model of Rock Material Based on Random Distribution
In RFPA, to simulate the nonlinear mechanical response in rock, Tang [7] established heterogeneous continuum model, in which the mechanical parameters (such as failure strength and elastic modulus) of the material are assumed to obey Weibull distribution. Its application contributes to successful simulation of the gradual rock failure process, including microcrack coalescence, the nucleation and growth of crack clusters, fault initiation and development, and discrete AE activity using a linear method.
In PFC, to simulate the heterogeneity of materials, the particle size is set within a certain range (from the minimum radius to the maximum radius). Moreover, the bonding strength between adjacent particles is usually assumed to follow a certain nonuniform distribution law, such as normal distribution. The overall mechanical properties of models are closer to the real rock by adjusting the heterogeneous distribution parameters, and the random occurrence of AE events can be also simulated [12, 13].
3.2.2. Heterogeneous Model Based on Microparticle Composition of Rock Material
In FDEM, Mahabadi et al. [14] presented a new heterogeneity modelling method considering the microscale characterization of rock. Minerals distribution, quantified microcrack density, and a series of micromechanical parameters of granite were obtained by conducting some fine tests, like grid micro-indentation, micro-scratch tests, and micro-CT scan, and this information was used as accurate input parameters and microstructure framework for numerical models. Liu et al. [15] established the heterogeneous uniaxial compression model with the above method and analyzed the evolution characteristics of AE activity of loaded granite samples further. The heterogeneity modelling method was thought to predict the mechanical response of geomaterials more accurately. Li et al. [16–18] systematically compared the advantages and limitations of different grain-mimicking methods, like UDEC-GBM (UDEC-grain boundary model), GB-FDEM (grain-based finite-discrete element method), cluster, and clump, and discussed intergranular and transgranular fracturing models and difficulties in parameter calibrations.
However, there may be some unavoidable problems in the above modelling methods based on microparticle composition of rock material. Firstly, mechanical parameters of the modelling method based on mineral distribution of rock and micromechanical tests are difficult to obtain by laboratory tests, including the micro-CT scan and nano-indentation test. Secondly, these mechanical parameters obtained by a small number of tests may not reflect the overall distribution in rock. Properties of mineral crack element, like mode I fracture energy and mode II fracture energy, need to be calibrated based on the loading curves of laboratory tests; therefore, having more mineral components in rock means more complex calibrations.
3.2.3. Heterogeneous Model with the Feature of Random Distribution Based on Microstructure of Rock Material
In order to establish a heterogeneous rock model and simulate AE activity in the process of rock macrofailure development, this study introduces the Weibull distribution which is widely used in cumulative damage statistics [19] into the open-source code, FDEM. In the algorithm, relevant mechanics properties of the rock model obey the nonuniform distribution.(1)Weibull distribution: general probability density function of two-parameters Weibull model can be shown as follows: wherein m defines the shape of the density function and n is the scale parameter of the variable. It can be assumed that and n = 1, where x and x0 represent mechanical parameters and their average values of crack element, respectively. It means that ratios of mechanical parameters to their average values obey Weibull distribution. The relationship between the probability density function and shape parameter, m, is shown in Figure 5. It can be seen that larger value of m corresponds to more concentrated distribution function and smaller value of m corresponds to more discrete distribution function. When m equals 1, the function represents the exponential distribution.(2)Weibull distribution of crack element parameters: in FDEM, the model parameters can be divided into triangular solid element parameters, quadrilateral crack element parameters, and element contact parameters according to their attributions, as shown in Table 1.
The sandstone is taken as the simulation object in this study. Figure 6 shows the polarized light microscopy observation of petrographic thin section of sandstone. It is thought that the micromechanical heterogeneity of sandstone is mainly attributed to the cementation tightness difference of clay matrix in the gap of adjacent mineral particles. Mineral particles have higher strengths and usually are intact before the clay matrix fails in loading test. Therefore, triangular solid element parameters and element contact parameters in rock model can be set homogeneously, and the values of elastic modulus (E), Poisson's ratio (ν), density (ρ), and friction angle (φ) equal the corresponding macromechanical parameters measured by laboratory test. According to Mahabadi’s study [20], the normal and tangential penalty values and which represent the contact stiffness can adopt 10 times and 1 times of the elastic modulus, respectively.
In terms of quadrilateral crack element parameters, the fracture penalty () is a coefficient and can be set the same for each crack element. The basic tensile strength (ft), cohesion (c), and internal friction coefficient () are not equivalent with macromechanical parameters of rock and can be set as the products of test values and the strength coefficient (coe). These parameters need to be multiplied by the heterogeneous coefficient (t) which follows Weibull distribution before assignment to a certain crack element. Mode I fracture energy () and mode II fracture energy () have properties of area, so they need to be multiplied by the square of distribution coefficient t. Therefore, the mechanical parameters of crack element with number i can be calculated as follows:
4. Micromechanical Parameters Determination of Rock Material
Same as the typical modelling method discrete element, several important microparameters in FDEM, such as penalty values (, , and ), fracture energy ( and ), strength coefficient (coe), and shape parameter m of density function, cannot be measured directly. Therefore, to determine these parameters, repeated calibration processes referring to existing laboratory tests are necessary before the formal simulation.
4.1. Model Description
In FDEM model, crack can only propagate along the boundary of a triangular element. In order to reduce the mesh dependence caused by the crack propagation mode, the little mesh size is favorable in spite of the increase of computational complexity. In uniaxial compression model shown in Figure 7, the sizes of 2D rock sample are set to 50 mm × 100 mm. The average edge size of Delaunay triangular finite element mesh is h = 0.8 mm. The integration time step is set as 2 × 10−9 s, and the constant axial moving velocity of two rigid platens at ends of the sample is set as 0.05 m/s. It should be noted that although the loading speed in the model is significantly higher than the actual loading speed used in the laboratory test, it ensures that the computation time of simulation process under the quasistatic loading condition is within acceptable limits [20]. Meanwhile, the viscous damping is also the key to make the nodal kinetic energy dissipate rapidly and ensure the quasistatic condition under the above loading speed.
4.2. Effect of m Value on the Macromechanical Performance and AE Response of the Model
For strength parameters of rock crack elements subject to Weibull distribution (m = 1.5, 3.0 and 6.0, resp.), corresponding stress-strain curves and AE activities of rock model under uniaxial compression are shown in Figure 8. For little m value, the mechanical properties of crack element in the model are highly discrete, local microdamage is easy to occur, AE activity increases gradually before the peak force, the overall strength of the rock model is low, and the slope of stress-strain curve decreases gradually before the peak. With the increase of m, the mechanical properties of the crack elements in the model are more uniform, local discrete damage occurrence is rare, AE activity is concentrated near the peak force, the overall strength of the rock model increases, and the slope of stress-strain curve shows a more brittle drop after the peak force.
(a)
(b)
(c)
4.3. Micromechanical Parameters of the Model
The stress-strain curve and AE response of sandstone specimen in laboratory tests are consistent with the progressive damage characteristics of rock with low m value. Further simulation comparisons show that the damage evolution characteristics of sandstone can be reflected well when the value of m is 2.0. Through several calibration simulations of mechanical parameters, a group of effective microparameters are obtained as shown in Table 2.
5. Simulation of Rock Failure Process and AE Event Evolution
5.1. Macromechanical Responses of Sandstone Failure
Simulated curve and test curve of stress to strain for sandstone in uniaxial compression test are shown in Figure 9. Figure 10 shows the real failure mode of sandstone sample after the test.
The comparison of loading curves shows that the numerical model fails to reflect the initial compaction stage because the natural pores in real rock are not considered, but its elastic modulus, Poisson’s ratio, and strength are very close to the test value after mechanical parameter calibration. The simulated lateral deformation curve of the model is close to the test curve at the beginning and then gradually deviates from the test curve due to occurrence of local failure when approaching the peak.
Figure 11 shows the axial stress-strain curve and AE response of heterogeneous model in uniaxial compression simulation. Before about sixty percent of peak force, AE activity in the model is relatively calm. With the increase of axial stress, the AE activity gradually becomes active and hits keep continuous growth. The stress curve begins to drop after its peak, and AE hits increase sharply and reach the maximum value, indicating that the model has a macrofailure at this moment. Compared with the AE activity monitored in the uniaxial compression test (as shown in Figure 12), the two AE records are well consistent in evolution law, which shows that the simulation can reflect the mechanical heterogeneity of sandstone well.
5.2. Evolution of Rock Damage and AE Activity
Figures 13 and 14 show the evolution process of rock damage and AE activity in uniaxial compression simulation. A-F corresponds to the different feature loading points in Figure 11. Each subgraph in Figure 13 illustrates the damage state and crack propagation of rock model corresponding to the current loading point wherein the blue, red, and green microcracks correspond to the tensile, shear, and mixed fracture modes, respectively. Each subgraph in Figure 14 shows the distribution of AE events accumulated in the former loading stage.
(a)
(b)
(c)
(d)
(e)
(f)
In the initial loading stage (O-A), the model remains intact, rock damages occur discretely, and AE events are mainly at low-energy level. With rising axial stress (A-B), rock damages gradually increase, the size of microcracks grows, and some local fractures propagate at the boundary of model which corresponds to a little of rock debris spalling. Meanwhile, some discrete high-energy AE events cluster in the rock, which may indicate the initiation of macrocrack. At the peak force point (C), there occurs obvious local microcrack coalescence at the lower right corner of model, and the macrocrack begins to appear. After that (D–F), with the sharp drop of loading curve, the oblique fracture belt rapidly extends and coalesces. It is characterized by the macroshear failure mode. Meanwhile, a large number of high-energy AE events occur along the macrofracture zone. Under the uniaxial compression, the inner damage in rock goes through the evolution process from discrete occurrence to concentrated rising, and the average magnitude of AE events has also changed from the low energy to the high energy. It indicates that the macrofailure of rock under loading is a progressive damage accumulation process with changing rates. The macrofailure mode of uniaxial compression model is shear along a single slope, with the shear angle of 67°, which is approximately equal to the theoretical shear angle 65.3° (=45°+φ/2).
5.3. Statistical Analysis of Rock Damage and AE Activity
On the other hand, ratios change of simulated microcracks in different modes to the total numbers and energy of microcracks with the stress-strain curve are also statistically analyzed as shown in Figure 15. During the whole process of uniaxial compression simulation, shear microcracks play the dominant role in both quantity and energy release; particularly, they sharply increase in postpeak stage. It has a satisfactory coherence with the real macrofailure mode of rock (as shown in Figure 10).
It is widely known that the AE event magnitudes against their cumulative frequency in fracture tests follow the Gutenberg-Richter relationship [21], which can be expressed as follows:where N is the cumulative number of events with magnitudes greater than M and a and b are constants. The constants a and b represent the mean AE activity level and relative proportion of events with low magnitudes, respectively. That means, the larger the b-value, the larger the proportion of low-energy events and vice versa. The relationship has been studied and confirmed by numerous researches, including laboratory tests and numerical simulations [3, 6, 9, 22, 23].
Figure 16 shows the evolution process of b-value during the present UCS test simulation. It can be observed that the seismic b-value ranges between 0.37 and 0.75, which is lower than the experimental values for granite ranging between 1.1 and 2.4 [6, 22], and the simulated values for granite ranging 0.76 and 1.80 [3]. The b-value experiences two significant drops. The first drop occurs at about 88 percent of the peak stress and the b-value decreases from 0.70 to 0.58. The second drop immediately occurs after the peak stress, and the b-value decreases from 0.53 to 0.38. The research of Lei et al. [24] indicated that the physical mechanism of b-value decreases can be attributed to the change of the dominant cracking mode from tensile mode to shear mode and more active interaction. In terms of the present simulation, the sharp increase of shear cracks number and energy ratio in peak force stage plays an important role in b-value decreases (as shown in Figure 15). The diversity of b-value range may mainly result from the property difference of simulated materials.
6. Conclusions
In this study, the Weibull random distribution of mechanical parameters of crack element is introduced into the heterogeneity modelling in the combined finite-discrete element method, in order to simulate the progressive damage process and AE activity. A set of micromechanical parameters of sandstone model under uniaxial compression are obtained through calibration simulations based on laboratory tests. Meanwhile, the development process of damage and fracture and the AE activity evolution of sandstone samples are analyzed.
It indicates that FDEM simulation with loaded heterogeneous rock model can reflect the development process of progressive rock damage, fracture, and acoustic emission activity in real rock well. Meanwhile, shear microcracks during the whole process of uniaxial compression simulation play a dominant role in both quantity and energy release, which has a satisfactory coherence with the real macrofailure mode of rock. The change of b-value also agrees with the increasing trend of high-energy events in the loading process. AE event energy based on the kinetic energy change of crack element in heterogeneous rock model provides a new tool for the analysis of rock damage and fracture evolution, and AE event magnitude can quantitatively reflect microcrack intensity.
Data Availability
The data used to support the findings of this study are included within this article.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The authors are grateful to the scholars who have contributed to this field of study. This work was supported by China Postdoctoral Science Foundation (Grant no. 2019M662568) and CRSRI Open Research Program (Grant no. CKWV2019741/KY).