Abstract
The study on the subsidence of backfill mining block has been a concern of many scholars. A mechanical model of plate subsidence is established by studying the roof of a filling mining area in Luo Iron Mine in this paper. The boundary conditions are given, and the Navier method is used to solve the problem. Based on the thin plate model, the subsidence distribution map of the roof of the underground plate area is obtained. Based on the basic calculation parameters, the influence of a different foundation coefficient, mining depth, length-width ratio of plate area, elastic modulus of roof rock, and thickness of roof on the subsidence of roof is studied. According to the deflection calculation formula obtained, the expression of the internal force and stress is deduced, and the distribution of stress and shear stress on the upper and lower surfaces of the roof is analyzed. The dangerous area of the roof can be obtained, which provides a theoretical basis for the daily maintenance of mine safety.
1. Introduction
Mineral resources are the important material basis of social development and construction. With the development and depletion of shallow mineral resources, more and more minerals are converted to underground mining. Filling mining technology can not only maintain stope stability but also effectively control surface subsidence and protect ecological environment, which has significant safety and environmental benefits [1]. In this paper, the roof of the filling mining panel in the Luoshan iron mine is taken as the research object. According to the theory of elastic-plastic mechanics and the characteristics of underground panel, the mechanical model of roof settlement is established to analyze the settlement law of the underground panel roof, which is expected to provide a theoretical basis for the daily maintenance of mine safety.
The Luohe iron mine is mined in two phases, one in the East and the other in the West. The design scale of phase I is 3 million t/a, which is divided into two stages of -560 m and -620 m. The stage height of -560 m is 120 m, the stage height of -620 m is 60 m, and the first stage is -560 m. The original mining method was sublevel caving without sill pillar. Due to the difficulty of land acquisition and relocation, the mining method is the filling method [2]. According to the characteristics of the underground mining panel, combined with the original solid mechanics theory and various research methods, the paper studies the deformation of overlying strata and surface settlement law in the process of mine filling, which can provide reference for optimizing mining and filling sequence and stope structure parameters, and realize low-cost, efficient, and safe mining [3].
2. Physical Model Description
The buried depth of the Luohe iron ore body is -382 m ~ -864 m. Mining is divided into two stages, first up and then down. Each panel is 126 m long along the strike of the ore body, which is divided into 7 rooms, each room is 18 m, and the mining method of “mining every other one” is adopted. There are 18 m pillars between each panel to support the roof and overburden. The width of each panel is 100 m, and the height is 120 m. The roof rock mass is layered. According to the theory of thin plate, when the ratio of the thickness to the characteristic dimension is about 1/80~1/5, it can be called a thin plate. According to the lithology and strata thickness of 13# geological section in the Luohe iron mine, a single panel and its roof can be regarded as a thin plate model.
The filling body can be regarded as an elastic foundation, and the supporting force is transferred to the basic roof through the direct roof [4]. Therefore, the basic roof can be regarded as a plate with the upper part bearing the load of the overlying strata and the lower part supported by the filling body and the surrounding pillars [5]. The schematic diagram of establishing the settlement mechanics model is shown in Figure 1.
(a) Model stereogram
(b) Model section
3. Mathematical Models
3.1. Basic Equations of Solid Mechanics Model
3.1.1. Equilibrium Differential Equation
When the interior of the object is in a state of equilibrium, it needs to meet the equilibrium conditions [6], and the equilibrium equation is shown in
or abbreviated as
where is the second-order stress tensor inside the body; are the components of the physical force in three coordinates; are the components of the displacement vector at any point in the body in three coordinate directions; and is the density of the body.
3.1.2. Geometric Equation
Because it does not involve the cause of deformation and the material properties of the object, they are generally applicable equations [7]. The geometric equation is as follows:
where is the strain tensor inside the body.
3.1.3. Strain Compatibility Equation (Saint Venant Equation)
In order to make the geometric equations not contradictory, the six strain components must satisfy certain conditions. Therefore, the strain compatibility equation is obtained:
3.1.4. Constitutive Equation
(1) Elastic Stage. There are two forms that can transform each other in this stage [8]: (a)Expression of the strain component by the stress component:(b)Expression of the stress component by the strain component
(2) Plastic Stage. The constitutive equation for elastoplastic materials is as follows:
where is the modulus of elasticity, Pa; is the shear modulus, Pa; is Poisson’s ratio; is the unit tensor; is the lame constant; is the volumetric strain, where ; SIJ is the partial stress tensor; is the bulk modulus; is the proportional factor greater than zero related to the loading history; is the total strain in the plastic stage; and is the total strain in the plastic stage; is called the average stress,
For the plastic region, it also needs to satisfy three equilibrium differential equations and six geometric and constitutive equations. The problem can still be solved by adding an increment and a uniform condition . In the plastic region, the constitutive equation is nonlinear.
3.2. Elastic Thin Plate Theory
Thin plate refers to the plate whose thickness is approximately 1/80~1/5 of the characteristic dimension l. For thin plate bending problems, the following assumptions should be introduced on the basis of continuous [9], uniform, and isotropic assumptions: (1)Any straight line segment perpendicular to the middle plane before deformation remains a straight line after deformation and is perpendicular to the middle plane after bending deformation, and the length remains unchanged. That is, , (2)The normal stress perpendicular to the middle plane of the thin plate is relatively small and neglected, (3)There is no in-plane expansion and shear deformation in the thin plate, i.e., .
3.2.1. Constitutive Relation
(1) Stress-Deformation Relationship. According to the hypothesis, the relationship between the strain component and the midplane deflection can be obtained as shown in equation (8), and the generalized Hooke’s law [9] can be written as the relationship between the principal stress component and the midplane deflection as shown in equation (9):
(2) Internal Force-Deformation Relationship. The relationship between the internal force and deformation can be expressed as
where is the bending rigidity of thin plate, expressed as ; and are bending moment per unit length of the thin plate cross-section; are the torque per unit length of the thin plate cross-section, respectively; and and are the transverse shear force per unit length of the cross-section.
(3) Bending Differential Equation of Thin Plate. After the underground stope is filled, the roof is located on the continuous elastic foundation and is ballasted by the load perpendicular to the slab [10]. When the deflection value is small, according to the Winkler foundation assumption, the reaction force of the filling body at any point under the roof can be expressed as proportional to the deflection of the point, so the load strength of each point of the plate is
where is the coefficient of elastic foundation.
Therefore, the differential equation of the roof bending surface can be expressed as
3.3. Solution of Mechanical Model of Roof Settlement
According to the above modeling analysis, formula (12) is obtained, which is the differential equation of the roof settlement mechanics model. The equation is a high-order partial differential equation, which requires a special solution to obtain the deflection and internal force expressions of the model [11]. Firstly, the boundary conditions of the model are given.
3.3.1. Boundary Conditions
The boundary conditions at and , and of the model are fixed, and the deflection and rotation angle are both 0, so the boundary conditions can be written as
3.3.2. Navier Method for Mechanical Model of Thin Plate
The Navier method is often used to solve the problem of plate deflection in mechanics, that is, double trigonometric series can be used to express the deflection of the thin plate:
where are any positive integers, is the undetermined coefficient, are length and width of the thin plate model, respectively, m.
From the deflection formula, bending moment, and torque formula, the stress components in the thin plate can be obtained as follows:
4. Study on Factors Affecting Roof Settlement
Through the assignment of each parameter, the settlement of different positions of the roof can be calculated [12]. According to the geological data of a mine, the values of basic parameters are shown in Table 1.
The distribution of roof settlement under this parameter is shown in Figure 2.
In Figure 2, the direction represents the advancing direction of the stope, and the direction represents the width of the panel area. It can be seen from the figure that the maximum settlement occurs at m and m. The maximum settlement is 23.3 cm. The roof boundary is limited by boundary conditions, so the displacement is zero. The following is to analyze the influence of different parameters on the roof settlement. It can be seen from the figure that the maximum roof settlement occurs in the middle of the stope. Therefore, based on the basic parameters, the influence of different factors on the maximum roof settlement is studied.
4.1. Foundation Coefficient
The foundation coefficient in the model is used to express the supporting force coefficient of the filling body to the roof after the panel is filled [13]. Different foundation coefficients represent different strength backfills. The variation curve of the maximum settlement of the roof with the foundation coefficient is shown in Figure 3.
The influence degree of different panel length to width ratios on the roof settlement curve is also different. When the ratio of length to width is 3 : 1, the influence of the foundation coefficient on settlement is small. When the length to width ratio is 1 : 1, the maximum settlement of roof is greatly affected by the foundation coefficient. The influence of the foundation coefficient is larger in the initial stage and gradually slows down with the increase of the foundation coefficient.
4.2. Mining Depth
With different mining depths, the load on the roof is different, so the study of different mining depths on the roof settlement is of great significance for mine work. Under different elastic moduli of the roof, the curve of the maximum roof settlement with mining depth is shown in Figure 4.
According to the theory of in situ stress, the load on the roof is equal to the weight of the overburden, so the deeper the mining depth is, the greater the is. It can be seen from the curve in the figure that the maximum roof settlement increases approximately linearly with the increase of depth, because the influence of lateral stress is not considered in the model. The influence of the roof elastic modulus on the maximum settlement of the roof shows that the smaller the elastic modulus is, the greater the growth rate of settlement is.
4.3. Aspect Ratio of Panel
The aspect ratio of the panel is an important parameter for mine design and evaluation of stope safety. In order to study its influence on roof settlement, the variation curves of the maximum roof settlement with panel width under different roof thicknesses are calculated, as shown in Figure 5.
According to the mining design of the mine, the length of a panel is 126 m, so when the stope length is unchanged and the stope width increases from 10 m to 120 m, the variation curve of the maximum roof settlement under different roof thicknesses is obtained. It can be seen from the figure that the maximum roof settlement increases rapidly with the increase of the stope width. When the ratio of stope length to width decreases to 1.5 : 1, the increase slows down. Comparing the curves of different roof thicknesses, it can be seen that the greater the thickness is, the slower the maximum roof settlement increases with the increase of stope width.
4.4. Panel Roof Thickness
The thickness of the roof has an important impact on the safety of the stope. The thickness of the roof depends on the thickness of the direct roof of the panel [1]. Therefore, the parameters of the underground stope can be understood according to the geological data to evaluate and predict the maximum settlement of the roof. In this model, when the roof thickness increases from 4 m to 20 m, the maximum settlement curve of the roof is shown in Figure 6.
It can be seen from the figure that the roof settlement decreases with the increase of the thickness of the roof and decreases rapidly in the initial stage. When the thickness increases to a certain value, the maximum settlement of the roof does not decrease significantly [14]. When the foundation coefficient is 20 MPa/m, it can be seen that the roof thickness has a great influence on the maximum settlement [15]. Therefore, the foundation coefficient is a sensitive parameter when studying the relationship between the maximum settlement and the thickness of the roof.
4.5. Elastic Modulus of Roof Rock
The elastic modulus is an indicator of rock strength [16]. The curve of maximum roof settlement is shown in Figure 7 when the elastic modulus increases from 12 GPa to 60 GPa.
It can be seen from the figure that the larger the elastic modulus is, the smaller the roof settlement is. When the elastic modulus increases linearly, the maximum roof settlement decreases nonlinearly. From the influence of different stope depths on settlement, we can see that the influence of stope depth is linear. In underground mining, due to different geological conditions, the elastic modulus of rock is different, so it is very important to understand the influence of the rock elastic modulus on the maximum roof settlement.
5. Stress Analysis
After analyzing the influencing factors of the roof settlement, the deflection formula can be substituted into the internal force formula to get the stress component and internal force of the roof [17]. According to the assumption of the model, the stress at is 0, so the stress changes at different will be analyzed below.
5.1. Stress at
It can be seen from Table 1 that the thickness of the thin plate is 10 m, so when m, it is the upper surface of the thin plate, and the plate is concave downward, so the stress in the upper part is compressive stress, and the stress in the lower part is tensile stress. According to the calculation, the stress distribution is shown in Figure 8.
(a) Stress in direction
(b) Stress in direction
(c) Shear stress in direction
The stress in and directions and shear stress in direction in the upper surface of the thin plate are, respectively, shown in Figures 8(a) and 8(b) showing that the maximum stress in the thin plate occurs in the middle of the plate, that is, the position of the maximum settlement. The maximum stress in the direction is 60 MPa, while the maximum stress in the direction is 138 MPa, because the model takes 126 m in the direction and 60 m in the direction. It can be seen that the panel parameters have a great influence on the roof stress. According to the boundary conditions, the stress around the thin plate is 0.
Figure 8(c) shows the shear stress in the direction in the thin plate. It can be seen from the figure that the shear stress at the four vertices is the largest, reaching 46 MPa, while the shear stress at the most central position of the plate is 0. The shear stresses at two adjacent vertices are equal in magnitude and opposite in direction. Therefore, it can be seen that the most dangerous part of the roof occurs in the center and four corners of the plate, which is prone to fracture and torsion.
5.2. Stress at
In the last section, the distribution of stress and shear stress on the upper surface of the roof is investigated. The symmetrical position of the surface is at , and the distribution of stress and shear stress is shown in Figure 9.
(a) Stress in direction
(b) Stress in direction
(c) Shear stress in direction
The three figures in Figure 9 show that the distribution of internal stress and shear stress is similar to that in Figure 8, but the directions are opposite. The stress in and directions is -62 MPa and -138 MPa, respectively. It means that the plate bears tensile stress. The magnitude and distribution of shear stress are the same, but the direction is opposite. The damage prone parts are also the center and four corners.
6. Study on Influencing Factors of Overburden Stability
Combined with the mechanical model of underground mining panel roof settlement established above, the numerical simulation method is used to study the influence of different mining depths, panel length to width ratios, roof thicknesses, and roof rock elastic moduli on roof settlement. The applicability of the settlement model can be verified by comparison [5].
6.1. Simulation Results of Influencing Factors of Maximum Roof Settlement
6.1.1. Influence of Mining Depth on Maximum Displacement of Roof
The deeper the mining depth is, the greater the pressure on the panel roof is, and the in situ stress increases correspondingly. Using the simulation method, the settlement curve of the maximum displacement of the panel roof with the calculation step length under the mining depth of 350 m, 450 m, 550 m, and 650 m is simulated, as shown in Figure 10.
6.1.2. Influence of Length to Width Ratio of Panel on Maximum Displacement of Roof
In underground mining, the ratio of length to width determines the exposed area of the stope when the panel length is fixed [18]. The larger the exposed area, the worse the stability of stope. Therefore, in the mathematical model and numerical simulation, the length of the panel is set as 126 m, and the maximum displacement of the panel roof is simulated when the panel width is 40 m, 60 m, 80 m, and 120 m (length to width ratio is 3 : 1, 2 : 1, 1.5 : 1, and 1 : 1, respectively). For each stope width, the results are shown in Table 2.
6.1.3. Influence of Panel Roof Thickness on Maximum Displacement of Roof
The maximum displacement of the roof is simulated when the thickness of the roof is 5 m, 8 m, 10 m, and 12 m, respectively. The results are shown in Table 2.
6.1.4. Influence of Elastic Modulus on Maximum Displacement of Roof
The variation of roof displacement with time step is simulated when the elastic modulus of roof rock is 15 GPa, 30 GPa, 45 GPa, and 60 GPa, respectively. The curve of the maximum displacement point of the roof is shown in Figure 11.
6.2. Comparative Analysis of Mathematical Model Calculation and Numerical Simulation Results
It can be seen from Figure 11 that when the mining depth increases linearly, the increase of the maximum displacement of the panel roof is nonlinear. When the mining depth increases from 550 m to 650 m, the maximum displacement of the roof increases from 0.34 m to 0.46 m, with an increase of 35%. It can be seen that in the process of deep mining, the pressure around the panel increases sharply, and the potential safety hazard also increases. As for the influence of the panel aspect ratio, it can be seen from Table 2 that the influence of the panel aspect ratio and roof thickness is as follows: with the increase of panel width, the exposed area of the panel increases, so the maximum displacement of the roof increases gradually; the greater the thickness of roof, the smaller the settlement.
Through the analysis and comparison of the calculation results of the mathematical model and the numerical simulation results, it can be seen from Table 2 that the influence law of the mining depth and other four factors on the maximum displacement of the panel roof is consistent. In order to verify the accuracy of the mathematical model, the difference between the two results is calculated, and the difference rate is shown in the table. It can be seen that most of the calculated results are within 15%, and only a few of them are 20% ~30%. Generally speaking, the mathematical model can be used to calculate the settlement of the panel roof.
7. Conclusion
(1)Based on the previous research and the theory of elastic-plastic mechanics, this paper analyzes the factors and laws of filling mining affecting overburden settlement. By establishing the mechanical model of panel roof settlement, the boundary conditions are given and solved by the Navier method; the deflection formula and internal force formula of panel roof settlement are obtained. After analysis, the settlement distribution map of the underground panel roof based on the thin plate model is obtained.(2)Based on the basic calculation parameters, the influence of different foundation coefficients, mining depths, panel length to width ratios, roof rock elastic moduli, and roof thicknesses on roof settlement is studied. According to the obtained deflection calculation formula, the expressions of internal force and stress are deduced, and the distribution of stress and shear stress on the upper and lower surface of the roof is analyzed. The dangerous area of the roof can be obtained, which provides a theoretical basis for the daily maintenance of mine safety.
Data Availability
All data used in this study can be obtained by contacting the corresponding author (Gang Huang) (email address: huanggang2016@whut.edu.cn).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the financial support by the National Nature Science Foundation of China under Grant 51804235 and the National Key R&D Plan under Grant Nos. 2018YFC0808405 and 2018YFC0604401.