Optimum Selection of Mining Plans for Pillars Containing Interlayers Based on Numerical Simulation

School of Civil Engineering, Shaoxing University, Shaoxing, Zhejiang 312000, China Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province, Shaoxing, Zhejiang 312000, China Zhejiang Collaborative Innovation Center for Prevention and Control of Mountain Geologic Hazards, Shaoxing, Zhejiang 312000, China School of Resources and Civil Engineering, Northeastern University, Shenyang, Liaoning 110819, China


Introduction
Several mined-out areas will be left behind in metal mines when using room and pillar mining. is is mainly from the need to reserve enough ore pillars to ensure safe mining activities in these processes. ese pillars and goaves constitute a large-scale underground goaf group. When mining is complete, the remaining pillars in the goaf group often face difficult problems for many reasons.
First, the underground geological conditions are complex, and there are several faults. When these faults pass through the area where pillars are, the stability of the surrounding rock mass becomes relatively poor. For example, Sun [1] studied the influence of faults on coal seam mining through simulation experiments. Second, the stress state of the surrounding rock changes when the ore body in the mine room is removed. e surrounding rock becomes extremely unstable from its original stress balance state and is prone to collapse. For example, Zhi et al. [2] analyzed the impact of coal mining on the farmland using the probability integral method. ird, there are often small structural surfaces in pillars, which are widely developed and will cut and damage them. For example, Liu and Zhang [3] adopted numerical simulations to analyze the instability of coal and rock roadways. In addition, some ore bodies may produce branched ore bodies with interlayers due to the particularity of their occurrence conditions. ese branches and interlayers often have weak layers, which aggravate the difficulty of mining pillars. For example, Hu et al. [4] studied the influence of weak layers on rock pressure and cracks through numerical simulation and model testing.
Mining pillars is therefore a relatively difficult and important task. An improperly selected mining approach for ore pillar can readily cause damage to the surrounding rock, which could cause safety accidents. For example, Zhao et al. [5] used ANSYS to analyze the impact of three different mining methods on the failure of rock mass using the trenchless, horizontal slice landfill method and the digging method. Ren [6] used the probability integral method to predict and analyze the displacement and deformation of the well body caused by the mining process. A suitable pillar extraction plan determines the safety of a mining enterprise. e optimization of the pillar mining plan is the first and necessary approach before pillar mining.
Numerical simulations have obvious advantages when optimization of pillar mining schemes. Using numerical simulations provides the stress field and ground pressure evolution law to provide the failure area and range of the surrounding rock and pillar. Scholars have performed abundant research in this area, such as using numerical simulations to study ground pressure. Wu et al. [7] established a numerical simulation model to study the evolution law of mining stress from footwall to fault. Qiao et al. [8] used FLAC3D to study the evolution of surrounding rock pressure in deep mining. Sun et al. [9] analyzed the changing law of the vertical stress field through numerical simulation. Guo and Lang [10] used UDEC to study the evolution of rock mass fractures at different depths. Singh et al. [11] estimated the stress on the rib with numerical modeling in the FLAC 3D software. Kumar et al. [12] studied these observed facts about the warning value of an AWTT in the laboratory using FLAC 3D. e above research shows that numerical simulations have good effects when studying the evolution law of the ground pressure, and can be used to consider the displacement law of coal seams and evolution law of stress.
Numerical simulations have been applied to study the damaged area and range of pillars and the surrounding rock, and Yang et al. [13] used FLAC3D to study the failure characteristics of the plastic zone under different coal pillar widths during the mining process. Yang and Company [14] used FLAC3D to numerically analyze the force characteristics, deformation, and failure range of coal pillars. Cao and Co [15] used numerical simulation to analyze the deformation and failure of surrounding rocks with different advancing degrees of the working face. Shen-Feng et al. [16] used FLAC3D to analyze the influence of section coal pillars on the stability of surrounding rock. Zhang et al. [17] analyzed the damage evolution process of tunnel surrounding rock based on the discrete element method. Xu et al. [18] discussed the surrounding stress conditions and deformation failure modes of coal pillar roadways through numerical calculations.
e above research shows that numerical simulations can help study failure characteristics, failure range, and deformations in surrounding rock. Both finite and discrete element approaches are suitable to study the deformation and failure of pillars and surrounding rocks.
Numerical simulations have been also used to study the mining sequence of pillars, and Shan et al. [19] used Map3d boundary element simulation software to study the two vertical and four horizontal mining sequences of ore body backfill mining. Huang et al. [20] used FLAC3D to study the influence of different mining sequences on the stability of coal pillars. Deng et al. [21] took the mining sequence as the research object and used numerical simulation to analyze the stress distribution characteristics under conditions of center-to-two-wing mining and two-wing-to-center mining. Ji et al. [22] used PHASE2D to study the overburden failure of inclined ore bodies under different mining sequences. Zhang et al. [23] used PHASE2D to analyze the change law of the stress and displacement of the surrounding rock under different mining sequences. Cai et al. [24] simulated and analyzed three feasible mining sequence schemes for mine slopes. Sears et al. [25] used numerical simulations to compare the mining sequence and the changes in the in situ stress field. e above research shows that numerical simulations have obvious advantages when studying the mining sequence of ore pillars. Among them, finite element analyses are widely used, such as the FLAC3D software, but some other software is also used, such as Map3d and PHASE2D.
COMSOL Multiphysics is a general-purpose software platform, based on advanced numerical methods to model and simulate physical field problems, which plays an important role in the analysis of geotechnical engineering applications, such as tunnels, excavations, surrounding rock stability, and support structures. Zhao [26] simulated the influence of mining disturbances and crustal stresses on rock deformation characteristics in COMSOL. Zhang et al. [27] analyzed stress in the surrounding rock mass around three chambers and the displacement change of several key monitoring points after excavation and evaluated the stability of the surrounding rock mass in COMSOL. erefore, the COMSOL software and finite element theory can be used to interpret the stress fields and mechanisms of surface collapse.
In summary, this paper first proposes two technical plans for the mining pillar to obtain optimal mining plans for the fourth middle section (+905 m level) of the Hongling leadzinc mine. Numerical simulation analysis was performed in the COMSOL Multiphysics software on two pillar mining plans. e analysis provided stress field and evolution law before and after stoping, as well as the influence on the pillar, roof, and the influence range. e optimal mining plan is finally determined in conjunction with a comprehensive analysis of factors, such as the degree of stress concentration, the stress state of the pillar, and the weak marble layer.

Engineering Background
Hongling lead-zinc ore is located in Chifeng City, China, as shown in Figure 1(a), which is the polymetallic mineral deposit shown in Figure 1(b). e length of the mining area is 5.7 km, the altitude ranges from 890 to 1362 m, and the vegetation is considered developmental. e mineral deposit has a length of 1350 m and is divided into two ore bodies, called 1# and 2#. e mineral deposit strike is north by east 60°, with the dip direction at the northwest and a dip angle ranging from 70°to 80°and surrounding rock of slate and skarn. us, the ground surface is amenable to collapse. e ore bodies are divided into eight sublevels at 50 m each. e upper four sublevels have already been mined and are recovering pillars. e fifth and sixth sublevels are mining rooms.
ere are large areas of collapse in the ground surface between the 1# and 2# ore bodies during pillar recovery, as shown in Figure 1(c). e collapse areas are close to industrial factory buildings with straight-line distances of less than 200 m, so there is a very high potential risk of collapse.
Due to the bifurcation of the 1# and 2# ore bodies, the weak marble layer will gradually be destroyed when mining the pillars, which causes the surface to collapse.
According to microseismic monitoring and analysis in the literature [28], the reasons for the collapsed rock mass are as follows. Failures are first produced within the +955 m level, which extended through the hanging wall to the ground surface, before moving to the footwall and forming the outline of the collapsed the ground surface. Meanwhile, multiple shear failures occur in the footwall surface, which extends downward in an arc-shaped path to the +905 m level to form the outline of the footwall destruction. In this process, many shear failures are produced in the ore body and hanging wall rock, which gradually extends to the surface. ese form the hanging wall destruction outline and ultimately cause collective large-scale failures. e footwall haul road in the +905 m level and the haulage way in the +955 m level already collapsed when the 5107 pillars recovered, as shown in Figures 1(d)-1(f ). e collapses not only influenced the mining scheduling but also caused high potential safety hazards.
According to the moment tensor analysis in the literature [28], it is concluded that the chain failure of the roadway is caused by the following reasons. Before recovery, tensile cracks existed in the footwall footing of the +855 m level and extended upward to the +905 m level, which did not cause surrounding rock damage. However, after the recovery of the 5107 pillars, the footwall haul road was damaged and collapsed due to cut-through cracks, which are composed of upward cracks and the newly generated tensile cracks in the mined-out roof region. ese cracks continued to extend upward and converged with the slanting shear cracks, which collectively formed a triangular failure in the footwall rock. Finally, the failures caused the tensile and shear cracks in the haulage way to extend and connect, which formed the haulage way chain failure.
When mining pillars, the rock mass around them is constantly deformed and damaged, which has caused surface and roadway collapse, and is very dangerous. For the overall stoping of the fourth middle section (+905 m level) pillars, selecting the optimal pillar extraction plan to recover the maximum number of residual pillars while maintaining the stability of the surface and surrounding rock is an urgent problem in mines. erefore, correctly choosing the mining plan of the pillar is important.

Proposed Technical Plan for Pillar
Mining. Two pillar mining plans are proposed based on the occurrence characteristics of pillars in the fourth middle section (+905 m level). First, the ore body is divided into four mining panels, namely, A, B, C, and D. Each panel contains 4-5 mined-out areas, and ore pillar mining is performed based on the panelbased mining method. e two proposed mining plans are as follows: Mining plan 1: From the middle of the ore body to the two wing retreat type stopes. e plan is to retreat from the position of the 13th exploration line in the middle of the ore body to the two wings, one panel at a time with a stoping sequence of C-B-D-A.
Mining plan 2: Backward mining from one side of the ore body to the other. e plan is to stope the ore body from the position close to the 0th exploration line to the 25th exploration line with a stoping sequence of A-B-C-D. e ore body model diagram of the fourth middle section is shown in Figure 2(d).

Establishment of the Mechanical Model of the Goaf Group.
e CAD software is used to establish the surface contour map and the three-dimensional ore body model map based on the geological data, vertical exploration line map, and middle section plan provided by the mine, as shown in Figure 3. e boundary size of the calculated model is 1400 m long and 1000 m wide. e surface is an irregular mountain with the highest altitude of 493 m and the lowest of 217 m. e constraint conditions of the calculation model are the surrounding rolling support and fixed bottom surface constraints. e model is a free surface and is only subject to its own weight. e resulting surface contour map is shown in Figure 2(a). ere are mountains on the upper part of the model with an elevation difference that fluctuates to reflect their overlying shape. e three-dimensional model of the ore body is shown in Figure 2(b). e location of the marble interlayer is shown in Figure 2(c). e contours of the mine room, pillars, and roof are connected based on the midsection plan, while their angles and contours differ.

Mechanical Calculation Parameters and Strength eory.
e COMSOL Multiphysics software was selected as the mechanical calculation software, and numerical simulations were performed based on the proposed mining plans. e multistep mining method is adopted in the simulations, and the stress field was iterated during the calculations. A highconfiguration SGI parallel machine was selected for parallel calculation due to the large ore body model and the amount of calculation. In the calculations, the stress evolution law on the roof was analyzed before and after pillar stoping, as well as the affected pillar and range. e results of the stress analysis provide the optimal stoping plan. e mechanical parameters of the rock mass are shown in Table 1.
e Mohr-Coulomb strength and yield criterion are used to judge the rock mass failure. e Mohr-Coulomb yield criterion is expressed in the form of principal stress as In the formula, σ 1 and σ 3 are the maximum and minimum principal stresses, respectively, and c and φ are the cohesion and internal friction angle, respectively. When f s > 0, the rock mass undergoes shear failure. e determination of the tensile failure of rock mass is determined from the relationship between the tensile stress σ t and the tensile strength R t When σ t ≥ R t , F ≥ 0, the rock mass is damaged under tension.

Comparative Analysis of Stoping Schemes
Based on Stress Evolution Law  Table 2.
In mining panel C, the left roof of the 4104 goaf was under tension, and the stress state changed from compressive at −1.75 MPa to tensile at 0.74 MPa with an increase of 142%. In mining panel B, the stress state of the right pillar of the 4100 goaf changed from compression to tension and increased 126% from −1.75 MPa to 0.46 MPa. In mining panels D and A, there would have stress effects on the surrounding rock.

Analysis of the Evolution Law of the Maximum Principal Stress.
e maximum principal stress curves in the middle section on mining plan 1 are shown in Figure 4. e negative values on the curves represent compressive stress, and the positive values represent tensile stress. e stress concentration area, stress change, increase, and increase rate of each panel during mining are shown in Table 3.
For mining panel C, the roof of the 4104 goaf was under pressure, and the compressive stress increased by 145% from −6.6 MPa to −16.2 MPa. e stress influence range on the left side of the panel was 18 m, which affected the roof of the 4104 goaf, and the stress influence range on the right side of the panel was 38 m, which affected the left pillar and part of the roof of the 4109 goaf, as shown in Figure 4(a). For mining panel B, the compressive stress of the roof of the 4104 goaf was reduced, which reduced by 74% from −16.2 MPa to −4.2 MPa. e compressive stress of the right pillar in the 4100 goaf was increased, which increased by 167% from −3.6 MPa to −9.6 MPa. e stress influence range on the left side of the panel was 43 m, which affected the right pillar and part of the roof of the 4100 goaf, as shown in Figure 4 Table 4.
For mining panel A, the compressive stress of the roof of the 4101 goaf was reduced by 162% from −1.62 MPa to 0 MPa. After pillars of the panel were mined, the compressive stress was released, and the tensile stress value becomes almost zero, as shown in Figure 5 The marble interlayer       e stress concentration area, stress change, increase, and increase rate of each panel during mining are shown in Table 5.
For mining panel A, the compressive stress of the roof of the 4101 goaf increased by 47% from −9.       goaf was most obvious. e compressive stress value reduced by 60% from −9.6 MPa to −3.84 MPa, and the compressive stress value of the right pillar in the 4109 goaf reduced by 52% from −7.95 MPa to −3.8 MPa. is showed that after the stoping of the middle section, the high-stress concentration on the left and right pillars of the 4109 goaf had been released, effectively alleviating the degree of stress concentration on pillars. Mining panel D would have a certain stress effect on the surrounding rock.

Optimization of Mining Plan Based on Stress Evolution
Law. is section summarizes the impacted pillars, impacted roofs, and stress range of mining plans 1 and 2, as shown in Table 6. According to mining plan 1, the panel C was first mined, and the stress influence range of 18 m and 38 m was generated on both sides of the panel, respectively. en, panel B was mined, and panel A had a stress influence range of 43 m. Mining panel D for the third time, there was only a stress effect on the surrounding rock. Panel A was finally mined and also only affects the surrounding rock. According to mining plan 2, the first mining area was panel A, which produced a stress influence range of 25 m on panel B. en, panel B was mined, which produced a 26 m stress influence range on panel C. Mining panel C for the third time, there was produced a stress influence range of 38 m on panel D. Finally, panel D was mined and only had a stress effect on the surrounding rock compared with mining plans 1 and 2. e advantage of mining plan 1 was that the impact of the pillars, roofs, and the stress impact range were all in panels C and B, which would be the first to be mined. Mining in panels D and A in the mining plan 1 only had an impact on the surrounding rock and had a small impact on the middle section. erefore, from the analysis of the stress evolution law, mining plan 1 had certain advantages over mining plan 2.

Optimization of Mining Plan Based on the Degree of Stress
Concentration.
e location and scope of the stress impact of the two mining plans on the pillars differ. An analysis was performed for the two stoping plans from the perspective of stress concentration. e minimum principal stress curve on the roof when the fourth middle section of the pillar was not mined is shown in Figure 7. is indicates there were two high-stress concentrations in panel C, one high-stress concentration in panel D, and two small-stress concentrations in panels A and D.
Stoping in mining plan 1 followed the sequence of C-B-D-A. Panel C was mined first, and the two high-stress concentrations were released, which effectively reduced the stress on pillars and roofs. en, panel B was stoped, and a high-stress concentration in the middle section was   removed. After panels D and A were mined, 4 weaker stress concentration areas were released. Mining plan 1 realized sequential mining and released high-stress concentrations in the middle section, which gave a better stoping plan. Stoping in mining plan 2 followed the sequence of A-B-C-D. Panel A was mined first, which had the least obvious stress concentration. e recovery of the panel A will not only alleviate the stress concentration, but also aggravate the highstress concentration of the panel B. Mining panel B could alleviate the high-stress concentration in panel B. However, after recovery, the stress concentration in the panel C increased, which could alleviate the stress concentration in panel C and impact stress concentration on panel D. When mining panel D, it had an impact on the surrounding rock. erefore, mining plan 2 cannot effectively reduce the high-stress concentration in panels B and C, but can aggravate the stress concentration. erefore, this mining plan was not the optimal stoping plan. Based on the above analysis, the optimal plan for mining the fourth middle section from the perspective of stress concentration is mining plan 1.

Optimization of Stoping Scheme
Based on the Force State of the Pillar. Reference [29] and Figure 8 show that there was stress concentration in the hanging wall of the 4200 goaf of the 2# ore body, and the surrounding rock was prone to damage. e right pillar in the 4200 goaf and the 4201 goaf would undergo X-shaped conjugate shear failure. e damage first occurred on both sides of the pillar and expanded to the middle area of the pillar in a circular arc shape, causing the destruction of the entire pillar. e 1# ore body was relatively narrow and thin, and the stress concentration during the stoping process was not obvious. Near the 13th exploration line, due to the existence of weak marble layers, the stress concentration on the surrounding rock was highest, and the surrounding rock was prone to damage. erefore, priority should be given to mining from the vicinity of the 13th exploration line.
Reference [28] and Figure 9 show that after mining fourth middle section, the surrounding rock of the hanging wall in the 4108 goaf suffered pressure failure. Pillars in the 4108 goaf bore compressive and shear forces and were the main pressure-bearing pillars in the middle section. e left pillar of the 4106 goaf and the 4107 goaf and the right pillar of the 4109 goaf were the secondary confined pillars in the middle section. When all the mining rooms were completed, the stress concentration on these pillars was particularly obvious, and they should be mined first.
Based on the above analysis, mining should be given priority to the 13th exploration line attachment from the perspective of the strength of the pillar, and mining plan 1 is the optimal stoping plan.

Optimization of the Mining Plan considering Weak Marble
Layers. It was known that there are weak marble layers between the 1# and 2# ore bodies. Large-scale subsidence occurred on the surface in the bifurcation area of 1# and 2# ore bodies based on previous work [30] shown in Figure 10. Figure 10 shows the minimum principal stress cloud nephogram of the stoping for fourth middle section. e marble interlayer hanged in the air after roofs and pillars of the +905 m level were recovered, caused damage to the ground surface, and resulted in a large amount of tension-shear stresses. Meanwhile, partial damage appeared in the top goaves of orebodies 1# and 2# and the plastic damaged areas in the +855 m level were extending.
If the marble interlayer area was mined first, it easily aggravated the stress concentration of the rock mass at the  interlayer. As a result, the collapse of the ground surface would intensify, resulting in more serious collapse accidents. erefore, mining panel A first would cause the continuous collapse of the surface in the interlayer area. ese collapses would be detrimental to the subsequent recovery of panels B, C, and D, which would also pose a significant safety hazard.
us, mining plan 1 is the optimal stoping plan from the perspective of weak marble layers.

Discussion
is article first proposes two technical plans mine pillars. Mining plan 1 retreats from the middle of the ore body to the two wings, while mining plan 2 retreats from one side of the ore body to the other. en, the article uses numerical simulations to analyze the two mining plans from four aspects: the stress evolution law, the degree of stress concentration, the stress state of ore pillar, and weak marble layer, and determines that the optimal mining scheme is plan 1. e determination of the optimal mining plan provides mechanical analysis and safety guidance for the pillar mining of Hongling lead-zinc mine, which can effectively reduce the hidden safety hazards during the mining of the pillar, and has important engineering significance. At the same time, the stress evolution analysis method provided in the article can provide technical guidance and theoretical reference for the safe  mining of pillars of the same type of complex goaf groups. irdly, the optimization of the multi-angle mining plans provided in the article can obtain more accurate analysis results, which ensures the optimality and accuracy of the selection of the mining plan. e two plans were developed by dividing panels during mining and did not specifically address the mining procedures or forms within a single room. is work does not consider technical and economic indicators such as the organization of process operations and the ability to recover residual ore. erefore, the optimal stoping plan is based only on the analysis of ground pressure evolution law and surrounding rock failure.
As the ore body was too long, for numerical simulations, the established mechanical calculation model had a boundary length of 1400 m and a width of 1000 m. Only three lithological rock masses such as ore body, surrounding rock, and weak interbed were established. e feature where the strength of the same lithological rock mass changes at different locations was not considered. At the same time, other small faults in the ore body were ignored and simplified. e establishment of each pillar in the mechanical calculation model was based on the site plans and sections. erefore, the spatial shape of pillars for each mine room differs to restore the actual outline of the goaf on-site as much as possible.

Conclusion
is article proposes two technical plans for mine pillars and analyzed them from four perspectives such as the law of stress evolution, the degree of stress concentration, the state of the pillar, and the weak layer of marble. Each analysis determined that mining plan 1 was the optimal plan. us, mining plan 1 is the optimal stoping plan.
According to the mining plan 1, the high-stress concentration area in the middle section could be successively released, which effectively alleviated the surrounding rock stress. It may also avoid aggravating the collapse of the marble interbedded area, which would happen with mining plan 2. e method of analyzing the mining plan from multiple angles provided in the article can obtain more accurate analysis results, which ensures the optimality and accuracy of the selection of the mining plan. It provides mechanical analysis and safety guidance for pillar mining and has important engineering significance. At the same time, it can also provide technical guidance and theoretical reference for the safe mining of ore pillars of the same type of complex goaf groups.
Data Availability e data that support the findings of this study are available from the corresponding author, G Hu, upon reasonable request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.  [30].