At the laboratory scale, locating acoustic emission (AE) events is a comparatively mature method for evaluating cracks in rock materials, and the method plays an important role in numerical simulations. This study is aimed at developing a quantitative method for the measurement of acoustic emission (AE) events in numerical simulations. Furthermore, this method was applied to estimate the crack initiation, propagation, and coalescence in rock materials. The discrete element method-acoustic emission model (DEM-AE model) was developed using an independent subprogram. This model was designed to calculate the scalar seismic tensor of particles in the process of movement and further to determine the magnitude of AE events. An algorithm for identifying the same spatiotemporal AE event is being presented. To validate the model, a systematic physical experiment and numerical simulation for argillaceous sandstones were performed to present a quantitative comparison of the results with confining pressure. The results showed good agreement in terms of magnitude and spatiotemporal evolution between the simulation and the physical experiment. Finally, the magnitude of AE events was analyzed, and the relationship between AE events and microcracks was discussed. This model can provide the research basis for preventing seismic hazards caused by underground coal mining.

1. Introduction

An increasing number of underground engineering and physical experiments have focused on investigating the failure or fracture properties of rocks [1, 2]. An understanding of the characteristics of AE in the process of crack propagation in rock materials helps to prevent underground engineering hazards, such as rock burst [3, 4], mining tremor [5], roadway rib spalling [6], and similar events. AE signals can be triggered when distortions in rock materials occur, such as generation of microcracks [7], pore closing or collapse [8], or fault sliding [9, 10]. Traditional rock deformation detectors (e.g., strain gauges) are only able to detect the overall deformation of rocks. Along with the improvement of experimental conditions, modern laboratory techniques allow for the detection of microcracks or local macroscopic fractures [11].

Among the numerical methods that have been used to study the deformation and cracks in rock materials, the finite element method (FEM) or Fast Lagrangian Analysis of Continua (FLAC) cannot generate actual cracks, and the majority of investigations have focused mainly on studying the plastic or damaged zones [12, 13]. The DEM can generate actual cracks and identify the types of cracks [14, 15], but few studies have used this method to perform a quantitative analysis of the magnitude and location of the crack source. Although the detection of AE events during the deformation process in rock materials is a powerful tool for the quantitative analysis of cracks, the AE characteristics cannot be derived directly from numerical simulations [16], and therefore a new module using an independent subprogram is needed. Several researchers have contributed to the monitoring of AE events using numerical simulations. Based on the FEM/DEM method, Lisjak et al. [17] proposed a new model to simulate microseismic events, which described the fracture process in rock materials using three parameters: -value, AE rate, and -value. Heinze et al. [18] developed a poroelastoplastic continuum model, which divided and quantified rock failure into a prefailure stage, a massive failure stage, and a macroscopic failure stage. Through sensitivity analysis, Tang et al. [19, 20] evaluated the strength and fracture processes of rock materials and applied the Weibull distribution rule of AE events to detect rock fractures and deformations.

In DEM, if each generated microcrack is considered an individual AE event, the magnitude of the AE event caused by the hypocenter of microcrack inversion is almost the same for each event. At the laboratory scale or at the field scale, the magnitude of AE events is generally in line with the exponential distribution [17]. Therefore, most previous studies have been unable to describe the spatiotemporal distribution, magnitude, and other behaviors of AE events simultaneously. Deformation and slip of rocks can cause quick release of energy and indicate the initiation, propagation, and coalescence of cracks as well [21]. The seismic moment can be used to represent the magnitude of the seismic source [22], which can be obtained from AE monitoring. Thus, in DEM, the magnitude of the seismic source can be determined by calculating the particle moment tensor, thereby establishing a natural congruent relationship between AE events and occurrence of microcracks.

Because the stress and movement of unit particles are already known, in this study, a new module was developed incorporating an independent subprogram to calculate the moment tensor of rock materials as cracks developed, so that the magnitude of AE events could be established. An algorithm for identifying the same spatiotemporal AE event was developed. This new model was defined as the DEM-AE model. Physical experiments were conducted for a quantitative analysis of the spatiotemporal relationship of the AE events as cracks were generated in argillaceous sandstones under different confining pressures. A validation analysis was performed on the model. A comparison of the magnitude and spatiotemporal relationship of AE events in the numerical simulations and the physical experiments demonstrated that the simulations were in good agreement with the experimental data. Finally, in this study, we also analyzed the magnitude of AE events under different confining pressures in numerical simulations and discussed the relationship between AE events and the number of microcracks.

2. Simulation Mechanism for Using PFC

A common representation of an explicit DEM is the Particle Flow Code (PFC). PFC theory assumes that macroscopic rock materials are composed of many bonded microscopic particles, which is a mechanical behavior of representing macroscopic rock materials by the relative motion of microscopic particles [23, 24]. PFC has the following advantages: () it overcomes the theoretical assumption of traditional continuum-based models [14] and () PFC is more suitable for simulating nonlinear deformation in rock materials and for detecting the number and evolution of cracks [25]. Two bond models are proposed based on the type of bond between particles. Based on previous work, the parallel bond model (PBM) was deemed more suitable for simulation of rock materials [26, 27]. A simulation of the basic mechanical behaviors of rock materials through PBM is shown in Figure 1.

2.1. Governing Equation of Crack Generation

In PBM, the movement of particles causes changes in force and moment. The resultant force and moment can be divided into normal- and shear-directed components, which can be expressed as follows:where , , and are the total, normal-, and shear-directed force of the parallel bond, respectively; , , and are the total, normal-, and shear-directed moment of the parallel bond, respectively; and are the normal- and shear-directed unit vectors, respectively.

and are initialized as 0. Each subsequent related displacement and rotation will create an increment in elastic force and moment that is added to the current value. The increment in elastic force and moment is calculated according towhere and are the increments of the normal- and shear-directed force, respectively; and are the contact normal stiffness and shear stiffness, respectively; and are the increments of the elastic force that is added to the current value, respectively; and are the increments of the normal- and shear-directed moment, respectively; and are the increments of the rotational angle that is added to the current value, respectively.

Other parameters , , and are the area, moment of inertia, and polar moment of inertia of the parallel bond cross section, respectively, which are determined based on

In PBM, the maximum tension and shear stress between particles are subject to

The radius of the cross section is described bywhere is a radius multiplier (usually equal to 1.0), and are the radii of A and B, respectively (Figure 5), and and are the tensile strength and shear strength, respectively.

In the normal and shear direction, the relationship between stress and displacement of the contact points depends on and , respectively [28]. According to Figure 2, the contact behavior of normal stress and displacement can be described as follows:

If , then . At this point, tensile microcracks develop, marked with a red short line. is the residual tensile strength.

The contact behavior of shear stress can be described as follows:where is the shear displacement, and are the cohesive strength and friction angle, respectively, and and are the elastic and total component of the incremental shear displacement, respectively.

If , then ; or else, if , then . The “sign” function here is a mathematical symbol that allows to be positive. At this point, shear microcracks develop, marked with a blue short line.

2.2. Algorithm of AE Events in DEM
2.2.1. Principle of Calculating the Magnitude of AE Events

Due to the growth of new cracks, friction motion in rock materials will create AE signals. At the laboratory scale, in order to estimate the damage zone and the crack location of rocks, signal monitoring is generally carried out by using an AE instrument [29, 30]. At the field scale, a microseismograph is commonly used for monitoring the fracture zone affected by the excavation of underground chambers as well as rock burst [31]. As for monitoring of the seismic source, AE events and microseismic events have similar mechanisms but a different range of frequencies and both events often precede or accompany cracks during the detection process [16, 32]. In order to conduct a quantitative analysis of the seismic source, namely, to detect the magnitude and location of AE events, the seismic moment is often utilized as the basic approach.

The scalar seismic moment is defined aswhere is the scalar seismic moment; is the shear modulus; is the current displacement in the calculation.

In DEM, since the stress and movement of particles can be derived directly from calculations, it is easier to determine the scalar seismic moment based on the change of contact force between particles as new cracks develop. As shown in Figure 3, particles at both sides of the new microcrack are defined as source particles (Particles A and B). After the growth of new microcracks, the contact force changes due to the movement of source particles. If the event is composed of a single microcrack, its zone of action centers on the microcrack, and the radius of the affected zone is the diameter of the bigger source particle (see Figure 5). If the AE event is composed of several microcracks, the geometric center of all the microcracks will be the spatial position of the AE event. In this way, the seismic moment tensor in DEM is the summation of all the variables and their quantity of contact forces multiplied by the corresponding arm of the force (the distance between the location of the contact point and the center of the microcrack): where is the scalar seismic moment in the calculation; and are the component of the contact force and the corresponding arm of the force.

The maximum scalar moment of the moment tensor is

The magnitude of AE events is calculated as per the empirical equationwhere is the magnitude of AE events.

2.2.2. Algorithm for Identifying Space and Time of AE Events

Only when new cracks occur at other contact points of the source particles (Particles A and B) are they likely to be the same AE event, because the stress in DEM can be transferred to other adjacent particles. The zone affected by the AE event will also extend gradually along with the coalescence of cracks.

The definition of excitation time of an AE event is shown in Figure 4. Since the propagation velocity of the crack is two times as fast as the shear wave velocity of rocks [33], the excitation time and affected area of the unit particles can be calculated. Assuming that the excitation time of Crack 1 for spreading to the boundary of its affected area is (see Figure 5), the calculation of the moment tensor is updated for each time step within the period of . If no new cracks occur within or within the affected area of Crack 1, this AE event contains only one microcrack, and the total excitation time is . On the contrary, if Crack 2 (its excitation time for spreading to the boundary of the affected area is ) develops both within and within the affected area of Crack 1, the affected area of this AE event will be superposed on the affected area of Crack 2, and the total excitation time is . By analogy, if Crack 3 comes into being within , the affected area and magnitude will also be superposed. As for this AE event in Figures 4 and 5, its affected area is the superposition of the action ranges of Cracks 1, 2, and 3, and its total excitation time is .

As shown in Figures 4 and 5, even if Crack 4 comes into being within the affected area of Crack 1 in terms of space, the excitation time for the development of Crack 4 is after , so Crack 4 and Crack 1 cannot be regarded as the same AE event. The excitation time when Crack 5 occurs is within , but the location where Crack 5 occurs is not within the affected area of Crack 1, so Crack 5 and Crack 1 cannot be regarded as the same AE event, either.

3. Laboratory Tests

3.1. Experimental Procedure

The rock materials used in this study were argillaceous sandstones from the Shanxi Formation at Shanxi Province, China, created to  mm standard specimens with the core drilled as per the International Society for Rock Mechanics (ISRM) standard (see Figure 6). The specimens were divided into four groups for the confined compression pressure tests of 0.5, 2.0, 4.0, and 6.0 MPa. Both the upper face and the lower face of the specimens were carefully polished and flatted.

As illustrated in Figure 7, the physical compression test with confining pressure was performed on a universal servo-control testing machine (MTS815.02, America) by imposing a constant speed of 0.2 mm/min until the ultimate failure of the specimens. The AE system was utilized to detect the magnitude and location of AE events. In this study, the AE system was composed of six commercially available piezoelectric sensors (Physical Acoustic Corporation type-Nano30), preamplifiers, and a MISTRAS Micro-II Digital AE data acquisition system. The AE sensors were attached to the surface of the hydraulic cylinder using a coupling agent and were further fixed to record AE events data of the specimens constantly. If the coordinates of specimens were as shown in Figure 7, the sensor coordinates were (25, 0, 20), (−25, 0, 20), (12.5, −21.65, 50), (−12.5, 21.65, 50), (12.5, −21.65, 80), and (−12.5, 21.65, 80) with the unit in mm. The AE record threshold was set to 40 dB in order to avoid the possibility of environmental disturbance. Both the servo-control testing machine and the AE system were executed simultaneously to obtain the correlation of mechanical behavior and AE signals of specimens.

Table 1 lists the experimental results for all the specimens. As for each group of specimens, the lateral pressure and corresponding peak deviatoric stress were obtained by linear regression, as shown in Figure 8. The gradient and intercept of the -axis were used for the calculation of cohesion c and internal friction angle , with the calculation equations as follows [34]:

According to Table 1 and Figure 8, Young’s modulus showed almost no change in increment as the confining pressure increased from 0.5 MPa to 6 MPa, and the average Young modulus was 21.34 GPa. Regarding the calculation results of (11) and (12), the internal friction angle was 19.36°, and the cohesion was 6.13 MPa.

3.2. Calibration of Micromechanical Parameters

The discrete element model must be calibrated to the associated microscopic parameters in order to describe the macroscopic mechanical behaviors of rock materials (i.e., deviatoric stress and Young’s modulus). Although the calculation equations for parameters such as cohesion, stiffness, and friction angle are included in the manual of the PFC software [35], there is no fixed theory that validates the microscopic parameters of rock materials [27]. Therefore, in order to validate the model, the following steps were adopted for parameter calibration:(1)Young’s modulus and strength parameters (i.e., cohesion and internal friction angle) were calibrated first. The calibration standard was subjected to the cohesion and the internal friction angle obtained from physical experiments, rather than using only the uniaxial compressive strength of argillaceous sandstones as the basis.(2)Initial microscopic properties were assigned in accordance with the relationship between microscopic parameters and macroscopic mechanical behaviors of rock materials.(3)In the DEM, a servo system was used for compression tests with the confining pressure set at 0.5, 2.0, 4.0, and 6.0 MPa, respectively. In this case, trial and error was applied to determine microscopic parameters, until Young’s modulus reached 21.34 GPa and the peak deviatoric stress was roughly consistent with the envelope line resulting from physical experiments. Microscopic parameters of argillaceous sandstones after calibration are shown in Table 2, and the deviatoric stress-strain curve of different confining pressures is displayed in Figure 8.(4)The established DEM-AE model was compared with the results of the physical experiment (mainly including the spatiotemporal location of AE events and magnitude), and the comparison results are presented in Section 4.

4. Results and Discussion

4.1. Time Distribution of AE Events

After ensuring that the macroscopic parameters of the argillaceous sandstone specimens in the physical experiments were consistent with those of the numerical simulations, the contrastive analysis of the DEM-AE model was carried out. Figure 9 shows the comparison of AE events for the physical experiments and the DEM-AE model, in which the -axis represents the relative axial strain and the -axis represents the ratio of cumulative AE events.

Under the different confining pressures, the results showed good agreement for the ratio of cumulative AE events in the physical experiments and in the DEM-AE model. It is worth noting that, in the physical experiment, cumulative AE events demonstrated linear growth at the initial loading stage; therefore, AE events could be detected even at this loading stage in the physical experiment, causing deviation from the numerical simulation. The detailed reasons are described in Section 4.4.

With the different confining pressures, the number of AE events prior to loading to the peak deviatoric stress was very small, especially in the numerical simulation. There were nearly no AE events at the linear elastic stage, followed by a small number of microcracks in the specimens. In the DEM-AE model, only a small number of AE events were monitored as well. Most AE events developed after the peak and cumulative failure events increased exponentially at the postpeak stage. In particular, as for specimens with lower confining pressure ( and 2.0 MPa), the increase of cumulative AE events did not stop and start again after the peak deviatoric stress, while such a phenomenon was observed for the specimens with higher confining pressure ( and 6.0 MPa). Additionally, the higher the confining pressure was, the more the intermittent stops occurred because a transitory secondary resistance phenomenon took place after the simulated specimens were damaged due to the high confining pressure. However, the postpeak curves in the physical experiment had better brittle-to-ductile transition; therefore, no intermittent stops of AE events occurred.

4.2. Spatial Distribution of AE Events

Based on the microscopic parameters after calibration, the failure modes of the argillaceous sandstones under different confining pressures can be seen in Figure 10. The DEM-AE model showed good agreement with the failure modes derived from the physical experiments. With the increase in confining pressure, the residual resistance of rock specimens increased, leading to a gradual increase in cracks (see Figure 11). It is noteworthy that tensile cracks always played a dominant role and were 2~2.5 times as abundant as shear cracks. The number and magnitude of AE events were also on the rise under the different confining pressures.

The temporal evolution of AE events at the stages of 50% peak deviatoric stress, peak deviatoric stress, and postpeak stage is shown in Figures 1215. At the point of 50% peak deviatoric stress, there were a small number of AE events in the physical experiments, about 7%~10% of the cumulative total of AE events, while in the numerical simulations, there were almost no AE events prior to the 50% peak deviatoric stress.

At the peak deviatoric stress, the AE events in the physical experiments and the numerical simulations accounted for 20%~35% of the cumulative total, while no obvious macroscopic cracks occurred in the specimens during the loading procedure. Until loading to 90% of postpeak strength, the physical experiments and numerical simulations resulted in relatively obvious cracks among the specimens with the confining pressures of 0.5, 2.0, and 4.0 MPa, accounting for 40%~55% of the cumulative total of AE events. It is worth noting that there were obvious macroscopic cracks in specimens with the confining pressure of 6.0 MPa only until 80% of postpeak strength.

4.3. Magnitude of AE Events

During monitoring of the spatial distribution of AE events, the magnitude of AE events can be calculated in the DEM-AE model as well. As shown in Figure 16, when the magnitude of AE events was smaller than −6.0 or larger than −4.5, there was not much difference in the number of AE events under four different confining pressures. Most of the AE events took place in the range −6.0~−4.5 and were on the increase along with the increase in the confining pressure. Based on the magnitude of AE events in Figure 16, it can be seen that AE events conformed to a normal distribution, and the peaks were within −5.25~−5.0 (as shown in Figures 16(a)16(d)). Additionally, with the increase in confining pressure, the concentration area of peaks tended to shift to the larger position.

Figure 17 shows the relationship between the AE events and the number of microcracks. When the confining pressure was 0.5 MPa, the number of AE events containing only one microcrack accounted for over 90% of the total number, which decreased gradually along with the increase of confining pressure, but the total number remained over 85%. As for different confining pressures, the number of microcracks showed an exponential decrease along with the increase of AE events. Under the confining pressure of 0.5 MPa, a maximum of eight microcracks were generated by only one AE event. When the confining pressure reached 6.0 MPa, the maximum number of microcracks generated by one AE event increased to 15, which occurred only once.

4.4. Discussion

As described in Section 4.1, under the different confining pressures, the postpeak stage in the physical experiments showed better brittle-to-ductile transition; however, this phenomenon was not observed in the numerical simulations (see Figure 8). As the numerical simulation was loaded to peak deviatoric stress, massive cracks were coalescing and interlocking among most of the particles in the simulated specimens, leading to an abrupt drop in peak deviatoric stress. This progression was also demonstrated in the AE model so that more than 70%~80% of the AE events were observed only after the peak deviatoric stress. In the physical experiments, however, hard minerals that exist in argillaceous sandstones were able to resist slippage and instability at the postpeak stage; therefore, the postpeak stage showed better brittle-to-ductile transition.

In addition, because of the closing of initial cracks or holes in the physical experiments, the AE events of specimens could be monitored even at the initial loading stage, which was not observed in the numerical simulations. Nevertheless, AE events at the initial loading stage accounted for only a small proportion of the entire loading process.

5. Conclusions

In this paper, a DEM-AE model was presented based on the movement of particles during the loading process. The fundamental principle of the process was the calculation of the scalar moment of the moment tensor and the determination of the magnitude of AE events. An algorithm for identifying the spatiotemporal location of the same AE event was put forward. With this method, it was possible to analyze and quantify the spatiotemporal location and the magnitude during the process of crack initiation, propagation, and coalescence. This model was applied to the compression test with different confining pressures of argillaceous sandstones. A systematic contrastive analysis was conducted with the results from the physical experiment and the numerical simulation, leading to a good agreement between both cases in terms of spatiotemporal distribution as well as the magnitude of AE events. Furthermore, this model was used to analyze the magnitude of AE events and discuss the relationship between AE events and microcracks.

With this model, an authentic simulation of AE events in physical experiments can be reproduced using numerical simulations. In the near future, this model will be applied to the production process at the field scale. The performance of this model makes it possible to evaluate and prevent seismic hazards caused by roadway excavation and coal mining.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


Financial support for this work is provided by the National Natural Science Foundation of China (no. 51474208), the National Key Research and Development Program of China (2016YFC0600904), and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). The financial support provided by China Scholarship Council (CSC, Grant no. 201606420013) is also gratefully acknowledged.