The principle of ETM is energy iteration. According to the elastic energy absorbed by the collision point of rockfall at the moment of collision, a pair of impulses is applied to the collision point to release energy until the migration stops continuously. As a result, the collision between the irregular rockfall and the working face floor, as well as the collision between the rockfall, is mainly a multi-point collision, as shown in Figure 2 (Tang et al., 2014).

Figure 2 shows the multi-point collision of multiple blocks b = b 0 , b 1 , b 2 , ⋯ , b m . c = c 0 , c 1 , c 2 , ⋯ , c n is the collision point.

In each iteration step, multi-point collisions can be regarded as simultaneous collisions, and the calculated impulses can be applied to the contact points. The interaction of multi-point collisions is simulated by cyclic operation until all relative average velocities are non-negative.

ETM adopts the Stronge assumption to represent the energy dissipation in the normal direction. The work W r done by the normal component of the collision impulse in a collision can be expressed as (Bi et al., 2016): W r = − e n 2 W m c (1) where e n is the normal impact recovery coefficient; W m c is the work done by the normal impulse at the maximum compression point where the relative normal velocity changes.

According to the normal and tangential components of collision, the local coordinates system t i , n i , q i is defined at the contact, as shown in Figure 3.

The following expression can be obtained: t i = n i × p i s k u i s k (2) q i = t i × n i (3) Where n i is the normal vector, t i and q i are tangent vectors.

The work Δ W i s k of the contact force is a function of the relative velocity before and after the impulse is applied: Δ W i ( s k ) = 1 2 ( u i ( s k + 1 ) + u i ( s k ) ) T K i − 1 ( u i ( s k + 1 ) − u i ( s k ) ) (4)

In the ETM method, the change of relative velocity along the normal direction is calculated first, and then the friction is considered. The change of relative normal velocity is: u ~ i s k + 1 = 0 u ~ i , n s k + Δ u ~ i , n s k 0 T (5)

The impulse p i s k can be decomposed into normal and tangential components as: p i , n s k = n i T p i s k n i (6) p i , q s k = p i s k − p i , n s k (7)

The condition of static friction is: p i , q s k ≤ f p i , n s k (8) Where f is friction factor.

The impulse should be calculated as follows p i s k = p ~ i , n s k n i + f p ~ i , n s k q i (9) where p ~ i , n s k = Δ u ∼ i , n s k n i T K i n i + f q i

The Monte Carlo method is based on the probability model. According to the model’s process, The results of simulation experiments are regarded as the approximate solution of the problem, which relax the requirement of the idealization of the theoretical model. Therefore, it can generate a more realistic simulation model, and this method is often used as a relatively exact solution to verify or validate approximate analytical solutions.

Therefore, in this paper, Monte Carlo random sampling technique combined with ETM (MCETM) can be used to simulate rockfall’s three-dimensional random migration process with the arbitrary shape when the parameters are random. Firstly, the random numbers of collision recovery coefficient and friction coefficient, which obey normal distribution, are generated. Then the standard normal distribution function and quantile are calculated, which are computed by continued fraction advancing an iterative method based on second-order expansion.

The random number of the normal distribution X ∼ N μ , σ 2 is x i , and the random number of the standard normal distribution is represented by r i , it can be expressed as x i − μ / σ ∼ N 0,1 , that is: x i = μ + σ r i (10)

For a given number of Monte Carlo simulations N, the sample values x i i = 1,2 , ⋯ , N of each simulation sample are independently and identically distributed, and the expected value is x i = μ , and the variance is D x i = σ 2 . The sample mean value μ ^ = ∑ i = 1 N x i / N is the estimated value of the unknown theoretical mean value μ . When N is sufficiently large N → ∞ , the following conditions are guaranteed with a given probability of 1 − α . 1 N ∑ i = 1 N x i − μ ≤ μ 1 − α / 2 σ N (11) where α is a given significance level.

The absolute error of estimator μ ^ = 1 / N ∑ i = 1 N x i is δ = μ ^ − μ , and given probability 1 − α satisfies the following conditions. δ ≤ μ 1 − α / 2 σ N (12)

The relative error of estimator μ ^ = 1 / N ∑ i = 1 N x i is ε = μ ^ − μ / μ , and given probability 1 − α satisfies the following conditions. ε ≤ μ 1 − α / 2 σ N (13) where μ and σ are the mean value and standard deviation of the independently distributed random variable x i , respectively. The specific simulation test flow is shown in Figure 4.

The migration process and the shape of coal/rock blocks that have stopped migration show that the shape of rockfall is generally ellipsoidal. Therefore, the transformation of the ellipsoid equation to the rockfall model to generate irregular random shapes is more in line with the actual working face. The specific modelling process is as follows: the coordinates of any point are represented by r, θ and ϕ . Where r is the distance from that point to the center of the ball, θ and ϕ are respectively latitude and longitude of the point.

The first is to determine the θ and ϕ of each vertex of the rockfall. Then, two random variables x i and x j are determined as intermediate variables by the following equation. x i = γ δ i − 0.5 + i x j = γ δ j − 0.5 + j (14) where x i ∈ i + γ δ i − 0.5 , i − γ δ i − 0.5 , x j ∈ j + γ δ j − 0.5 , j − γ δ j − 0.5 ; δ i and δ j are random variables on 0,1 ; γ ∈ 0,1 , The change amplitude of x i and x j are controlled by i ∈ 1 , n , j ∈ 1 , m ; the number of vertices of the rockfall is m ⋅ n + 2 , and n ≥ 1 , m ≥ 3 .

The θ and ϕ can be expressed as: θ = π n x i + 1 φ = 2 π m x j (15)

Plus θ = 0 and θ = π , any vertex of the rockfall can be denoted by: x = r 1 − η δ r cos ⁡ θ y = r 1 − η δ r sin ⁡ θ ⁡ cos ⁡ ϕ a z = r 1 − η δ r sin ⁡ θ ⁡ sin ⁡ ϕ b (16) where a ∈ 0,1 ; b ∈ 0 , a ; δ r is a random variable on 0,1 ; η ∈ 0,1 is the value controlling the change amplitude of r. This method can simulate the random shape of polyhedral rockfall with free sphericity and roundness. The rockfall model established by this method is shown in Figure 5.

This paper takes a steeply dipping coal seam as an example in Gansu Province. Contour lines and drilling information of the working face floor can be shown in Figure 6. The buried depth of the working face is 270 m, the north side of the working face is 201 working face for mining, the west side is the mine field boundary (30 m), the east side is the main transportation lane, and the south side is the north return air lane. The average coal seam dip angle is 35°, and the dip angle is relatively stable. The coal seam thickness is between 1.8 and 2.8 m, with an average of 2.3 m.

The contour lines and drilling information are extracted, and the drilling information is regarded as elevation points, and a working face floor mesh model is established by using cass software. As shown in Figure 7.

The model is a 3D mesh model of 4,500 m × 2,500 m × 500 m. It covers the geographic information of the entire working face floor. And the inclination of the coal seam is about 25°–35°, The upper limit of the inclination belongs to the SDCS. Rockyfor3D software is used to analyze the inclination of working face, and the results are shown in Figure 8.

It can be seen from Figure 8 that the purple area appears in the lower left corner of the figure, which is the area most prone to rockfall hazards. The inclination of this area is about 35°, which belongs to the SDCS. It is consistent with the result of contour measurement and calculation in Figure 6. Capture the region to establish a 3D model of the working face floor, as shown in Figure 9.

In this paper, hexahedral rockfall is used for modeling, and the model of working face floor is shown in Figure 9. The collision recovery coefficient of the working face floor is 0.84, the friction coefficient is 0.4, the density of the rockfall is 2,500 kg/m³, and the model size is 0.1 m × 0.2 m × 0.2 m. Rockfall is formed when coal or rock blocks leak from the roof of coal seam. The monitoring point is above the rockfall’s first falling point, the elastic modulus is 2.6 GPa, the Poisson’s ratio is 0.21, and the rockfall model is shown in Figure 10.

The rockfall moves freely from the face end of the working face at a height of 20 m, and collides with the working face floor, and then bounces and slides along the working face floor. The trajectories simulated by the ETM method and the Rockyfor3D software are compared. Figures 11, 12 show the inclination trajectories and lateral offset trajectories of rockfall along the working face floor, respectively. In Figure 12, X represents the motion track, Y represents the lateral offset track, and Z represents the height. It can be seen that the inclination trajectory and lateral migration trajectory of rockfall simulated by the ETM method are consistent with the trajectory of the Rockyfor3D simulation.

In this paper, the model of the working face floor is shown in Figure 9, and the rockfall model is modelled by the modelling method of rockfall with irregular shapes. The lithology of rockfall is sandstone, and its shape is shown in Figure 5. The collision recovery coefficient and friction coefficient were separately regarded as random variables to analyze the influence of parameters on the randomness of rockfall movement. Let the random parameters conform to the normal distribution. Their mean values are μ e = 0.875 and μ f = 0.534 (The values are based on the average values of the measured rockfall parameters in the mining area with SDCS in Gansu Province and the accuracy required by the simulation test) respectively. The importance of other parameters is the same as those of the corresponding parameters in the verification experiment. The collision recovery coefficient and friction coefficient are equivalent random variation parameters. In the experiment, when one parameter changes, the other parameters are fixed, so the calculation formula of the variation curve is the same. MCETM is used to study the random influence of each parameter. Other parameters are regarded as deterministic parameters. MCETM can usually meet the requirements for large samples ( N ≥ 30 ). To ensure the simulation accuracy, the simulation times N = 100. Taking the characteristic quantity of movement of rockfall after colliding with the working face floor for the first time as an example, this paper analyzes it.

MCETM is used for numerical simulation, and the variation curves of rockfall speed, angular velocity and energy interval width with the coefficient of variation are shown in Figure 13. In the pictures, γ represents the coefficient of variation, e represents the coefficient of collision recovery, and f represents the coefficient of friction. It can be seen from Figure 13 that with the increase of the variable coefficient of parameters, the interval width of kinematic characteristic quantities such as velocity and angular velocity shows an upward trend. This is because the impact of the randomness of the collision recovery coefficient on the interval width of kinematic characteristics such as velocity and angular velocity is far more significant than that of the randomness of the friction coefficient, and the randomness of friction coefficient has relatively little impact on the transport speed of rockfall.

Figure 14 shows the variation curves of rockfall velocity, angular velocity and energy variance with the coefficient of variation. It can be seen from Figure 14 that the dispersion of the characteristic quantity of rockfall migration is greatly influenced by the randomness of the collision recovery coefficient. In the picture, σ ² represents the variance curve. The randomness of the friction coefficient has a significant influence on the angular velocity of rockfall. The reason is that when rockfall slides along the floor of the working face, the work done by friction is converted into the rotational kinetic energy of rockfall, which increases its angular velocity. The randomness of the friction coefficient affects the angular velocity by affecting the friction force. The influence of the randomness of friction coefficient on the dispersion of kinematic characteristics of rockfall is still less than that of the randomness of collision recovery coefficient. The randomness of the collision recovery coefficient always occupies a dominant position in the process of rockfall migration.

The influence of randomness of parameters on the migration trajectory of rockfall is reflected in the trend direction of the working face and the lateral migration along the strike direction of the working face. The lateral migration of rockfall can be described according to the migration index proposed by Azzoni et al., as shown in Figure 15.

Combined with the working face of SDCS, the offset ratio is defined as: η = D 2 L (17) where D is the offset of the falling point of rockfall; L is the working face length.

The offset ratio of rockfall is shown in Figure 16. In Figures 16–22, γ(e) represents the variable coefficient of collision recovery coefficient, and γ(f) represents the variable coefficient of friction coefficient. It can be seen from Figure 16 that the migration offset ratio of rockfall increases with the increase of the variation coefficient. Furthermore, the impact curve of the randomness of the collision recovery coefficient on migration ratio is smoother, the overall curve is higher than that of the friction coefficient, and the randomness of the collision recovery coefficient has a more significant impact on migration and migration ratio of rockfall. Therefore, when the coefficient of variation is small, the randomness of the collision recovery coefficient and friction coefficient has less influence on the offset ratio, which indicates that the smaller the randomness of the parameters, the less the offset ratio is affected by it.

To compare the influence of parameters randomness on the migration of rockfall, the MCETM is used to simulate the variation curves of the upper and lower boundaries of the position coordinates of rockfall with time when the coefficient of variation of collision recovery and friction coefficient take different values, as shown in Figures 17–20. It can be seen from Figures 17–20 that the randomness of the collision recovery coefficient and friction coefficient has different effects on the tendency migration and lateral migration of rockfall, and with the increase of parameter variation coefficient, the randomness of parameters has a more significant influence on the migration of rockfall.

Figures 21, 22 are comparative diagrams of changes in inclined position coordinates and strike position coordinates of rockfall when the coefficient of variation of collision recovery and the friction coefficient is the same. It can be seen from Figures 21, 22 that the randomness of the collision recovery coefficient has a more significant influence on the migration of rockfall than the friction coefficient, and the randomness of the collision recovery coefficient still dominates the migration of rockfall.