Abstract
In view of the randomness of rockfalls shape and irregularity of the bottom floor of working face in steeply dipping coal seams (SDCS), it is difficult to accurately simulate rockfall movement, and it is consequently unable to effectively protect against multirockfalls. Therefore, a method for generating random shape rockfalls based on ellipsoid equation is proposed, and a 3D grid model of real bottom floor of working face is established based on the geographic information system data. In order to verify the accuracy and feasibility of the method and 3D model, the trajectory simulated by Rockyfor3D software is compared, and the proposed method and 3D model prove to be effective in simulating rockfall movement more accurately. Then the proposed method and 3D grid model are applied to solve the problem of multirockfalls protection in numerical simulation, and the main factors affecting the structural stress response of protective netting are analyzed by taking three incident modes of parallel heights, ladder parallel, and the same trajectory. In the simulation, it is found out that the trajectory of irregular rockfalls is greatly affected by the shape of rockfall and working face floor; during the process of multiple rockfalls colliding with the protective netting, the peak stress on the protective netting is inversely proportional to both the time interval between each rockfall and the distance between each rockfall. The findings presented in this research contribute to rockfall prediction and protection against rockfall hazards.
1. Introduction
In recent years, with the development of SDCS mining and combining mining technology, the number of steep coal mines in China has increased year by year, and therefore the coal production, as well as the economic benefits, has increased [1–7]. At the same time, the rockfall hazard, a disaster that occurred in SDCS mining, has become increasingly prominent [8–19]. Rockfalls in the long walls of a SDCS also have a colloquial name, “flying gangue” hazards [20]. Because of the large difference in height between the lower end and the upper end of the mining working face, the kinetic energy of a coal/rock block falling off from the parent body would be significantly increased after acceleration; the collision of the coal/rock blocks with people or equipment would occur instantly, and thus it features complicated motion and may lead to serious accidents. During mining, coal/rock blocks may be separated from coal walls or rest on the top beam of a support or stay still after the bottom plate is slipped and destroyed. The coal/rock blocks, once subjected to external disturbance such as coal swing or support moving, will start to move from static state. Since the natural-rest-angle of a coal/rock block is smaller than the coal seam inclination angle, it is difficult to stop after the first collision with the bottom plate and will slide downwards along the working face and may collide with people or equipment. The coal seam structural information and the main movement modes of rockfall are shown in Figure 1 [12]. As can be seen from Figure 1, the coal seam is on the left and goaf is on the right. In addition, Figure 1(a) shows the primary collision slip mode of rockfall, and Figure 1(b) shows the multiple intermittent collision and slip mode of rockfall.

(a)

(b)
Therefore, studies on the mechanism of rockfalls motion would fundamentally improve the protection against rockfalls that occurred in steep seam mining. While the research on steep seam mining mainly focuses on mining equipment, especially hydraulic supports and shearers [21–30], and technical management of longwall working faces, only few studies have been carried on the motion mechanism of rockfalls in steep seam mining [9, 10]. In 2015, Tu et al. comprehensively review the fully mechanized coal mining technology in steep coal seam in China. Reference [8] and different rockfalls prevention devices applied in practice are compared. In 2017, a rockfall control principle is proposed by Wu et al. [11]: depending on the occurrence locations of rockfalls, rockfalls can be dealt with through technical management. In 2018, based on the Kalman filter principle, the parameters of rockfall control are studied and the protective netting is proposed and applied in practice by Wu et al. [12, 13]. In 2019 and 2020, based on nonprobability interval analysis, the motion model of rockfalls is established under uncertain environment. With theoretical analysis and numerical simulation, the kinematic characteristics of rockfalls in steep seam mining on longwall mining working face are investigated, and the kinematic mechanism of rockfalls is demonstrated by Jing et al. [14]. Based on the grating to identify the dangerous source, automatic devices for preventing rockfalls are designed by Liu et al. [15–17]. The index scale method and analytic hierarchy method are used to calculate the weight, the expert scoring method is used to determine the membership function, and the safety evaluation of rockfall hazards is carried out by the fuzzy comprehensive evaluation method by Liu et al. [18]. And the dynamic Bayesian network model is used to evaluate the dynamic threat level for rockfalls along the working face by Liu et al. [19]. In 2021, a method for risk assessment and for determining the principles of protective systems is provided in underground steep coal seams by Wu et al. [20]. Hu and Luo do a case study about rockfall hazards and failure characteristics in steeply dipping coal seam [21, 22].
Although the above devices and technological managements have improved rockfalls prevention to a certain extent, they still cannot fundamentally enhance the protection against the rockfalls, without considering the kinematic characteristics of rockfalls. Therefore, it is of great significance to explore the mechanism of rockfalls movement for establishing the safety protection system of longwall working face in SDCS mining area.
Most of the rockfall and working face floor models used in the above research are idealized models and in lack of research on three-dimensional movement of irregular random rockfall. The innovation of this research is: the previous two-dimensional space rockfall trajectory simulation is upgraded to three-dimensional stepping space rockfall trajectory simulation; the former regular shape rockfall modeling is converted into irregular random shape in model establishment; the idealized simple working face modeling is changed to the working face modeling based on geographic information system. The innovation of the current work contributes to the prediction of rockfall movement and therefore the protection against rockfall disasters. The research is given as follows: the modeling method of irregular rockfall and the modeling method based on geographic information system are put forward; according to the proposed principles, a group of simulation tests of rockfall movement are carried out to find out the most dangerous rockfall model; the protection simulation test is carried out for the most dangerous rockfall model.
2. Basic Principles
2.1. Energy Tracking Method
Energy tracking method (ETM) is based on the principle of energy iteration. The collision between irregular rockfalls and working face floor and the collision between rockfalls are mainly multipoint collisions, as shown in Figure 2 [31, 32].

There are a lot of blocks in the figure with multipoint collision. is the collision point, n is the number of all collision points, and m is the number of all collision bodies.
ETM uses the “Stronge assumption” [33], to express the energy dissipation in the normal direction. Let be the work done by the normal component of the impact impulse in a collision. The relational expression can be as follows:where is normal impact recovery coefficient, and is work done by normal impulse at the maximum compression point where the relative normal velocity changes.
According to the normal component and tangential component of collision, the local orthogonal coordinates system is defined as at the contact, as shown in Figure 3.

The relational expression can be as follows:where is normal vector of contact point, and and are tangent vectors of contact points, respectively.
Mirtich proposes that the work of the contact force is a function of the relative velocity before and after the impulse is applied [31], which can be expressed as
In ETM, the change of relative velocity along the normal direction is first calculated, and the change of relative normal velocity is
The normal and tangential components decomposed by impulse can be expressed as
Define the condition of static friction aswhere is friction factor.
To ensure that the relative change of normal velocity is , the impulse must be calculated as follows:where .
2.2. Modeling of Irregular Random Shape Rockfalls
Corresponding to the problem, randomness of rockfall shape, based on the model shown in Figure 3, an innovative model is established which can realize the randomness of rockfall shape.
The ellipsoid equation is used to transform into an irregular random rockfall model. Establishing a space coordinate system, the coordinates of any point in the space are represented by r, , and . Among them, r is the distance from the point to the center of the sphere, is the latitude of the point, and is the longitude of the point.
First determine the and of each vertex of the rockfall. Take two random variables and as intermediate variables, which can be determined by the following formula:where , , and are random variables on , , that is, to control variation amplitude of and , , , the number of vertices of the polyhedron is , and among them .
Then and can be expressed as
Adding the two points and , any vertex on the polyhedron can be expressed aswhere , , is a random variable on , is the value that controls the magnitude of change in r.
For a polyhedral rockfall, only the parameters , , , , γ, and need to be given. Both and are random numbers on , which control the randomness of the rockfall shape. Use the two parameters and to control the overall shape of the rockfall and then control the sharpness of the rockfall. When and are close and small, the shape of the rockfall is relatively slender. When is small and is close to 1, the shape of the rockfall is relatively flat. When and are close, and both are close to 1, the rockfall shape is close to spherical. With this method, using the model established in this section, the rockfall with arbitrary sphericity and roundness of polyhedral random shapes can be simulated. The maximum diameter of the rockfall model is the longest diagonal length r of the polyhedron, and its value is usually a fixed value.
2.3. Modeling Base on Geographic Information System
This article takes a test mine with SDCS in Gansu Province as an example. Contour lines and drilling information of the working face floor of the mine can be shown in Figure 4. Figure 4 is a horizontal projection drawing of contour lines of a super thick coal seam floor in Gansu province. The curve lines in the drawing are contour lines, and the points are drilling positions. The original drawing is CAD engineering information map, which has coded geographic information of coal seam working face and can be used to establish a three-dimensional model of the working face after extraction.

Extract contour and borehole information from Figure 4. Using CASS software, the drilling information is regarded as an elevation point, and a mesh model of the working face bottom plate is established, as shown in Figure 5.

The model is a 3D mesh model of 4500 m × 2500 m × 500 m. It covers the geographic information of the entire coal seam floor, the inclination of the coal seam is about 25°–35°, and the upper limit belongs to the SDCS. To select the location of the working face bottom plate that is most prone to rockfall hazards for numerical simulation experiments, Arcgis and Rockyfor3D software are used for model rasterization and face inclination analysis, and the results are shown in Figures 6 and 7, respectively.


It can be seen from Figures 6 and 7 that the purple area (marked in Figure 7) is the area that is most prone to rockfall hazards; the inclination angle of the working surface in this area is about 35° which belongs to the SDCS. The detailed bottom floor of working face of the area is presented in a 3D model, as shown in Figure 8.

2.4. Numerical Simulation Principle of Rockfalls Impact Protection Netting
The problem of rockfalls collision protection structure belongs to the study of structural dynamics. This article adopts the explicit analysis algorithm in HyperMesh/LS-DYNA software. Numerical simulation of the process of rockfalls collision protection structure is analyzed. The structural dynamics equation iswhere , , and are the mass, damping, and stiffness matrices of the protective netting structure system, respectively, is the speed response of the falling rock collision protection netting, and is the load array of the protective netting.
The central difference method is an explicit algorithm commonly used in this field. During the operation, the time can be divided into multiple steps. The acceleration at time iswhere is the external force of the protective netting at time , and is the internal force of the protective netting.
Let ; the velocity and displacement at time can be obtained by the following formula:
To ensure the stability of this method, a small time step is required and needs to meet the following formula:where is the maximum natural vibration frequency of the protective netting, and is the maximum time step that satisfies the stability condition.
3. Verification of ETM Model
This article takes a hexahedral rockfall as an example. Verify the multipoint collision between the rockfall of the polyhedron and the bottom of the working face. The bottom plate model of the working face is shown in Figure 8. The collision recovery coefficient of the working floor is 0.84, the friction coefficient is 0.4, the rockfall density is 2500 kg/m3, and the model size is 0.1 m × 0.2 m × 0.2 m. The coal/rock block falls from the roof to form rockfall. The rockfall moves from a height of 20 m above the end of the working face and collides with the bottom of the working face in free fall. Then bounce and roll motions occur along the bottom of the working face. The trajectories simulated by ETM and the Rockyfor3D software were compared. The motion trajectory and lateral offset trajectory of rockfall along the working face inclination are obtained as in Figures 9 and 10, respectively. It can be seen that the rockfall tendency trajectory and lateral movement trajectory simulated by ETM are consistent with Rockyfor3D simulation results.


4. Simulation Results
4.1. Movement of Irregular Random Shape Rockfall
The shape of the rockfall can be expressed by sphericity and roundness, as shown in Figure 11 [34].

The sphericity calculation formula iswhere S is the sphericity of the rockfall, and A, B, C, respectively, represent the length of the three axes of the rockfall.
The roundness calculation formula iswhere P is the roundness of the rockfall, represents the radius of the inscribed circle of the corner of the rockfall model, N represents the number of corners, and R represents the maximum inscribed circle radius of the rockfall model.
Field measurement found that most of the rockfalls present irregular random shapes, and the roundness of rockfalls is mainly concentrated in , and the sphericity is mainly concentrated in . Rockfall with extremely small roundness is not easy to move. After being formed, it mostly stays on the bottom of the working face and is not easy to slip off. And neither the sphericity nor the roundness will show extremely large values. Therefore, the irregular random shape rockfall production method in this paper is used to generate the roundness in . Numerical simulation test is performed on nine rockfalls with sphericity in , as shown in Table 1, where r = 0.2 m.
The working face model in this paper is taken from an SDCS test mine in Gansu Province of China. The bottom plate model of the working face shown in Figure 9 is adopted. It is assumed that, during the mining process of the working face, the coal cutting trajectory is consistent with the contour line of the floor. The collision recovery coefficient of the working face bottom plate is 0.84, and the friction coefficient is 0.4. The density of rockfall is 2500 kg/m3, and the rockfall moves freely from 10 m above the working face. In order to study the movement law of rockfall, it is assumed that there are no equipment and other obstacles at the lower end of the working face.
ETM is used for numerical simulation experiments. All rockfall models fall in the same location. A comparison chart of rockfall trajectories as shown in Figure 12 is obtained. It can be seen from Figure 12 that the greater the sphericity of the rockfall, the stronger the bounce ability after a collision and the easier it is for a large-scale bounce motion to form a simpler parabolic motion trajectory. Rockfall with small sphericity will bounce many times on the working surface. Therefore, the uncertainty of the movement trajectory is increased. In engineering practice, rockfall with small sphericity is more unstable on the working face, and it is easy to slip with vibration. When the sphericity of rockfall is very small, the motion of rockfall will show sliding and rolling motion that fits the floor of working face, which leads to less intense movement. In the case of the same sphericity, the smaller the roundness, the higher the complexity of the rockfall motion trajectory, and the more difficult it is to protect.

(a)

(b)

(c)
The change of rockfall velocity is shown in Figure 13. It can be seen from Figure 13 that the larger the sphericity of the rockfall, the greater the initial velocity of the first bounce after the free fall, and the less the energy loss caused by collision. During the entire movement process, the peak velocity that can be reached is greater. In the case of the same sphericity, the smaller the roundness, the smaller the velocity change of rockfall after colliding with the working face floor. Its ability to convert kinetic energy into potential energy is lower, and the movement takes place longer on the working surface. Therefore, it should be focused on.

(a)

(b)

(c)
The energy changes of three groups of irregular rockfalls are shown in Figure 14. It can be seen from Figure 14 that the total kinetic energy of rockfall is composed of translational kinetic energy and rotational kinetic energy. The changing trend of the total kinetic energy of rockfall is similar to that of the center of mass velocity. It can be seen that the translational kinetic energy of rockfall accounts for a larger proportion of the total kinetic energy during the movement process. The roundness and sphericity of rockfall mainly affect the quality of rockfall by changing the shape of rockfall. And then it affects the total kinetic energy in the process of rockfall movement. The change law of total kinetic energy is affected by the shape of the bottom of the working face and the shape of the rockfall. The complexity of the field situation directly affects the occurrence of rockfall disaster. Therefore, as important technical means in the simulation process, the generation methods of irregular random shape rockfalls and working face can make the numerical simulation test effectively relate to the actual situation.

(a)

(b)

(c)
4.2. Process of Rockfall Collision Protection Netting
In the field practice, the setting of the protective netting in SDCS depends on the distribution of rockfall trajectory and the area that may be endangered. This article will adopt the principle of graded protection based on field experience. The protective netting is set at the position where the bounce height is small and the kinetic energy is small. Using HyperMesh combined with LS-DYNA software, the process of rockfalls collision protection netting is simulated numerically. Take the second set of No. 1 rockfall model as an example: set the protective netting between the second and third collisions between the rockfall and the working face floor, which is 40 m away from the lower end of the working surface. The direction is perpendicular to the floor of the working face. The rockfall enters the protective netting at a speed of 15.4 m/s and 47.33° with the horizontal plane of the protective netting. In order to explore the collision between many rockfalls and the protective netting, taking two falling rocks A and B as an example, the following three falling rocks incidence modes are used to simulate: parallel height, parallel steps, and the same trajectory. The shape and parameters of the protective netting refer to [13]. Under the three different incident modes, the rockfalls touch the netting, drop to the lowest point, and then contact and separate from the netting at three moments, as shown in Figure 15. It can be seen from Figure 15 that, with the change of time, the mesh of the middle protective netting is constantly stretched under tension. The mesh continues to slide, and the degree of tension becomes greater and greater, and the middle part of the protective netting gradually begins to tighten. When the middle part of the protective netting is fully tightened, the rockfall starts to rebound. The motion of rockfall and protective netting depends on the posture and contact position of those and the incident speed of rockfall.

(a)

(b)

(c)
As shown in Figures 16 and 17, these are curve of stress distribution and limit stress changing with time of two rockfalls A and B under three different incident modes. It can be seen from Figure 16 that, in the three modes, the connection between the mesh surface and the constraint and the mesh surface corners are the locations where the stress is the greatest. It can be seen from Figure 17 that the peak stress of the protective netting is related to the continuity of rockfall incidence. When the rockfall is parallel and incident at the same time, the stress peak value of the protective net is the largest. When two rockfalls are incident one after the other in a step parallel or in the same orbit and in the same direction, the peak stress of the protective netting will be reduced, and it is inversely proportional to the time interval between two rockfalls. In the actual protection, we should pay more attention to the impact of multiple rockfalls on the protective netting.

(a)

(b)

(c)

Figure 18 shows the velocity curve of rockfalls during the collision. It can be seen from Figure 18 that the gangue blocking effect of the protective netting is affected by the incident position and interval time of the rockfalls. When two rockfalls enter the protective netting in a parallel step pattern, the touch screen points are far away and the touch screen has a certain interval time. The falling speeds of the two rockfalls when the rebound is separated are 3.17 m/s and 3.89 m/s, respectively; it is the minimum separation speed in the three modes. Therefore, the protective netting has the most obvious protective effect against rockfalls incident in this mode. In addition, in the three modes, the protective netting has the best protective effect on the rockfall which first touches the netting, and the results are the same as the actual engineering results.

(a)

(b)

(c)
5. Conclusion
(1)In the mining of SDCS, the movement of rockfalls is a complex dynamic process. This paper proposes a method for generating rockfalls with irregular random shapes. The rockfall models with arbitrary sphericity and roundness based on ellipsoid can be generated. Combining with ETM, the movement process of rockfalls with irregular random shape in three-dimensional space can be accurately simulated, including the speed and energy changes at any moment. The simulated trajectory is highly consistent with the trajectory simulated by Rockyfor3D software.(2)The trajectory of irregular rockfalls is greatly affected by the shape of rockfall and working face floor. The bigger the sphericity of rockfall is, the stronger the jumping ability of rockfall is, and the smaller the change of velocity in the collision process is. In addition, the smaller the roundness is, the more complex the movement trajectory is. During the collision process between rockfall and protective netting, the joint between mesh and constraint and the corner of mesh are stressed the most, and fracture is most likely to occur. During the process of multiple rockfalls colliding with the protective netting, the peak stress on the protective netting is inversely proportional to the time interval between each rockfall and the distance to the netting.(3)This paper uses the contours of the working face bottom plate and drilling information for modeling, combines with ETM, simulates the three-dimensional movement process of random-shaped rockfalls, and opens up the interface between the geographic information system data model and ETM self-programming. The perfect combination of three-dimensional visualization modeling of the real working face floor of the SDCS and the self-programming of ETM is realized. The matching degree between numerical simulation test and engineering practice is improved, which has an important theoretical and practical significance for predicting rockfall hazards in the working face.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 National Natural Science Foundation of China (NNSFC) through Grant nos. 51604213 and 51634007 and the Scientific Research Project of the Department of Education of Liaoning Province (Grant no. L2020041).