Research on Optimization and Numerical Simulation of Layout Scheme of Mining Approach in Downward Slicing and Filling Method

Abstract

The stability of the filling roof—as an important bearing unit in the stope of the access stope in the filling mining method—is of great significance to guarantee the safe and efficient production of the mine. Arrangement of the mining approach in downward cemented filling stope is the key factor affecting the stability of the filling body roof. Based on a combination of laboratory tests, theoretical analysis, and numerical simulation, the influence of different mining approach arrangements on the stability of the filling body roof is analyzed. The weak filling surface is formed between adjacent mining paths. The mechanical strength of the weak filling surface is significantly reduced by laboratory experiments. The relationship between the distribution of the weak filling surface, azimuth angle, and the stability of the filling roof is further studied by numerical simulation. The results show that, when the upper and lower layered mining approaches are arranged vertically or nearly vertically, the areas of stress concentration and the plastic zone in the numerical simulation results are the smallest, and the stability of the filler roof is the best.

1. Introduction

The downward slicing and filling method was introduced to China in the 1960s [1,2]. Under the background of high-efficiency and green mining, it develops rapidly and is widely used in metal mines [3]. There are often weak structural planes such as bedding, joints, weak interlayers, and fracture zones in an underground engineering rock mass. The existence of a weak plane not only weakens the overall strength of rock mass but also is a potential risk factor for overall engineering instability.

The most common panel mining method in downward layered cemented fill mining is roadway mining. Ore bodies in the layered area are usually divided according to a certain width of the roadway and then filled by interval or continuous mining; however, no matter which mining method is adopted, due to the different filling sequences, there must be weak filling surfaces between adjacent mining paths. When the lower stratified orebody is mined, the filling weak surface is arranged in the filling roof parallel to the gravity direction. In the actual mining process of ore bodies, production technicians usually reinforce the support under the weak filling surface.

The occurrence of filling weak surfaces in exposed roofs is different from the layout of an upper- and lower-stratified mining approach. When the upper- and lower-stratified mining approaches are arranged vertically and staggeringly, the angle between the filling weak surface and the direction of the mining approach is a 90° location relationship; when the upper- and lower-stratified mining approaches are arranged parallel, the angle between the filling weak surface and the direction of the mining approach is a 180° location relationship. The research shows [4,5,6,7,8] that in the field of geotechnical engineering, the occurrence of weak surfaces in the engineering structure body is different, and the deterioration degree of the stability of the engineering structure body is also different. Even though existing research [9,10,11,12,13,14] has used numerical simulation to simulate and optimize the layout of upper and lower layers for mining access, the distributions of the stress and plastic zones in the two different mining and filling sequences were obtained through numerical simulations by Jiang N et al. [9]; the reasonable mining face width was determined combined with numerical simulation [10]; and field investigation and numerical simulation were combined to study the mechanisms of roof collapse by Zhong M et al. [11]. Research on the intrinsic mechanisms of the layouts of upper and lower layers is still lacking, which puts engineering scheme practice far ahead of scientific theory. Cao ZQ et al. [15] studied the relationship between the rectangular section size and stope stability using FLAC^3D based on practical engineering and analyzed the distribution of stresses, displacement, and the plastic zone.

For the direct shear tests on cemented paste backfill, Gao T et al. [16] prepared a two-layered cemented paste backfill (CPB) with layering angles of 5°, 10°, 15°, 20°, 25°, and complete CPB to conduct direct shear experiments. Koupouli et al. [17] performed direct shear tests on CPB–CPB interfaces in order to assess the frictional shear strength parameters. The action of compression–shear, measured using 30°, 45°, and 60° variable angle shear tests (VAST) was conducted by Chen T et al. [18] to investigate the mechanical properties. Researchers have studied the strength, deformation, and failure of the direct shear tests on CPB [19,20,21,22,23,24].

Under this engineering background, combined with the direct shear test of the weak surface of the filling body in the chamber, the stability of the filling roof under different layered mining approaches is numerically simulated and analyzed using a medium interface unit in FLAC^3D, the finite difference numerical simulation software. Finally, the optimal layered layout scheme of upper and lower layers is determined, which provides a theoretical basis for mines adopting the downward horizontal layered cemented filling method.

2. Materials and Methods

2.1. Laboratory Test on Mechanical Parameters of Filling Weak Surfaces

2.1.1. Sample Preparation of Layered Paste Filler

The test sample was layered full of tailings paste filler. The aggregate was taken from the whole tailings of a copper mine in Jiangxi Province. The tailings were placed in dry and ventilated places to dry and then fully compacted to no obvious coarse particles. The chemical mineral composition and proportions are shown in

Table 1. Chemical composition and proportion of Total Tailings.

Chemical Composition	SiO2	CaO	MgO	Al2O3	Fe	Cu	S	C	Others
Percentage of mass (%)	50.73	16.55	3.14	4.33	5.59	0.04	1.30	1.78	16.54

. P.O32.5 mineral Portland cement was used as cementing agent to prepare filling slurry, with a lime–sand ratio of 1:4 and a mass concentration of 74% according to the actual filling ratio of the mine; the density of the cemented paste backfill was 2.02 g/cm³.

In actual stope, the filling weak surface was arranged parallel to the direction of gravity, so the filling weak surface was designed to be distributed horizontally in the layered paste filling body. Sample pouring dimensions were 150 mm × 150 mm × 150 mm, the sample pouring process was divided into two steps, the first step was to pour half-height filler in the filling mold, and the second step was to pour the other half of the filler in the filling mold after 24 h of rest until the filler had the basic strength. When there was certain bond strength between weak surfaces of layered paste filler, the sample was de-molded to obtain layered paste filler with horizontally weak surface. According to the production environment of mining access and different layering times, the specimens were kept in a constant temperature and humidity curing box for 28 days, with a curing temperature of 20 °C and an air humidity of 41%.

2.1.2. Direct Shear Test of Layered Paste Filling

The direct shear test was conducted with RMT-150C rock mechanics experimental system, as shown in Figure 1, and the inner wall size of the direct shear test box was 150 mm × 150 mm × 150 mm. The direct shear specimen was placed in the shear box, and a height of 10 mm was reserved between the upper and lower shear boxes for shearing, as shown in Figure 2. The direct shear test adopts the force (large) displacement control method of the test system. According to the standard for test methods of engineering rock mass, the axial force is controlled by force stroke, and the loading rate is 500 N/s; shear force is controlled by displacement; and the loading rate is 5 × 10⁻⁵ m/s. A total of four groups of direct shear tests of filling weak surfaces under normal stress (0.5 MPa, 1.0 MPa, 1.5 MPa, 2.0 MPa) were conducted. The test loading process was vertical preload → horizontal preload → vertical load → horizontal load → specimen failure.

In the direct shear test, the position of the shear failure surface and the filling weak surface basically coincide, and the shear failure surface was neat and regular, and the shear failure surface was smooth and flat. It can be concluded that when the layered paste filler was subjected to large shear load as a whole, shear slip is most likely to occur at the weak filling surface. Layered paste filler had fewer other failure types near the weak surface during shear-slip process, which fully indicates that the mechanical strength of the weak surface is much lower than that of surrounding filler.

The shear stress–strain curve divides the shear failure process of filling weak surfaces into two stages, strain-softening stage and strain-hardening stage, as shown in Figure 3a. The peak shear stress of filling weak surface increases with the increase in normal load and the strain reaching the peak shear stress migrates backward. Based on Moore-Coulomb strength criterion, as shown in Figure 3b, the parameters of shear strength of filling weak surfaces are obtained by linear fitting with least squares method, with cohesion of 0.105 MPa and internal friction angle of 20.49°.

2.2. Research on Numerical Simulation of Direct Shear Test

2.2.1. Establishment of Numerical Model

The numerical model for direct shear test is cubic. The size of the model is set at 100 times the size of the test specimen in the laboratory. The dimension length is 15 m × 15 m × 15 m; the cross-sectional area of excavation is 3 m × 4 m. An interface unit is added to half the height of the numerical model to simulate the filling of weak surfaces by filling samples. The model is divided into 76,245 grid units and 78,693 nodes. Among them, the dimension of line segments on the interface boundary is 0.3 m. Other line segments are dimensioned by 0.3~0.5 linear gradient according to the distance from the interface. In Figure 4, the numerical model is shown.

A displacement constraint is added to the normal direction of each surface of the material group below the interface and a uniform load is applied to the top surface of the material group above the interface. At the same time, the shear rate of 5 × 10⁻³ m/s is used to control the sliding of material groups above the interface towards the x-axis in the positive half direction.

2.2.2. Calibration of Interface Meso-Parameters

The material units with the strongest mechanical properties on both sides of the interface are selected to calculate the equivalent stiffness of the interface. The calculation formulas are as follows:

$$k_{n0}=k_{s0}=10\times \left[\frac{(K+\frac{4}{3}G)}{\Delta z_{\min }}\right]$$

(1)

where K is the volume modulus of filling material; G is the shear modulus of filling material; and ∆_Zmin is the smallest grid size in the normal direction of both sides of interfacial, which equals 0.3 m. The calculation of two modulus parameters K and G are as follows:

$$K=\frac{E}{3(1-2v)}$$

(2)

$$G=\frac{E}{2(1+v)}$$

(3)

It can be calculated by combining the above formulas. Initial Normal Stiffness k_n0 and Initial tangential stiffness k_s0 are equal to 2.14 × 10¹⁰ Pa/m.

2.2.3. Calibration of Meso-Parameters

In order to determine reasonable stiffness parameters, the numerical simulation results are closest to the laboratory test results. Initial normal stiffness and shear stiffness of 0.001-, 0.01-, and 0.1-times are, respectively, set to simulate the shear behavior of interface changing with stiffness parameters under four normal stresses (0.5 MPa, 1.0 MPa, 1.5 MPa, 2.0 MPa).

The shear stress–shear strain curves during interfacial shear are similar to those of direct shear tests in laboratory. They can be divided into two stages, strain softening and strain hardening, without changing the stiffness parameters and normal stresses. Under the same normal stresses, the peak shear stress and peak shear strain at the interface show a non-linear upward trend with a decrease in stiffness parameters; at the same stiffness parameter, the peak shear stress and peak shear strain at the interface also show a non-linear upward trend with an increase in normal stress applied to the model. It is shown that the interfacial micro-stiffness parameter will affect the overall strength of the interfacial in the model but will not change the development trend of the interfacial shear stress–shear strain curve. The fitting curves of peak shear stress and normal stress under different stiffness parameters are shown in Figure 5. After adjusting, Pearson’s correlation coefficient R² = 1.

In order to comprehensively compare and analyze the numerical simulation test results under different stiffness parameters with the laboratory test results, the peak shear stress, cohesion, and internal friction angle of the test results are extracted and plotted as a pie-column diagram for comparison, as shown in Figure 6 and Figure 7.

Combining the above numerical simulation results for analysis, the initial normal stiffness, k_n0, and initial tangential stiffness, k_s0, of the decomposition surface are chosen to be 2.14 × 10⁷ Pa/m, for which the interfacial cohesion is 0.105 MPa and the internal friction angle is 20.49°. The results of numerical simulation direct shear tests are closest to those of laboratory tests, which can effectively simulate the mechanical characteristics of filling weak surfaces.

2.3. Optimum Simulation of Layout Scheme of Mining Approach

The different layout schemes of mining approach are mainly different from the occurrence and exposed area of weak surface of upper-stratified filling relative to lower-stratified filling approach. The stability of filling roofs is different under different layout schemes of mining approaches.

Establishment of Numerical Models for Different Mining Routes

Wang XJ et al. [25] fully considered the engineering practice—such as the load of the mining rock, the dip angle of the ore body, the staggered arrangement of the adjacent layered roadway, and the contact between the backfill and the surrounding rock—then used FLAC^3D to simulate the orthogonal calculation of roof stability under the influence of multiple factors.

In order to comprehensively analyze the stability of filling roofs under different layered mining approach arrangement schemes, the excavation schemes of the mining approach—without considering the filling weak surface and six layouts of the mining approach with 30°, 45°, 60°, 90°, and 180° between the filling weak surface and the central axis of the approach direction–are designed. The model grid unit is generated by mixing tetrahedron and trihedron units. In the vicinity of excavated tunnels and the filling roof of tunnels, the mesh size is controlled to 0.3 m, while in other locations, the mesh size is controlled to 0.5 m. The excavation models under different schemes are shown in Figure 8.

The filling weak surface in the model obeys the linear Coulomb shear strength criterion and other materials obey the Mohr-Coulomb strength criterion. The mechanical parameters of the model material are shown in

Table 2. Physical and mechanical parameters of the model.

Physical and Mechanical Parameters	Density (kg/m3)	Elastic Modulus (GPa)	Cohesion (MPa)	Internal Friction Angle (°)	Poisson’s Ratio	Tensile Strength (MPa)	Normal Stiffness (Pa/m)	Shear Stiffness (Pa/m)
CPB	2000	0.65	1.15	45	0.20	0.64	/	/
The weak surface	/	/	0.125	20	/	/	2.14 × 107	2.14 × 107
Rock mass	3500	40	22	36	0.30	16	/	/

.

3. Results

By reason of the weak filling surface being distributed in the filling roof at multiple angles, the mechanical state of the filling roof is studied by an interception section, which is located in the X–Y plane of the model coordinate axis and at the lower surface of the filling roof in the Z direction of the model height, when processing and analyzing the calculation model results.

Stress Field and Plastic Zone Distribution

Since the numerical model calculation is completed, the stress state nephogram and plastic zone distribution of the profile for the six mining route layout schemes can be shown in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14, where the purple area in the plastic zone distribution represents the tensile plastic zone, and the red area represents the shear plastic zone.

4. Discussion

By comparing and analyzing the stress state and plastic zone distribution of the filling roof after excavation of the mining approach, the optimum mechanical state of the filling roof after excavation of the mining approach is selected, excluding the situation without weak filling surface.

The optimal structure parameters were determined based on the calculation and analysis of displacement variation rule, surrounding rock stress distribution and plastic zone by Peng, G [26]. According to the calculation results of the model, the mechanical behavior characteristics of the backfill roof after excavation of the mining approach under different schemes are consistent, which is mainly reflected in the tensile stress concentrations and tension plastic zones of the roof directly above the mining approach, and the compressive stress concentrations and shear plastic zones of the two sides of the mining approach. Skrzypkowski K [27,28] presented some issues with laboratory and spatial numerical modeling of cemented paste backfill, which should be taken into account; thus, the above numerical simulation results combined with indoor experiments can infer that the existence of weak filling surface increases the distribution area and breadth of tension stress concentration area of the filling roof to a certain extent. When the angle between upper and lower layered mining approaches is arranged at sharp angles, the tension stress concentration area of a filling roof mostly appears near the weak filling surface, which has a certain partition effect on the tension stress concentration of the filling roof. When the upper and lower layered mining approach is arranged vertically, the size of the area of tensile stress concentration and the degree of change is the smallest, the distribution characteristics of the stress field of filling weak face have the least influence on the filling roof, and the distribution characteristics of stress field are closest to those of the plan without the filling weak face. The distribution of the plastic zone of the filling weak face has the least influence on the filling body roof, and the distribution of the plastic zone is the closest to that of the plan without the filling weak face.

By comparing and analyzing the stress field and plastic area distribution of the filling body roofs of the six schemes, it can be seen that the mechanical state of the filling body roof is closest to that of the non-filling weak surface scheme after excavation of the upper- and lower-layered mining approach; the area and degree of variation in the area of tensile stress concentration and the area of plastic zone concentration are minimal; and the mechanical state of the filling body roof is least affected by the filling weak surface. This arrangement scheme is the best one, considering the filling weak surface arrangement scheme.

5. Conclusions

Due to the production-process characteristics of the downward slicing and filling method, there are many vertical distribution weak filling surfaces in the roof of the filling body in the mining approach. Indoor tests and numerical simulation were adopted, and numerical simulation comparisons and analysis are carried out on different upper- and lower-layered mining approach layouts considering filling weak surfaces. The following conclusions are drawn.

(1)

The mechanical strength of filling weak surface is much lower than that of surrounding filling body;

(2)

The shear failure process of filling weak surfaces can be divided into two stages: strain-softening stage and strain-hardening stage. With an increase in normal load, the peak shear stress of a filled weak surface rises, and the strain migrates back to the peak shear stress;

(3)

Comparing the numerical simulation results of various schemes, it can be seen that the upper- and lower-layered mining method is the best one when the arrangement is vertical; thus, it is advised that the real engineering mining route arrangement be as similar to the vertical layout as possible. After the mining approach excavation, the mechanical condition of the roof is most similar to the design without covering the weak surface.

This research combined together the direct shear test of weak surfaces of a filling body in the chamber and the stability of filling roofs under different layered mining approaches, which was numerically simulated and analyzed by means of a medium interface unit in FLAC^3D. The optimal layered layout scheme of upper and lower layers is determined, which can provide a theoretical basis for mines adopting the downward horizontal layered cemented filling method.