Simulating the Failure Mechanism of High Slope Angles Under Rainfall-Mining Coupling Using MatDEM

1. Introduction

Mineral resources serve as the foundation for national industrial development [1,2]. To meet the demands of modern industry and infrastructure construction, an increasing number of mineral resources are being utilized [3]. According to a report released by the Organization for Economic Cooperation and Development in 2023, global raw material consumption has increased by 2.25 times from 1970 to 2023 [4]. It is expected to further grow to 167 billion tons by 2060, with the use of metal minerals increasing from 9 billion tons to 20 billion tons [5]. Although China is a big country in terms of mineral resources, its per capita share is very small, and its per capita reserves of major metals are less than one-fourth of the world's per capita reserves [6]. Furthermore, a significant amount of mineral resources need to be imported, and the external dependence on certain resour- ces exceeds 50% and continues to increase year by year [7]. At present, 80% of China's ore resources produced from metallurgical mines are obtained through open-pit mining [8]. With the high-intensity excavation of metal mines, most large- and medium-sized open-pit mines have transitioned to deep concave mining, and some mines have further transitioned from deep concave mining to underground mining. Hence, many large rock slopes with high and steep vertical depths of over 400 m and slope angles exceeding 40° are formed [9,10,11]. Due to the impact of open-pit mining and environmental changes, slope safety has gradually emerged as a significant challenge.

In the natural environment, the most common external influences on mining areas are rock mass erosion caused by rainfall infiltration and rock mass destruction resulting from mining activities [12,13,14]. Hence, the analysis of the influence of rainfall infiltration and mining unloading on the stability of rock masses has always been a research focus. The properties of the rock mass itself will be weakened by rainfall, and the connection between the rock particles will also be weakened. Meanwhile, it will also create a water wedge effect and dissolution effect, thereby compromising the stability of the rock mass. Moreover, the infiltration of free water into the rock mass generates pore water pressure and triggers expansion and strengthening phenomena, resulting in the development of cracks within the rock mass [15,16].

In recent years, extensive research has been conducted by many scholars on the stability of rock slopes under two conditions: unloading from open-pit excavation and infiltration caused by rainfall [17,18,19]. Dintwe et al. monitored the stress distribution and failure mechanism around the roof pillar during underground mining using FLAC3D numerical simulation. It was discovered that the stress generated near the slope of the roof pillar becomes a load, subsequently affecting the deformation process of the roof pillar [20]. Khaboushan and Osanloo developed a set of classified integer programming models to determine the optimal transition depth between open-pit mines and various underground mining methods [21]. Furthermore, Jiang et al. established a 3D numerical model of underground mining in an open-pit mine and verified its reliability using field monitoring data [22]. Vyazmensky et al. used the finite/discrete element modeling (FEM /DEM) method to analyze the development of step failure induced by slope block caving in a large open-pit mine, and the results demonstrate that the conversion from open-pit to underground mining in Parabora mine is caused by caving [23]. Besides, Amoushahi et al. utilized the finite element shear strength reduction method and limit equilibrium method to conduct deterministic and probabilistic analyses of slope stability. The results from all analyses indicate that the current mining slopes at LAB Chrysotile meet the acceptable design criteria limits [24].

The effect of rainfall on landslides is a dynamic process. Atmospheric precipitation infiltrates into the sliding body, increasing the water content and bulk weight of rocks and soil, thereby softening them and reducing their shear strength [25]. Herein Ali et al. examined how the spatial variability of saturated hydraulic conductivity and triggering mechanisms affect the risk of rainfall-induced landslides [26]. Mondini et al. proposed a deep learning-based strategy that correlates rainfall with the occurrence of landslides, indicating the possibility of effectively predicting large-scale rainfall-induced landslides [27]. Thus, this study offers the potential for predicting slippery slope operations based on rainfall measurements and quantitative weather forecasting. Moreover, Sun et al. used a numerical simulation method that combines seepage flow and stability to explore the formation mechanism of landslides. They found that with the increase in rainfall intensity and duration, the maximum slope suction continuously decreased, and the north highwall became unstable after 7.2 hours when the rainfall reached 16 mm/h [28]. Based on numerical calculations and the Monte Carlo method, Deliveris et al. discussed the significance of slope geometry, rainfall intensity, soil properties, and other related parameters in relation to extremely deep slopes in a lignite mine. It was found that the most influential factors were rainfall intensity, soil permeability, and slope height [29]. Li et al. discussed the influence of rainfall infiltration and excavation unloading on the stability of a mine slope, primarily considering two factors: the angle of excavation slope and the intensity of rainfall. The results indicate that an increase in slope angle (45°, 50°, and 55°) accelerates the process of slope instability, while an increase in rainfall intensity expands the range of landslide deformation [10,11].

Although the aforementioned scholars have utilized similar physical experiments, numerical simulations, and machine learning methods to conduct extensive research on the stability of open-pit slope mining and the mechanism analysis of rainfall-induced landslides. Nevertheless, there have been few studies on the coupling between rainfall infiltration and excavation unloading. Based on the excavation status of most open-pit mines in China, a considerable proportion of slope angles are relatively high (generally defined as angles greater than 45°). Coupled with frequent heavy rainfall in the southern region, this poses a serious threat to mine production and regional security. In previous studies, it has been shown that an increase in slope angle accelerates slope instability [9,10,11]. Nevertheless, the impact of a high angle on mine slope stability has become a challenging problem due to the simultaneous effects of open-pit to underground mining and rainfall infiltration.

To address these challenges, this research investigates the deformation evolution characteristics and the effects of slope angles (45°, 55°, and 65°) on the rock mass in mine slopes under rainfall-mining coupling utilizing the discrete element MatDEM numerical calculation method. The innovation of this study lies in the model established by the van Genuchten model (VGM), which discusses the displacement, stress change rules, failure characteristics, and energy conversion modes of overlying strata under variable rainfall conditions during open-pit to underground mining. Accordingly, this study can provide theoretical guidance and technical support for the safe open-pit to underground mining of high-steep rock slopes worldwide that are affected by heavy rainfall.

2. Materials and Methods

2.1. Engineering Background

Herein the engineering study is based on the slope of Huangniuqian in the mining area of Dexing Copper Mine, located in Sizhou Town, Dexing City, Jiangxi Province (Figure 1). Dexing Copper Mine is one of the key copper mines in China, with large mineral reserves and deep-lying ore bodies. However, its ore exhibits high selectivity and contains numerous components suitable for comprehensive utilization. Currently, under the jurisdiction of Jiangxi Copper Co., LTD., it is approximately 20 km away from the center of Dexing City. The study area is located between 29°01'40"N and 29°02'30"N, as well as between 117°41'26"E and 117°42'23"E. It belongs to a subtropical monsoon climate characterized by high temperatures and year-round rainfall. The site geological occurrence extends from over 60 m below ground to up to 1000 m above the surface. The exploitation of mineral resources in the past 30 years has transformed Dexing Copper mine from open-pit mining to underground mining, resulting in the formation of some projects with slope angles exceeding 45°. These projects are susceptible to destabilization and damage caused by external factors such as heavy rainfall. The Huangniuqian slope has an overall height of approximately 540 m and a length of around 400 m. It strikes between 144° to 152° and fault dips between 234° to 242°. Additionally, there are several man-made soil platforms ranging from 0.5 to 2.5 m on the surface of the slope, with the main mining area located below and surrounded by transportation lines. In addition, we investigated the precipitation situation in China in 2023 (Figure 2), and the rainfall in Jiangxi Province is approximately 800-1600 mm, with Dexing City receiving close to 1600 mm of precipitation. Therefore, extreme construction weather conditions such as heavy rainfall and high slope angles make the transition from open-pit mining to underground mining more challenging.

2.2. MatDEM Numerical Simulation Method

To better analyze the stress distribution of rock mass during the transition from open-pit to underground mining under the coupling effect of rainfall and excavation, MatDEM numerical simulation software was used to simulate high slopes with angles of 45°, 55°, and 65°. In this study, the implicit computing format is adopted, which can utilize Graphics Processing Unit (GPU) to perform high-performance parallel comput- ations at a faster speed than other numerical simulation software [30]. MatDEM adopts the matrix discrete element calculation method, which can perform multi-field and fluid- structure coupling analyses, and can be applied to industrial analysis, landslide disaster simulation, rock failure, and more, with a broad application prospect [31,32].

2.2.1. Basic Principles of Discrete Element

In practical engineering, rock mass is influenced by various factors and undergoes different changes. To facilitate the numerical simulation, the establishment of the model needs to be based on some assumptions. Firstly, the model material is considered to be an ideal elastic-plastic body obeying the Mohr-Coulomb criterion. Secondly, the under- ground geological structure is ignored, so that the rock mass material is homogeneous. Finally, the rainwater infiltration capacity of the rock slope is considered to be uniform [33,34,35]. Accordingly, based on the contact model with linear elasticity in MatDEM and according to the Mohr-Coulomb criterion:

$$F_{smax}=F_{s0}-\mu F_n$$

(1)

where F_smax is the maximum shear force, F_s0 is the adhesion between particles, μ is the friction coefficient, and F_n is the normal pressure.

$$K_n={\displaystyle \frac{k_{n1}{k}_{n2}}{k_{n1}+k_{n2}}}$$

(2)

$$K_s={\displaystyle \frac{k_{s1}{k}_{s2}}{k_{s1}+k_{s2}}}$$

(3)

where K_n is the normal elastic modulus, K_s is the tangential stiffness, as well as k₁ and k_n are the properties of the two particles themselves.

For the energy property, the energy consumption of the object is represented by the conversion of mechanical energy into kinetic and heat energy in the numerical model, which is calculated using the damping force. The expression is as follows:

$$F_v=-\eta x$$

(4)

where F_v is the global damping force, η is the damping coefficient, and x is the particle velocity.

Each potential energy is defined as:

$$E_e={\displaystyle \frac{K_n{X}_n^2+K_s{X}_s^2}{2}}$$

(5)

$$E_g=mgh$$

(6)

where E_e is the elastic potential energy, E_g is the gravitational potential energy, X_n is the normal displacement, X_s is the tangential displacement, m is the particle mass, g is the gravitational acceleration, and h is the relative height of the plane.

The form of energy conversion during interparticle fracture is as follows:

$$E={\displaystyle \frac{F_{smax}^2-{F^{\prime }}_{smax}^2}{2K_s}}$$

(7)

To correlate the macroscopic parameters with the microscopic parameters of the particles, quantitative processing is required in the establishment of the model. Here, the specific scaling equation is as follows:

$$K_n={\displaystyle \frac{\sqrt{2}Ed}{4(1-2v)}}$$

(8)

$$K_s={\displaystyle \frac{\sqrt{2}(1-5v)Ed}{4(1+v)(1-2v)}}$$

(9)

where E is Young's modulus, d is the particle radius, and v is Poisson's ratio.

$$X_b={\displaystyle \frac{3K_n+K_s}{6\sqrt{2}{K}_n(K_n+K_s)}}{T}_u{d}^2$$

(10)

where X_b is the fracture displacement, and T_u is the tensile strength.

$$F_{s0}={\displaystyle \frac{1-\sqrt{2}{\mu }_i}{6}}{C}_u{d}^2$$

(11)

where μ_i is the coefficient of internal friction, and C_u is the compressive strength.

$$\mu ={\displaystyle \frac{\sqrt{2}\left\{{[{(1+{\mu }_i^2)}^{0.5}+{\mu }_i]}^2-2\right\}}{2\left\{{[{(1+{\mu }_i^2)}^{0.5}+{\mu }_i]}^2+1\right\}}}$$

(12)

2.2.2. Unsaturated Seepage Model of Rock Mass

Based on previous studies on numerical simulation of seepage flow, it is believed that each rock and soil particle can be analyzed as a collection of discrete particles and pores, as well as water flow is transmitted between these particles [36]. In particular, van Genuchten proposed an equation for the soil water content-pressure head curve [37,38]. The closed-form expression for relative hydraulic conductivity can be derived from this equation, which can be used to accurately predict unsaturated water conductivity. The van Genuchten model (VGM) is expressed as follows:

$$K=K_s{S}^{0.5}{[1-{(1-S^{{\displaystyle \frac{n-1}{n}}})}^{{\displaystyle \frac{1}{n}}}]}^2$$

(13)

$$S={\displaystyle \frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}}=={\displaystyle \frac{1}{{[1+(\alpha {\left|H_p\right|}^n\left)\right]}^{\frac{1}{n}}}}$$

(14)

where K is the hydraulic conductivity coefficient corresponding to saturation; K_s is the saturated hydraulic conductivity; S is the saturation; θ_s and θ_r are the saturated water content and initial water content, respectively; H_p is the pressure head; and α and n are fitting parameters.

Equations (13) and (14) can be used to calculate the saturation of each particle in the VGM, and Equation (15) can be used to calculate the expansion coefficient.

$$K_d={\displaystyle \frac{(1-\frac{1}{n})}{\alpha \frac{1}{n}({\theta }_s-{\theta }_r)}}{S}^{0.5-n}[{(1-S^n)}^{-{\displaystyle \frac{1}{n}}}+{(1-S^n)}^{{\displaystyle \frac{1}{n}}}-2]$$

(15)

Thus, the particle microscopic transfer coefficient (three-dimensional calculation) at each time step in the numerical simulation is:

$$d_i={\displaystyle \frac{2K_d\left({\theta }_i\right)}{{\sum }_{c=1}^n{I}_i}}$$

(16)

where θ_i is the water content at the current time, and I_i is the distance from the particle center to the contact point.

The microscopic transfer coefficient of particles to surrounding particles is:

$$d_{ij}={\displaystyle \frac{2d_i{d}_j}{d_i+d_j}}$$

(17)

The particle water content is:

$${\theta }_i^{t+∆t}={\theta }_i^t+∆t\left({\displaystyle \frac{1}{V_i}}{\sum }_{j=1}^n{d}_{ij}{\displaystyle \frac{∆{\theta }_{ij}}{L_{ij}}}\right)$$

(18)

where L_ij is the particle spacing based on the particle centers. Update seepage calculation time step.

$$∆t\le {(∆z)}^2{\displaystyle \frac{Fn}{K_d}}\mathrm{˄}∆t\le {\displaystyle \frac{V_i}{{\sum }_{j=1}^n\frac{d_{ij}}{L_{ij}}}}$$

(19)

The advantage of VGM is that it requires calibration of only the hydraulic para- meters, and the boundary conditions can be adjusted for rainfall or water infiltration. As illustrated in Figure 3, using measured data and curve fitting, the hydraulic parameters in seepage simulation were calibrated as follows: θ=0.15, θ_s=0.368, θ_r=0.102, K_s=0.0922 cm/s, fitting parameters: n=2, α=0.00335.

2.2.3. Model Parameters and Boundary Conditions

Herein the mechanical parameters of the rock masses should be determined before numerical simulation. In the MatDEM numerical software, the model is formed by the random accumulation of particles of various sizes, and there is no direct correspondence between the micro-mechanical parameters of the DEM model particles and the macro- mechanical properties of the rock [39]. Besides, the particle micromechanical parameters obtained by MatDEM through macro- and micro-mechanical transformation equations cannot meet practical requirements. Accordingly, it is necessary to calibrate the micro- mechanical parameters of the particles to match the macroscopic mechanical behavior of the model material. In this study, the material training box module of MatDEM software is utilized. The module integrates uniaxial compression, tensile, and other unit tests, and adjusts the microscopic parameters of the model by automatically testing five macro- scopic mechanical properties of the material (Young's modulus, Poisson's ratio, compre- ssive strength, tensile strength, as well as internal friction angle) through continuous modification of the proportionality coefficient, so as to achieve the convergence of each mechanical property to the set value. After four rounds of parameter training, the error in the rock masses parameters can be reduced to 2%. The macro- and micro-mechanical parameters of the materials after training are presented in Table 1 and Table 2.

In the simulation process, displacement constraints are set on the front, back, left, right, and bottom boundaries of the numerical model, and the top boundary is kept free (Figure 4). The model takes into account the effect of gravity, with a g value of -9.8 m/s². The Neumann-type boundary conditions with known water fluxes are used to simulate unsaturated flow [40]. In other words, given the rainfall intensity q, the slope surface is designated as the infiltration boundary. When the volumetric water content of adjacent particles at the upper boundary forms a hydraulic difference, the seepage channel opens, and the water migrates downward. Considering that each particle occupies a certain equivalent continuous area, the q is converted into the discharge rate for each particle on the boundary, which is expressed as follows:

$$Q_i=d_i{\displaystyle \frac{∆{\theta }_i}{L_i}}{\displaystyle \frac{qS{C}_i}{{\sum }_{i=1}^B{C}_i}}$$

(20)

where Q_i is the discharge rate of particle i at the boundary, S is the surface area of the boundary region, B is the number of particles in the boundary region, and C_i is the cross- sectional area of particle i. By substituting Q_i into Equation (18), the infiltration process of boundary particles under constant rainfall intensity can be realized. Before the seepage simulation began, each particle was assigned initial water content information, and zero flux boundary conditions were applied to both sides of the model. After the seepage calculation, the slope's strength parameters are adjusted based on the particle volume water content information.

2.2.4. The establishment of Numerical Models and Simulation Schemes

As depicted in Figure 5, the transition from open-pit to underground excavation is mainly divided into two stages: open-pit slope cutting and underground room and pillar mining. The length, width, and height of the numerical model established by MatDEM are 160 cm, 30 cm, and 120 cm, respectively. Considering the effect of high slope angles, the open-pit slope is divided into three numerical models: 45°, 55°, and 65°. At this stage, the rainfall intensity is 10 mm/h, and the rainfall duration is 1 hour. The mining stage of the underground room and pillar method revolves around the mining sequence of the first mine rooms (1# to 6#), the second pillars (1# to 3#), and the final continuous pillars (1# to 2#). Detailed mining sizes are shown in Figure 5. In this stage, the rainfall intensity (time) in the mine room, pillar, and continuous pillar is 20 mm/h (1 hour), 40 mm/h (0.5 hours), and 40 mm/h (0.5 hours), respectively. Furthermore, a 3D numerical simulation model is established based on the actual engineering conditions, and 30 measurement points are set up above the mining area in the model to record its displacement changes.

3. Results and Discussion

3.1. Numerical Simulation Results: Stress Distribution Characteristics

Figure 6 illustrates the stress distribution after the underground mining stage at different high slope angles. It can be seen that after the mining stage of the mine room, the stress is mainly concentrated in the upper failure zones of 1# and 2#, while the stress concentration is less pronounced in the upper parts of the other two groups of mine rooms. When the pillar mining stage is completed, part of the roof of the mine room collapses due to its inability to withstand the subsidence force, leading to caving, as well as the stress becomes concentrated in the goaf [41]. After the end of the continuous pillar mining stage, a large-scale collapse occurs; the fracture extends upwards, the failure area develops significantly, and stress concentrates in this area. By comparing the stress distr- ibution at different stages and under various slope angles, it can be found that the slope angle has little influence on the stress distribution during the mining stage of the mine room. Nevertheless, the stress concentration area in each failure area differs during the continuous pillar mining stage. Herein the 55° model has the smallest area, while the 65° model has the largest.

To analyze the stress variations in different locations, Figure 7, Figure 8 and Figure 9 shows the stress values recorded at each measurement point after the underground mining stage of different mine rooms and pillars. Here, the initial stress value refers to the value after the open-pit mining stage. During the mining stage, the stress changes significantly from No. 3 to No. 28 measuring points. In other words, the stress above the goaf is concentrated and deformed. With the mining of the mine room, the stress values at each measuring point are gradually reduced, the overall strength is decreased, and the supporting force is mainly provided by the pillars and the continuous pillars. Besides, the stage of stress change varies according to the different mine rooms. For example, the measuring point above the mine room of the 45° slope angle model has the largest difference between its initial value and those of 1# and 2# (Figure 7). In the 55° slope angle model, the stress value at the measuring points above the mine room decreases the most at 3# and 4#, and changes little at 5# and 6# (Figure 8). Meanwhile, the measurement points above the mine room in the 65° slope angle model decrease the most at 5# and 6#, and remain almost unchanged at 1# and 2# to 3# and 4# (Figure 9). In this research, the maximum mean stress change is 43.2 kPa above pillar 1#. In the underground mining stage, the stress value continues to decrease, and the maximum average stress change is 46.67 kPa, which is located above continuous pillar 1#. Moreover, the maximum average value of stress variation in the continuous pillar mining stage is 32.7 kPa, which is located above mine room 1#.

Compared with the stress variations at different high slope angles, the overall stress variation pattern is similar. The difference lies in the varying reduction values of stress. In the mining stage, the reduction value for the 65° model is higher than that of the 45° and 55° models. Furthermore, during the mining stage of the pillar, the change in the 55° model is minimal, with changes remaining roughly consistent throughout the continu- ous pillar mining stage. Among them, the overall maximum stress difference for the 65° model is 51.5 kPa, for the 55° model it is 43.3 kPa, and for the 45° model it is 48.61 kPa. Accordingly, this fully demonstrates that the 65° slope angle model exhibits a greater conversion of mechanical energy into other forms during the numerical simulation of transitioning from open-pit to underground mining under the coupled effect of rainfall and mining activities, resulting in the highest degree of damage.

3.2. Numerical Simulation Results: Displacement Evolution Rule

Figure 10 shows the Z-direction displacement variation cloud map of different high -slope angle models in numerical simulation. Taking the negative direction of the Z-axis as the downward displacement, it can be found that the main displacement occurs above the mine room during the mining stage, forming an arc-shaped area. It indicates that the vertical crack bends from the edge of the mine room to the arc crack, nevertheless, the displacement area is not extensive, and there is no large-scale damage. After entering the stage of pillar mining, different displacement areas tend to become interconnected. At this time, cracks start to link up, and the settlement displacement of rock strata results in separation cracks. As continuous pillar mining is completed, a large displacement area forms above the goaf, and the cracks fully propagate and connect, causing the roof to collapse.

By comparing the numerical simulation displacement diagrams of different slope angles, it can be found that there is little difference in the displacement area during the mining stage of the mine room, but there is a difference during the mining stage of the pillar. Among them, in the 45° slope angle model, the displacement areas above the first two groups of mine rooms are connected, and there is a sinking trend above continuous pillar 1#. Moreover, in the 55° slope angle model, the displacement areas are located above the mine rooms, and there is no trend of mutual connectivity. The displacement area of the 65° slope angle model is basically connected after the end of the pillar mining stage, and subsidence occurs above continuous pillars 1# and 2# (Figure 10). After the end of the continuous pillar mining stage, the displacement areas under different slope angles are clearly distinct. Herein, the displacement area of the 45° slope angle model has extended to the edge of the upper surface; the 65° model has the most significant impact on the overlying rock displacement, as well as the 55° model primarily affects the upper part of the goaf. This shows once again that different slope angle models have different influences on the subsidence of overlying rock under numerical simulation. A 55° slope angle has the least influence, followed by a 45° slope angle, and a 65° slope angle has the most influence on the stability of overlying rock.

To better analyze the displacement in the numerical simulation, Figure 11 depicts the vertical displacement recorded by 30 monitoring points at the end of each stage. It can be seen that the displacement during the mining stage is very small and relatively stable. The displacement of different slope angle models in the mining stage is all less than -0.00435 mm. When it comes to the pillar mining stage, subsidence will occur in this area every time a pillar is mined. After the pillar mining stage, the average sinking distance of pillar 1# under the three slope angles is -0.390 mm, the pillar with the largest displacement is 2#, with an average displacement of -0.422 mm, while the pillar with the smallest displacement is 3#, with an average displacement of -0.154 mm. It is worth noting that the displacement of pillar 3# at 45° and 55° slope angles is less than -0.065 mm. However, the 65° slope angle model has a displacement of -0.334 mm at point 3# after the completion of the pillar mining stage. Furthermore, during the mining stage of pillar 1#, the 65° slope angle model for pillar 2# already shows a slight subsidence pheno -menon, while the 45° and 55° slope angles do not. It is verified that the 65° slope angle model is much stronger than the 45° and 55° models in the pillar mining stage. Finally, after the continuous pillar mining stage ends, the goaf subsidence is significant, with mine rooms 1# to 4# experiencing the greatest subsidence and the most severe damage. Among them, the maximum displacement is at the 65° slope angle, measuring point No. 13 (above mine rooms 2# and 3#), which is -0.605 mm. The minimum displacement is -0.41 mm at measuring point No. 17 (in the continuous pillar 2# area), with a 55° slope angle.

3.3. Discussion on Numerical Simulation of Energy Conversion

In DEM numerical simulations, the conversion of energy can be utilized to analyze the failure of the model [42,43,44,45]. Figure 12 illustrates the overall energy variation of the numerical simulation. Taking the slope angle of 45° as an example, the gravitational potential energy of the model decreases continuously with the progress of mining, and the gravitational potential energy, elastic potential energy, and kinetic energy convert to each other. There is damping in the numerical model, leading to damping heat, friction heat generated by particle displacement and friction between particles, as well as energy released from inter-particle fractures [46], which results in part of the mechanical energy being converted into heat energy. Accordingly, the heat energy continuously increases while the mechanical energy decreases. However, the total energy in the system remains constant throughout the simulation, indicating that the MatDEM numerical simulation satisfies the rule of conservation of energy.

Figure 13, Figure 14 and Figure 15 show the changes in mechanical energy and elastic potential energy at different high slope angles during various mining stages. Table 3 records the energy loss values at different stages. In detail, both the elastic potential energy and mechanical energy decrease during the open-pit mining stage. During the mining stage, the energy loss of each mine room varies, with rooms 1#~4# experiencing a higher reduction in mining energy, while rooms 5#~6# show a lesser loss. This indicates that the energy loss in rooms 1#~4# has a greater impact on the stability of the numerical model, although the model remains stable. The decrease in elastic potential energy is analogous to the reduction in mechanical energy. When the model enters the stage of pillar mining, it begins to collapse, and the proportion of mechanical energy reduction is much greater than that of elastic potential energy reduction. Mining pillar 2# causes a significant energy reduction, indicating that pillar 2# plays a decisive role in the stability of the model. Furthermore, the continuous pillar mining increases the damage degree; cracks develop rapidly, mechanical energy decreases significantly, and the top of the goaf collapses entirely.

The influence of slope angle on energy reduction was analyzed and compared. In the open-pit mining stage, the steeper the slope angle, the greater the reduction in elastic potential energy and mechanical energy. However, the overall energy reduction value does not differ significantly. Here, the difference in energy between a 45° and a 55° slope angle is 1.16 J, and between a 45° and a 65° slope angle is 1.4277 J. In the mining stage, the energy reduction values of different slope angles vary; at 45°, it is 0.602 J more than at 55°, and 0.2801 J less than at 65°. It shows that different slope angle models have little difference in the mining stage. However, when entering the pillar mining stage and the continuous pillar mining stage, the slope angle effect on mechanical energy reduction becomes more pronounced as the model begins to collapse. The overall rule is that the mechanical energy reduction value for the 65° model is the greatest, followed by the 45° model, and the least is the 55° model.

4. Conclusions

In this research, MatDEM numerical simulation software is utilized to establish an open-pit to underground mining model based on VGM under rainfall conditions, as well as to investigate the rock mass deformation under different high slope angles (45°, 55°, and 65°) under the coupling action of rainfall and mining. Furthermore, combined with the numerical simulation results of stress distribution, overlying rock failure patterns, and energy conversion analysis, the main conclusions are as follows.

(1)

In the three groups of high slope angle models, the stress during the mine room mining stage is mainly concentrated above 1# and 2#. After the completion of the pillar mining stage, the upper part of the mine room collapses, and the stress becomes concen- trated in each goaf. Additionally, after the continuous pillar mining stage ends, a large- scale collapse occurs, and stress concentrates in the deformed area. On the other hand, the high slope angle has little influence on the stress distribution during the mine room mining stage. In the continuous pillar mining stage, the stress distribution area of the 55° model is smaller than those of the 45° and 65° models. Thus, it is suggested that the open -cut slope angle should be controlled at 55° in practical engineering, which is conducive to avoiding the risk of slope instability.

(2)

The numerical simulation results show that the vertical displacement in the mining stage appears over the mine room, forming a small arc-shaped area. When in the pillar mining stage, different displacement areas show a tendency to interconnect. Here, the 55° slope angle model is relatively stable, with no subsidence occurring in any pillar. Continuous pillar 1# subsidence under the 45° slope angle model. The displacement area of the 65° slope angle model is basically connected, and continuous pillars 1# and 2# both show signs of subsidence. When the continuous pillar mining is completed, a large area of displacement occurs above the goaf, and the roof collapses. Hence, the support work for the overlying rock should be done well during the underground mining using the room-and-pillar method, especially in the stage of continuous pillar extraction.

(3)

Comparing the energy conversion under different slope angles at each stage, the greater the slope angle, the greater the reduction of elastic potential energy and mech -anical energy during the open-pit mining stage. Subsequently, mining rooms 1# to 4# have a significant impact on the stability of the model. In the stage of pillar mining, the proportion of mechanical energy reduction is much greater than that of elastic potential energy reduction, and pillar 2# plays a decisive role in the stability of the model. Finally, in the continuous pillar mining stage, the high slope angle reduction effect of mechanical energy becomes obvious after the model begins to collapse. The overall rule is that the mechanical energy reduction value of the 65° model is the largest, followed by the 45° model, and the smallest is the 55° model. This is consistent with the stress distribution and vertical displacement evolution results.