Wawu IV ore block belongs to medium and high mountain area. After denudation and strong erosion, it forms the geomorphic characteristics of vertical and horizontal gullies and overlapping mountains. The highest altitude in the area is 1710 m, which is located in Gaobaling in the east of the mining area, and the lowest is 574 m, which is located in Zhuyuan River in the southwest corner of the mining section. The maximum relative elevation difference is 1136 m, generally about 400∼800 m. Most of the peaks are round, a few isolated peaks are bamboo shoots, and most of them are cliffs. The valley is a “V” shaped valley with narrow valley bottom and valley slope of 40°∼80°. The whole mining area is covered with lush forests, and the steep terrain is conducive to drainage.

Phosphate rock is located in the northeast of Huangling anticline. The exposed strata in the area include Mesoproterozoic Shennongjia group (Ptsh) and Kongling group (Ptkn). Neoproterozoic upper Nanhua NANTUO formation (Nh₂n), Sinian Doushantuo Formation (Z₁d), Dengying Formation (Z₂dn), Paleozoic Cambrian (∈), Ordovician (O), Silurian (S) and Cenozoic Quaternary. The phosphorous rock series in this ore block is mainly composed of phosphorous mudstone and dolomite. The main ore bed in the ore section is Ph₁ ³ ore bed, which is stably distributed in the north of the ore section and occurs in layers. It is monoclinic structure as a whole with slight wavy fluctuation. The ore bed is inclined to the southeast as a whole, with a dip of 130°∼175°, and the dip angle is generally 18°∼35°. Affected by faults, the ore bed becomes steeper or even upright locally.

The mining area belongs to subtropical monsoon climate with four distinct seasons and abundant rainfall, with typical mountain climate characteristics. The annual maximum precipitation (83 years) is 1552.42 mm, the annual minimum precipitation (88 years) is 913.7 mm, the annual average rainfall day is 155 days, and the maximum daily precipitation (September 9, 83 years) is 109.7 mm. The precipitation is concentrated. The precipitation in rainy season (may October) accounts for 80.9% of the annual precipitation, the precipitation in dry season (November February of the next year) accounts for 9.4% of the annual precipitation, and the freezing period is from December to march of the next year.

The phosphate rock layers in ore block IV of Wawu mine are concealed underground, and more than two-thirds of the phosphate rock layers occur above the local erosion base level. The hydrogeological conditions of ore block IV belong to the karst water filled deposit of “the water filled rock stratum is mainly composed of dissolution fractures, the roof has direct and indirect water inflow, and the hydrogeological conditions are medium to complex.”

The surrounding rock of phosphorus ore bed roof in Wawu IV ore block is dolomite, with hard rock, high strength and good integrity. Ph₁ ³ phosphorene bed is dense, hard and complete. The engineering geological conditions of Wawu IV ore block belong to medium type. Because there are some differences between the mechanical parameters of rock and rock mass, it is necessary to weaken the mechanical parameters of rock to obtain the mechanical parameters of rock mass. Firstly, the mechanical parameters of rock are obtained through indoor experiments (Wu et al., 2020; Dong et al., 2021; Gao et al., 2021; Wang and Liu 2021), and then weakened by using empirical reduction method and Hoek-Brown strength criterion (Serrano et al., 2016; Shi 2016) according to the on-site geological conditions.

The selection of rock mass parameters in this paper is based on the analysis of the existing data of the mine, and referring to the above relevant literature and field mechanical tests, to comprehensively obtain the physical and mechanical parameters of roof dolomite, floor dolomite, filling body and ore deposit. The obtained physical and mechanical parameters of rock mass are shown in the Table 1.

According to records, earthquakes in Yichang City in history are mainly weak earthquakes, and there has been no destructive earthquake. According to the seismic ground motion parameter zoning map of China (GB18036-2001), the mining area belongs to the seismic basic intensity VI area, and the seismic basic acceleration is 0.05 g. According to the crustal stability grade and seismic index, the mining area belongs to the stable area. No magmatism has been found in the ore section. The metamorphism in the ore block only shows the shallow metamorphism of the basement strata of the Middle Proterozoic Shennongjia group.

In view of the complexity of stope mining technology, in order to facilitate modeling and calculation analysis, the following simplifications and assumptions are made:

(1) It is assumed that the rocks and fillings in the model are homogeneous and isotropic continuous media.

Select the drift of ore Section 4 of Wawu phosphate mine, take 32 m thick rock mass above and below the phosphate rock layer, and select the section passing through exploration line 23 for numerical analysis. The phosphate rock layer is subjected to the self-weight load of the overlying rock mass. The self-weight load is estimated by multiplying the thickness of the phosphate rock layer from the surface by the unit weight of the soil. The thickness from the surface is taken as 400 m, and the unit weight of the rock mass is taken as 27 KN/m³, equivalent to 400 m × 27 KN/m³ = 12420 KPa, applied to the upper part of phosphate rock layer, p = 12420 KPa. The selected simulation area is shown in Figure 1 and Figure 2.

The stope roof can be assumed to be simply supported beams at both ends Yang et al. (2015). According to the material mechanics, the stress at any point on the upper and lower surfaces of the neutral axis of the rock beam is: σ ( x ) = γ ( 2 x − l ) sin ⁡ α 2 ± γ x ( x − l ) cos ⁡ α h (1) Where, ɑ is the dip angle of ore body, (°). l is the span of rock beam, m. h is the height of rock beam, m. γ is the unit weight of rock mass, 10⁴ N/m³. The maximum tensile stress occurs on the lower surface of the neutral axis of the rock beam a x = L / 2 + h tan α / 6 and the maximum tensile stress is: σ max = 3 γ L 2 ⁡ cos ⁡ α 4 h − h γ ⁡ tan 2 ⁡ α ⁡ cos ⁡ α 12 (2)

The maximum allowable span of roof inclination is: L qy = [ 4 h σ t 3 γ ⁡ cos ⁡ α − h 2 ⁡ tan 2 ⁡ α 9 ] 1 / 2 (3)

The maximum allowable span of the vertical strike of the roof is: L s p = L q y | ( α = 0 0 ) = [ 4 h σ t 3 γ ] 1 / 2 (4)

The stress analysis of simply supported beam is shown in Figure 3.

For the stope mined by stage open stope and subsequent filling method, the width of the stope can be determined according to the safety guarantee and production capacity of stope mining Zhang et al. (2021). When the stope safety is guaranteed, the greater the width and height of the stope, the greater the production capacity of the stope Zhang et al. (2014). The calculation formula of room width is as follows: L = 4.5 ( β γ Η + σ t ) γ Η ( 1 − 0.8 β ) − 0.8 σ ′ t (5) Where, H is the mining depth of the ore body, m. β Is the lateral pressure coefficient, generally taken as 0.65. σ _t T is the tensile strength of top pillar rock stratum (reduced according to the tensile strength of rock mass, generally 0.4), MPa. Taking H = 400 m, γ = 29 kN/m³, σ _t = 2.1 MPa is substituted into Eq. 5, It can be obtained that the room width of the filling yard of phosphate rock in Ore Section 4 of Wawu phosphate rock area is about L = 9 m. Therefore, we have obtained the empty yard span of each scheme, as shown in Table 2.

The calculation result of simply supported beam theory is the maximum allowable span under corresponding ore and rock conditions, so the result is too large. In practical application, the calculation result of reasonable room width should be taken as reference, and the maximum cannot exceed the calculation result of simply supported beam theory. Therefore, the width of the simulated ore room is 12 and 9 m respectively.

PLAXIS program has strong application and can simulate complex engineering geological conditions, especially suitable for deformation and stability analysis. PLAXIS program can calculate two kinds of engineering problems: plane strain problem and axisymmetric problem. It can simulate soil, wall, slab, beam structure, contact surface between various elements and soil, anchor bolt, geotextile, tunnel, pile foundation and other elements. Analyze 2D and 3D deformation, consolidation, graded loading, stability analysis and seepage calculation in geotechnical engineering. PLAXIS uses the finite element strength reduction method to calculate the safety factor in the sense of soil mechanics.

The values of displacement, stress, strain and member internal force can be output to the table. The constitutive model of soil, linear elasticity, Mohr Coulomb, soft soil model, hardening model and soft soil rheological model can simulate the construction steps, and the post-processing of multi-step calculation is simple and convenient (Huang et al., 2017; Zhu et al., 2021b).

The lithology of the phosphate rock is mainly composed of dolomite and mudstone, which are elastic-plastic bodies. Therefore, it is applicable to the strength failure criterion of Mohr Coulomb. The mechanical model is shown in Figure 4. Based on Mohr-Coulomb strength criterion, the shear failure envelope ƒ^s = 0 on segment AB is defined as: f s = σ 1 − σ 3 1 + sin ⁡ ϕ 1 − sin ⁡ ϕ + 2 c 1 + sin ⁡ ϕ 1 − sin ⁡ ϕ (6)

Based on the tensile failure criterion, the tensile failure envelope ƒ ^t = 0 on segment BC is defined as: f t = σ t − σ 3 (7) Where, σ ₁ is the maximum principal stress, MPa. σ ₃ is the minimum principal stress, MPa. φ is the internal friction angle, (°). c is cohesion, MPa. σ _t is tensile strength, MPa. In Figure 4, the points on line AB or line BC are in critical failure state. The lower left half (area I) of the area surrounded by three line segments AB, BC and σ ₃-σ ₁ = 0 is in an undamaged state, the upper right half (area II) is in a damaged state.

This study mainly focuses on the numerical simulation and optimization of stope structural parameters Zhang and Chen (2016). According to 3.3 simply supported beam theory and reasonable chamber width calculation results, the maximum allowable span and reasonable span of each case are determined. Combined with the existing mining technical conditions and actual production conditions, the following four cases are put forward, and the stope stability under the four cases is compared. There are two options. The first scheme is to change the stope structural parameters (Case 1 and Case 3), and the second scheme is to increase waste rock filling (Case 2 and Case 4).

Case 1: The stope structural parameters are the pillar width of 5 m, the room width of 12 m and the ore body thickness of 8 m. No filling is carried out after one-time excavation.

Case 2: The stope structural parameters are pillar width of 5 m, room width of 12 m and ore body thickness of 8 m. For waste rock filling, three ore chambers are filled horizontally and longitudinally to form a large point column with the ore pillar. Two ore chambers are left at the periphery without filling, and then filled all the time according to the same method, as shown in Figure 5.

Case 3: The stope structural parameters are pillar width of 7 m, room width of 9 m and ore body thickness of 8 m. No filling is carried out after one-time excavation.

Case 4: The stope structural parameters are pillar width 7 m, room width 9 m and ore body thickness 8 m. The waste rock filling method is the same as Case 2. The schematic diagram of mining scheme is shown in Figure 6.

The finite difference software PLAXIS 2D is used for modeling. According to the elastic-plastic theory, the mining influence range after stope excavation is usually within 5 times of the roadway radius. In this paper, the model is established based on 3–5 times of the analysis area. The model is 270 m long, 122 m high and 8 m thick, with 42,763 nodes and 5,275 units. The horizontal displacement constraint in the X direction is adopted on the left and right sides of the model, the horizontal and vertical constraints are applied on the bottom of the model, the model surface is free constraint, show in Figure 7.

Because the original stress balance state is broken, the roof and floor of the stope are affected by horizontal stress, so they squeeze towards the goaf, resulting in displacement towards the goaf. Among them, the vertical displacement generated in the middle of the roof of the goaf is the largest. It can be seen from Figure 8 that the maximum displacement generated by the roof in Case 1 is 2.8 cm, in Case 2 is 2.0 cm, in Case 3 is 1.8 cm and the maximum displacement of the top plate in Case 4 is 1.0 cm. Therefore, the maximum displacement generated by the roof in Case 4 is the minimum. After filling the goaf, it can be seen from the Figure 8 that the overall displacement of the rock mass in the mining area decreases to a certain extent. The addition of filling body reduces the settlement displacement of surrounding rock in the mining area and the affected area of surrounding rock settlement, so that the influence of mining disturbance on the overlying strata of goaf roof and surface becomes smaller. This shows that filling the goaf after room mining effectively enhances the stability of the seam and reduces the possibility of geological disasters such as surface subsidence and mine collapse.

Due to the mining of stope, the stress of stope and surrounding rock is released, breaking the original stress state and redistributing the stress of surrounding rock. Stress concentration occurs at the top and bottom corners of the stope, and the compressive stress is mainly concentrated on the room and pillar. It can be seen from Figure 9 that the maximum principal stress of the whole goaf roof in Case 1 is about 36.84 MPa, in Case 2 is about 32.51 MPa, in Case 3 is about 28.98 MPa and in Case 4 is about 24.71 MPa. Therefore, the maximum principal stress of goaf roof in Case 4 is the minimum. After filling the goaf, the distribution of stress concentration area becomes less, and part of the tensile area is transformed into compression area; the maximum compressive stress of ore body and the maximum tensile stress of ore pillar are obviously reduced. This shows that the filling body bears part of the vertical compressive stress instead of the roof, floor and pillar, so as to release the concentrated stress in the goaf, effectively prevent the vertical deformation of the roof and floor of the goaf, reduce the possibility of pillar failure and collapse, and strengthen the stability of the surrounding rock of the seam.

From Figure 10, in Case 1, the plastic zone of the surrounding rock around the stope is asymmetric, mainly because the shear stress generated by the surrounding rock is greater than the shear strength of the rock mass, resulting in the plastic state of the surrounding rock. In Case 3, the plastic area of the stope is basically concentrated at the pillar. And the plastic area generated by Case 1 is larger than that generated by Case 3, and the stress of surrounding rock changes sharply. This is because the two ends of the stope have exceeded the mechanical strength of the rock mass itself, resulting in stress release. To sum up, the stope under the two working conditions are unstable. The plastic zone ratios for Case 1 and Case 3 are shown in Table 4.

In Case 2 and 4, the plastic area of the stope is basically concentrated at the top and bottom of the pillar, and the plastic area is small. Compared with the working condition of unfilled goaf, the range of plastic zone of phosphate rock layer is reduced to a certain extent, and the rock mass stress around the plastic zone is in the state of elastic stress, and the overall stability of phosphate rock layer is strengthened. It can be seen from Figure 11 that the plastic area generated by Case 2 is larger than that generated by Case 4. Compared with Case 1 and 3, the plastic area of Case 2 and 4 is smaller and the stope stability is higher. The plastic zone ratios for Case 2 and Case 4 are shown in Table 5.

The results of numerical simulation show that the maximum values of tensile stress and final displacement of each scheme usually appear in the middle of the area, and the increase of mine room span is very easy to cause the rapid increase of the maximum settlement in the mining area. It can be seen from the plastic zone that with the increase of mine room span, the plastic zone is also increasing. In Case 1 and Case 3, the failure area in the pillar appears the failure state of plastic zone penetration, so it is not suitable to use. Considering several factors such as roof stability and plastic area, Case 2 and Case 4 can ensure the safety of stope mining. The data of Case 2 and Case 4 show little difference. Considering the efficient production and stope stability, Case 4 is the most suitable, that is, the stope structure scheme of 7 m pillar and 9 m room waste rock filling.

The optimized stope structure parameters are applied to the stope with a buried depth of about 400 m. The edge pillar, stope pillar, roof and floor can remain stable, and only a small amount of ore caving occurs at the place with relatively developed joint fissures. Compared with the original stope design scheme, the optimized scheme increases the width of pillar and reduces the number of pillars, The stability of the stope has been greatly improved, and the recovery rate of the mine has been increased to about 79%, which has created good economic benefits for the mine.

Stope structure parameter design is a multi-factor and multi-objective optimization problem. The more factors and objectives, the more advantages of this method can be fully reflected. Therefore, it can be regarded as a general method to solve this kind of problem. The evaluation indexes set in this paper only include stability indexes and economic indexes, while other indexes in the actual production of the mine (mining and cutting ratio, stope production capacity, geological conditions, etc.) can also be taken into account; At the same time, the comprehensive evaluation method used in this paper is only the objective weighting method. If the objective weighting method is combined with the subjective weighting method, a more accurate comprehensive evaluation value can be obtained, and the stope structure parameters most in line with the current situation of mine production can be obtained accordingly.

Mine geological conditions are generally complex, and the simulation process of finite element analysis software fully reflects the complex and changeable geological conditions of the mine. The determination of ore body parameters should adopt the combination of numerical simulation and theoretical calculation. Therefore, based on the elastic-plastic constitutive model of Mohr Coulomb strength criterion, the reasonable span and critical span are calculated by using the simply supported beam theory and the calculation formula of room width. By improving the reliability of the calculation results, the data are closer to the engineering practice, so as to ensure the safety and efficiency of mining operation.

The generalization of geometric model involves complex geological bodies. This paper only studies and analyzes the two-dimensional model of phosphate rock mining without three-dimensional numerical simulation, which is slightly different from the actual situation. The model itself has experienced a generalization, so various assumptions are put forward in the simulation calculation to facilitate the simulation calculation. For example, the horizontal displacement constraints in the X direction are adopted on the left and right sides of the model boundary conditions, the horizontal and vertical constraints are applied on the bottom of the model, and the model surface is free constraints, which will be different from the actual situation of the mine; In the actual engineering filling, the filling body can not completely fill the ore chamber, and this factor is not considered in the simulation calculation. The above deficiencies need to be further discussed in the follow-up study.