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

Abstract

The safety production of gold, silver, copper, and other important metals is seriously threatened in the process of mining from open-pit to underground due to various factors such as infiltration caused by rainfall and unloading during mining. Furthermore, the current situation of open-pit mining in an increasing number of mines presents a high and steep terrain, which poses significant security risks. Accordingly, it is of great practical significance to investigate the failure mechanism of high-slope angles to ensure the long-term safe mining of mines, considering factors such as rainfall infiltration and excavation unloading. In this study, the slope failure of high-slope angles (45°, 55°, and 65°) under rainfall-mining coupling was analyzed using the discrete element MatDEM numerical simulation software. Herein, the stress distribution, failure characteristics, and energy conversion of the model were simulated under different slope angles to analyze the failure mechanism at each stage. The simulation results show that the damage scale is smallest at 55° and largest at 65°. This indicates that setting the slope angle to 55° can reduce the risk of slope instability. Moreover, the reduction of elastic potential energy during the mine room mining stage is similar to that of mechanical energy. During the pillar mining stage, stress is concentrated in each goaf, resulting in a greater reduction in mechanical energy compared to elastic potential energy. Finally, after the completion of the continuous pillar mining stage, stress becomes concentrated in the failure area, and the effect of the slope angle on mechanical energy reduction becomes evident after the complete collapse of the model.

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 resources 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 the 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 h 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 the 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, 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 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’s 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 the 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): 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 v2.02 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 a graphics processing unit (GPU) to perform high-performance parallel computations 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:

(1)$F_{smax}=F_{s0}-\mu F_n$

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

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

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

where K_n is the normal elastic modulus, MPa; K_s is the tangential stiffness, N/m³; and k_n1, k_n2, k_s1 and k_s2 are the properties of the four 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:

(4)$F_v=-\eta x$

where F_v is the global damping force, N; η is the damping coefficient, N/(m/s); and x is the particle velocity, m/s.

Each potential energy is defined as:

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

(6)$E_g=mgh$

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

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

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

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:

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

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

where E is Young’s modulus, MPa; d is the particle radius, m; and v is Poisson’s ratio, MPa.

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

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

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

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

(12)$\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\}}}$

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, number of particles is 76,425) is expressed as follows:

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

(14)$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}}}}$

where K is the hydraulic conductivity coefficient corresponding to saturation, m²/d; K_s is the saturated hydraulic conductivity, cm/s; S is the saturation; θ_s and θ_r are the saturated water content and initial water content, respectively; H_p is the pressure head, m; 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:

(15)$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]$

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

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

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

The microscopic transfer coefficient of particles to surrounding particles is:

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

The particle water content is:

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

where L_ij is the particle spacing based on the particle centers, m. The seepage calculation time step is updated as:

(19)$∆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}}}}$

The advantage of VGM is that it requires calibration of only the hydraulic parameters, 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 macroscopic mechanical properties of the material (Young’s modulus, Poisson’s ratio, compressive 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:

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

where Q_i is the discharge rate of particle i at the boundary, m³/s; S is the surface area of the boundary region, m²; B is the number of particles in the boundary region; and C_i is the cross-sectional area of particle i, m². 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 [10].

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 h. 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 h), 40 mm/h (0.5 h), and 40 mm/h (0.5 h), respectively. Herein, the rainfall infiltration process is simulated by setting the parameters of rainfall intensity and duration in the VGM model. Moreover, 30 monitoring points are established above the mining area in the 3D model to record the 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, and 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 distribution 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 show 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, 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 continuous 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. This 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, and 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 phenomenon, 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, 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 it 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°. This shows that different slope angle models have little difference in the mining stage. Nevertheless, 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.

Overall, according to the simulation results from the MatDEM numerical software, the stress distribution, displacement changes, and energy conversion under room-and-pillar mining with different slope angle models are summarized. Herein, the 55° slope angle model is the smallest in terms of stress distribution, displacement change, and mechanical energy reduction. Therefore, we should design the slope cutting angle as 55° where possible to reduce the instability of the slope in practical engineering.

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 concentrated 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. There is 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. 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. 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 mechanical 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.

In the future, based on this research, we will explore the characteristics of fracture evolution in slope rock masses under different slope height conditions. Additionally, we will continue to study the deformation rules of rock masses under the coupling effects of three scenarios: rainfall, excavation, and blasting.

Author Contributions

Conceptualization, Q.L., Y.W. and X.L.; methodology, Q.L.; software, Q.L.; validation, Y.W., X.L. and D.H.; formal analysis, Q.L. and S.J.; investigation, X.L.; resources, X.L.; data curation, B.G., X.L. and S.J.; writing—original draft preparation, Q.L.; writing—review and editing, X.L.; visualization, Y.W. and D.H.; supervision, B.G.; project administration, X.L.; funding acquisition, X.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Most of the data generated during this study are included in the article. For other datasets, please contact the corresponding author with reasonable requests.

Acknowledgments

We would like to express our gratitude to Tao Zhou and Jiawen Wang for their invaluable assistance in constructing the numerical model.

Conflicts of Interest

Author Yunmin Wang was employed by the company Sinosteel Maanshan General Institute of Mining Research Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Site engineering geological map.

Figure 2. Precipitation distribution in China in 2023.

Figure 3. Principle of seepage simulation.

Figure 4. Boundary conditions and direction of the seepage flow model.

Figure 5. Numerical simulation scheme and model establishment.

Figure 6. Cloud map of rainfall and excavation-induced coupling stress variation under different numerical models of high-slope angles.

Figure 7. Stress values of different monitoring points under a 45° slope angle model.

Figure 8. Stress values of different monitoring points under a 55° slope angle model.

Figure 9. Stress values of different monitoring points under a 65° slope angle model.

Figure 10. Cloud map of rainfall and excavation-induced coupling vertical displacement variation under different numerical models of high-slope angles.

Figure 11. Displacement values of different monitoring points under different slope angle models.

Figure 12. Energy change curve.

Figure 13. Energy change curves of different mining stages under a 45° slope angle model.

Figure 14. Energy change curves of different mining stages under a 55° slope angle model.

Figure 15. Energy change curves of different mining stages under a 65° slope angle model.

Table 1

Macroscopic mechanical parameters of the model material.

Modulus of Elasticity (E)/GPa	Poisson’s Ratio (λ)	Tensile Strength (Tu)/MPa	Compressive Strength (Cu)/MPa	Coefficient of Internal Friction (μi)
18.81	0.1507	24.36	187.5	0.8

Table 2

Microscopic mechanical parameters of the model material.

Normal Stiffness (Kn)/(KN/m)	Tangential Stiffness (Ks)/(KN/m)	Fracture Displacement (Xb)/m	Initial Shear Resistance (Fs0)/KN	Friction Coefficient (μ)
4.01 × 105	2.57 × 105	1.5 × 10−3	7.55 × 103	0.31

Table 3

Energy loss condition.

Slope Angle	Mining Stage	Elastic Potential Energy Reduction/J	Mechanical Energy Reduction/J
45°	Open-pit mining stage	1.5416	1.611
	Mine room 1#, 2# stage	1.8289	1.9852
	Mine room 3#, 4# stage	1.7495	1.9996
	Mine room 5#, 6# stage	1.6907	1.9871
	Pillar 1#	0.8179	1.1537
	Pillar 2#	1.1145	16.3499
	Pillar 3#	2.0565	6.2662
	Continuous pillar mining stage	11.6399	132.01
55°	Open-pit mining stage	1.7807	1.891
	Mine room mining stage	4.6668	5.3161
	Pillar mining stage	2.2785	16.4677
	Continuous pillar mining stage	11.1219	105.6281
65°	Open-pit mining stage	1.9139	2.2124
	Mine room mining stage	5.5492	6.5949
	Pillar mining stage	4.9981	27.4484
	Continuous pillar mining stage	12.0075	151.199

References

1. Hermesmann, M.; Müller, T.E. Green, Turquoise, Blue, or Grey? Environmentally friendly Hydrogen Production in Transforming Energy Systems. Prog. Energy Combust. Sci.; 2022; 90, 100996. [DOI: https://dx.doi.org/10.1016/j.pecs.2022.100996]

2. Ghisellini, P.; Ji, X.; Liu, G.Y.; Ulgiati, S. Evaluating the transition towards cleaner production in the construction and demolition sector of China: A review. J. Clean. Prod.; 2018; 195, pp. 418-434. [DOI: https://dx.doi.org/10.1016/j.jclepro.2018.05.084]

3. Luhar, S.; Luhar, I. Potential application of E-wastes in construction industry. Constr. Build. Mater.; 2019; 203, pp. 222-240. [DOI: https://dx.doi.org/10.1016/j.conbuildmat.2019.01.080]

4. Duc, N.H.; Kumar, P.; Long, P.T.; Meraj, G.; Lan, P.P.; Almazroui, M.; Avtar, R. A Systematic Review of Water Governance in Asian Countries: Challenges, Frameworks, and Pathways Toward Sustainable Development Goals. Earth Syst. Environ.; 2024; 8, pp. 181-205. [DOI: https://dx.doi.org/10.1007/s41748-024-00385-1]

5. Lucas, S.A.; Lee, C.C.C.; Love, E. Characterising Recycled Organic and Mineral Materials for Use as Filter Media in Biofiltration Systems. Water; 2019; 11, 1074. [DOI: https://dx.doi.org/10.3390/w11051074]

6. He, M.C.; Wang, Q.; Wu, Q.Y. Innovation and future of mining rock mechanics. J. Rock Mech. Geotech. Eng.; 2021; 13, pp. 1-21. [DOI: https://dx.doi.org/10.1016/j.jrmge.2020.11.005]

7. Xi, X.; Zhou, J.S.; Gao, X.Y.; Wang, Z.; Si, J.J. Impact of the global mineral trade structure on national economies based on complex network and panel quantile regression analyses. Resour. Conserv. Recycl.; 2020; 154, 104637. [DOI: https://dx.doi.org/10.1016/j.resconrec.2019.104637]

8. Li, Q.H.; Song, D.Q.; Yuan, C.M.; Nie, W. An image recognition method for the deformation area of open-pit rock slopes under variable rainfall. Measurement; 2022; 188, 110544. [DOI: https://dx.doi.org/10.1016/j.measurement.2021.110544]

9. Leong, E.C.; Rahardjo, H. Two and three-dimensional slope stability reanalyses of Bukit Batok slope. Comput. Geotech.; 2012; 42, pp. 81-88. [DOI: https://dx.doi.org/10.1016/j.compgeo.2012.01.001]

10. Li, X.S.; Li, Q.H.; Wang, Y.M.; Liu, W.; Hou, D.; Zheng, W.B.; Zhang, X. Experimental study on instability mechanism and critical intensity of rainfall of high-steep rock slopes under unsaturated conditions. Int. J. Min. Sci. Technol.; 2023; 33, pp. 1243-1260. [DOI: https://dx.doi.org/10.1016/j.ijmst.2023.07.009]

11. Li, X.S.; Li, Q.H.; Wang, Y.M.; Liu, W.; Hou, D.; Zhu, C. Effect of slope angle on fractured rock masses under combined influence of variable rainfall infiltration and excavation unloading. J. Rock Mech. Geotech. Eng.; 2024; 16, pp. 4154-4176. [DOI: https://dx.doi.org/10.1016/j.jrmge.2024.08.019]

12. Ebrahim, K.M.P.; Gomaa, S.M.M.H.; Zayed, T.; Alfalah, G. Rainfall-induced landslide prediction models, part ii: Deterministic physical and phenomenologically models. Bull. Eng. Geol. Environ.; 2024; 83, 85. [DOI: https://dx.doi.org/10.1007/s10064-024-03563-7]

13. Peruccacci, S.; Brunetti, M.T.; Gariano, S.L.; Melillo, M.; Rossi, M.; Guzzetti, F. Rainfall thresholds for possible landslide occurrence in Italy. Geomorphology; 2017; 290, pp. 39-57. [DOI: https://dx.doi.org/10.1016/j.geomorph.2017.03.031]

14. dos Santos, T.B.; Lana, M.S.; Pereira, T.M.; Canbulat, I. Quantitative hazard assessment system (Has-Q) for open pit mine slopes. Int. J. Min. Sci. Technol.; 2019; 29, pp. 419-427. [DOI: https://dx.doi.org/10.1016/j.ijmst.2018.11.005]

15. Galeandro, A.; Simunek, J.; Simeone, V. Analysis of rainfall infiltration effects on the stability of pyroclastic soil veneer affected by vertical drying shrinkage fractures. Bull. Eng. Geol. Environ.; 2013; 72, pp. 447-455. [DOI: https://dx.doi.org/10.1007/s10064-013-0492-5]

16. Dorn, R.I.I.; Walker, I.J.J. Dirt cracking as rock fracture-wedging process in the Mediterranean climate of Victoria, British Columbia, Canada. Catena; 2022; 210, 105920. [DOI: https://dx.doi.org/10.1016/j.catena.2021.105920]

17. Li, Q.H.; Wang, Y.M.; Li, X.S.; Tang, S.B.; Gong, B.; Jiang, S. Evolutionary characteristics of the fracture network in rock slopes under the combined influence of rainfall and excavation. Bull. Eng. Geol. Environ.; 2025; 84, 47. [DOI: https://dx.doi.org/10.1007/s10064-025-04084-7]

18. Li, Q.H.; Liu, W.; Jiang, L.L.; Qin, J.X.; Wang, Y.F.; Wan, J.F.; Zhu, X.S. Comprehensive safety assessment of two-well- horizontal caverns with sediment space for compressed air energy storage in low-grade salt rocks. J. Energy Storage; 2024; 102, 114037. [DOI: https://dx.doi.org/10.1016/j.est.2024.114037]

19. Fang, K.; Tang, H.M.; Li, C.D.; Su, X.X.; An, P.J.; Sun, S.X. Centrifuge modelling of landslides and landslide hazard mitigation: A review. Geosci. Front.; 2023; 14, 101493. [DOI: https://dx.doi.org/10.1016/j.gsf.2022.101493]

20. Dintwe, T.K.M.; Sasaoka, T.; Shimada, H.; Hamanaka, A.; Moses, D.N.; Peng, M.; Meng, F.F.; Liu, S.F.; Ssebadduka, R.; Onyango, J.A. Numerical Simulation of Crown Pillar Behaviour in Transition from Open Pit to Underground Mining. Geotech. Geol. Eng.; 2022; 40, pp. 2213-2229. [DOI: https://dx.doi.org/10.1007/s10706-021-02022-4]

21. Khaboushan, A.S.; Osanloo, M. A Set of Classified Integer Programming (IP) Models for Optimum Transition from Open Pit to Underground Mining Methods. Nat. Resour. Res.; 2020; 29, pp. 1543-1559. [DOI: https://dx.doi.org/10.1007/s11053-019-09551-z]

22. Jiang, N.; Zhou, C.B.; Lu, S.W.; Zhang, Z. Propagation and prediction of blasting vibration on slope in an open pit during underground mining. Tunn. Undergr. Space Technol.; 2017; 70, pp. 409-421. [DOI: https://dx.doi.org/10.1016/j.tust.2017.09.005]

23. Vyazmensky, A.; Stead, D.; Elmo, D.; Moss, A. Numerical Analysis of Block Caving-Induced Instability in Large Open Pit Slopes: A Finite Element/Discrete Element Approach. Rock Mech. Rock Eng.; 2010; 43, pp. 21-39. [DOI: https://dx.doi.org/10.1007/s00603-009-0035-3]

24. Amoushahi, S.; Grenon, M.; Locat, J.; Turmel, D. Deterministic and probabilistic stability analysis of a mining rock slope in the vicinity of a major public road—Case study of the LAB Chrysotile mine in Canada. Can. Geotech. J.; 2018; 55, pp. 1391-1404. [DOI: https://dx.doi.org/10.1139/cgj-2017-0298]

25. Shah, B.L.; Alam, A.; Bhat, M.S.; Ahsan, S.; Ali, N.; Sheikh, H.A. Extreme precipitation events and landslide activity in the Kashmir Himalaya. Bull. Eng. Geol. Environ.; 2023; 82, 328. [DOI: https://dx.doi.org/10.1007/s10064-023-03350-w]

26. Ali, A.; Huang, J.S.; Lyamin, A.V.; Sloan, S.W.; Griffiths, D.V.; Cassidy, M.J.; Li, J.H. Simplified quantitative risk assessment of rainfall-induced landslides modelled by infinite slopes. Eng. Geol.; 2014; 179, pp. 102-116. [DOI: https://dx.doi.org/10.1016/j.enggeo.2014.06.024]

27. Mondini, A.C.; Guzzetti, F.; Melillo, M. Deep learning forecast of rainfall-induced shallow landslides. Nat. Commun.; 2023; 14, 2466. [DOI: https://dx.doi.org/10.1038/s41467-023-38135-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37117189]

28. Sun, S.W.; Liu, L.; Hu, J.B.; Ding, H. Failure characteristics and mechanism of a rain-triggered landslide in the northern longwall of Fushun west open pit, China. Landslides; 2022; 19, pp. 2439-2458. [DOI: https://dx.doi.org/10.1007/s10346-022-01926-3]

29. Deliveris, A.V.; Theocharis, A.I.; Koukouzas, N.C.; Zevgolis, I.E. Numerical Slope Stability Analysis of Deep Excavations Under Rainfall Infiltration. Geotech. Geol. Eng.; 2022; 40, pp. 4023-4039. [DOI: https://dx.doi.org/10.1007/s10706-022-02135-4]

30. Scaringi, G.; Fan, X.M.; Xu, Q.; Liu, C.; Ouyang, C.J.; Domènech, G.; Yang, F.; Dai, L.X. Some considerations on the use of numerical methods to simulate past landslides and possible new failures: The case of the recent Xinmo landslide (Sichuan, China). Landslides; 2018; 15, pp. 1359-1375. [DOI: https://dx.doi.org/10.1007/s10346-018-0953-9]

31. Liu, W.; Li, Q.H.; Yang, C.H.; Shi, X.L.; Wan, J.F.; Jurado, M.J.; Li, Y.P.; Jiang, D.Y.; Chen, J.; Qiao, W.B. et al. The role of underground salt caverns for large-scale energy storage: A review and prospects. Energy Storage Mater.; 2023; 63, 103045. [DOI: https://dx.doi.org/10.1016/j.ensm.2023.103045]

32. Song, D.Q.; Du, H. Numerical Investigation of the Evolution Process of an Open-Pit Mine Landslide Using Discrete-Element Method. Int. J. Geomech.; 2023; 23, 04023054. [DOI: https://dx.doi.org/10.1061/IJGNAI.GMENG-7568]

33. Singh, M.; Singh, B. Modified Mohr-Coulomb criterion for non-linear triaxial and polyaxial strength of jointed rocks. Int. J. Rock Mech. Min. Sci.; 2012; 51, pp. 43-52. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2011.12.007]

34. Marino, P.; Santonastaso, G.F.; Fan, X.M.; Greco, R. Prediction of shallow landslides in pyroclastic-covered slopes by coupled modeling of unsaturated and saturated groundwater flow. Landslides; 2021; 18, pp. 31-41. [DOI: https://dx.doi.org/10.1007/s10346-020-01484-6]

35. Fatichi, S.; Or, D.; Walko, R.; Vereecken, H.; Young, M.H.; Ghezzehei, T.A.; Hengl, T.; Kollet, S.; Agam, N.; Avissar, R. Soil structure is an important omission in Earth System Models. Nat. Commun.; 2020; 11, 522. [DOI: https://dx.doi.org/10.1038/s41467-020-14411-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31988306]

36. Zhu, Y.; Liu, C.; Zhang, H.Y.; Zhao, C.; Wang, B.J.; Mao, M.Y.; Geng, H. Micro mechanism investigation of hydraulic fracturing process based a fluid-solid coupling discrete element model. Comput. Geotech.; 2024; 174, 106640. [DOI: https://dx.doi.org/10.1016/j.compgeo.2024.106640]

37. Hayek, M. Analytical solution for steady vertical flux through unsaturated soils based on van Genuchten-Mualem model. J. Hydrol.; 2024; 634, 131066. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2024.131066]

38. Younes, A.; Mara, T.; Fahs, M.; Grunberger, O.; Ackerer, P. Hydraulic and transport parameter assessment using column infiltration experiments. Hydrol. Earth Syst. Sc.; 2017; 21, pp. 2263-2275. [DOI: https://dx.doi.org/10.5194/hess-21-2263-2017]

39. Khakpour, M.; Mirghasemi, A.A. Macro-micro mechanical study of principal stress rotation in granular materials using DEM simulations of hollow cylinder test. Powder Technol.; 2023; 425, 118580. [DOI: https://dx.doi.org/10.1016/j.powtec.2023.118580]

40. Paluszny, A.; Zimmerman, R.W. Modelling of primary fragmentation in block caving mines using a finite-element based fracture mechanics approach. Geomech. Geophys. Geo-Energy Geo-Resour.; 2017; 3, pp. 121-130. [DOI: https://dx.doi.org/10.1007/s40948-016-0048-9]

41. Yan, B.X.; Jia, H.W.; Yang, Z.; Yilmaz, E.; Liu, H.L. Goaf instability in an open pit iron mine triggered by dynamics disturbance: A large-scale similar simulation. Int. J. Min. Reclam. Environ.; 2023; 37, pp. 606-629. [DOI: https://dx.doi.org/10.1080/17480930.2023.2233866]

42. Liu, C.; Xu, Q.; Shi, B.; Deng, S.; Zhu, H.H. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks. Comput. Geosci.; 2017; 103, pp. 12-20. [DOI: https://dx.doi.org/10.1016/j.cageo.2017.03.003]

43. Wang, L.; Zhu, Y.J.; Mei, C. Numerical modeling of cracking behaviors for the rock-inclusion composite under dynamic tensile loading. Comput. Geotech.; 2023; 157, 105325. [DOI: https://dx.doi.org/10.1016/j.compgeo.2023.105325]

44. Rálek, P.; Sosna, K.; Hokr, M.; Franek, J. In situ measurements of thermally loaded rock and evaluation of an experiment with a 3D numerical model. Int. J. Rock Mech. Min. Sci.; 2018; 105, pp. 1-10. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2018.01.042]

45. Kurosawa, K.; Genda, H. Effects of Friction and Plastic Deformation in Shock-Comminuted Damaged Rocks on Impact Heating. Geophys. Res. Lett.; 2018; 45, pp. 620-626. [DOI: https://dx.doi.org/10.1002/2017GL076285]

46. Sun, H.D.; He, H.X.; Cheng, Y.; Gao, X.J. Numerical simulation and dynamic experimental investigation on vibration control of a structural system with a stacked single particle-inerter damper. Earthq. Eng. Struct. Dyn.; 2024; 53, pp. 3649-3677. [DOI: https://dx.doi.org/10.1002/eqe.4191]