# Study on the Settlement Mechanism with Effect of Backfilling Method of an Underground Mine

Authors: Dequan Xuan; Zhiming Liu; Gang Huang; Jianhua Zhang; Tingting Jiang; Bo Ke
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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 thicknessh to the characteristic dimension L 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.Figure 1
Mechanical model of roof settlement.
(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
(1)∂σx∂x+∂τxy∂y+∂τxz∂z+X=0orρ∂2u∂t2,∂τyx∂x+∂σy∂y+∂τyz∂z+Y=0orρ∂2v∂t2,∂τzx∂x+∂τzy∂y+∂σz∂z+Z=0orρ∂2w∂t2,or abbreviated as(2)σij,j+Xi=0orρ∂2u∂t2,whereσij is the second-order stress tensor inside the body; X,Y,andZ are the components of the physical force in three coordinates; u,v,andw 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:
(3)εij=12ui,j+uj,i,whereεij 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:(4)εij,kl+εkl,ij=εik,jl+εjl,ik.


#### 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:(5)εij=1+υEσij−υEσkkδij(b)
Expression of the stress component by the strain component(6)σij=λεkkδij+2μεij=λΘδij+2μεij(2) Plastic Stage. The constitutive equation for elastoplastic materials is as follows:
(7)dεij=12GdSij+dλSij+1−2υ3EdσkkδijwhereE is the modulus of elasticity, Pa; G is the shear modulus, Pa; υ is Poisson’s ratio; δij is the unit tensor; μ,λ is the lame constant; Θ is the volumetric strain, where Θ=εx+εy+εZ; SIJ is the partial stress tensor; K is the bulk modulus; dλ is the proportional factor greater than zero related to the loading history; eij is the total strain in the plastic stage; and Θ is the total strain in the plastic stage; σm is called the average stress, σm=σx+σy+σz/3For 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 incrementdλ and a uniform condition fσij=0. In the plastic region, the constitutive equation is nonlinear.


### 3.2. Elastic Thin Plate Theory

Thin plate refers to the plate whose thicknessh 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,γzx=γxz=0, εz=0(2)
The normal stress perpendicular to the middle plane of the thin plate is relatively small and neglected,σz=0(3)
There is no in-plane expansion and shear deformation in the thin plate, i.e.,ux,y,0=vx,y,0=0.


#### 3.2.1. Constitutive Relation

(1) Stress-Deformation Relationship. According to the hypothesis, the relationship between the strain component and the midplane deflection w 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 w as shown in equation (9):
(8)εx=−z∂2w∂x2,εy=−z∂2w∂y2,εxy=−2z∂2w∂x∂y,(9)σx=−Ez1−υ2∂2w∂x2+∂2w∂y2,σy=−Ez1−υ2∂2w∂y2+∂2w∂x2,τxy=−Ez1+υ∂2w∂x∂y.(2) Internal Force-Deformation Relationship. The relationship between the internal force and deformation can be expressed as
(10)Mx=−D∂2w∂x2+υ∂2w∂y2,My=−D∂2w∂y2+υ∂2w∂x2,Mxy=Myx=−D1−υ∂2w∂x∂y,Qx=∫h/2h/2τxzdz=−D∂∂x∇2w,Qy=∫h/2h/2τyzdz=−D∂∂y∇2w,whereD is the bending rigidity of thin plate, expressed as D=Eh3/121−υ; Mx and My are bending moment per unit length of the thin plate cross-section; Mxy are Myx the torque per unit length of the thin plate cross-section, respectively; and Qx and Qy 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 qx,y 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
(11)Qx,y=qx,y−k⋅wx,y,wherek is the coefficient of elastic foundation.Therefore, the differential equation of the roof bending surface can be expressed as(12)∂4w∂x4+2∂4w∂x2∂y2+∂4w∂y4=qx,y−k⋅wx,yD.


### 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 atx=0 and x=a, y=0 and y=b of the model are fixed, and the deflection and rotation angle are both 0, so the boundary conditions can be written as
(13)wx=0=wx=a=0,∂w∂xx=0=∂w∂xx=a=0,wy=0=wy=b=0,∂w∂yy=0=∂w∂yy=b=0.


#### 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:(14)wx,y=∑m=1∞∑n=1∞Amnsinmπxasinnπxb,wheremandn are any positive integers, Amn is the undetermined coefficient, aandb 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:(15)σxx,y,z=12Mxh3z,σyx,y,z=12Myh3z,σxx,y,z=12Mxyh3z.


## 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
(1)∂σx∂x+∂τxy∂y+∂τxz∂z+X=0orρ∂2u∂t2,∂τyx∂x+∂σy∂y+∂τyz∂z+Y=0orρ∂2v∂t2,∂τzx∂x+∂τzy∂y+∂σz∂z+Z=0orρ∂2w∂t2,or abbreviated as(2)σij,j+Xi=0orρ∂2u∂t2,whereσij is the second-order stress tensor inside the body; X,Y,andZ are the components of the physical force in three coordinates; u,v,andw 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:
(3)εij=12ui,j+uj,i,whereεij 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:(4)εij,kl+εkl,ij=εik,jl+εjl,ik.


### 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:(5)εij=1+υEσij−υEσkkδij(b)
Expression of the stress component by the strain component(6)σij=λεkkδij+2μεij=λΘδij+2μεij(2) Plastic Stage. The constitutive equation for elastoplastic materials is as follows:
(7)dεij=12GdSij+dλSij+1−2υ3EdσkkδijwhereE is the modulus of elasticity, Pa; G is the shear modulus, Pa; υ is Poisson’s ratio; δij is the unit tensor; μ,λ is the lame constant; Θ is the volumetric strain, where Θ=εx+εy+εZ; SIJ is the partial stress tensor; K is the bulk modulus; dλ is the proportional factor greater than zero related to the loading history; eij is the total strain in the plastic stage; and Θ is the total strain in the plastic stage; σm is called the average stress, σm=σx+σy+σz/3For 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 incrementdλ and a uniform condition fσij=0. In the plastic region, the constitutive equation is nonlinear.


## 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
(1)∂σx∂x+∂τxy∂y+∂τxz∂z+X=0orρ∂2u∂t2,∂τyx∂x+∂σy∂y+∂τyz∂z+Y=0orρ∂2v∂t2,∂τzx∂x+∂τzy∂y+∂σz∂z+Z=0orρ∂2w∂t2,or abbreviated as(2)σij,j+Xi=0orρ∂2u∂t2,whereσij is the second-order stress tensor inside the body; X,Y,andZ are the components of the physical force in three coordinates; u,v,andw 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:
(3)εij=12ui,j+uj,i,whereεij 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:(4)εij,kl+εkl,ij=εik,jl+εjl,ik.


## 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:(5)εij=1+υEσij−υEσkkδij(b)
Expression of the stress component by the strain component(6)σij=λεkkδij+2μεij=λΘδij+2μεij(2) Plastic Stage. The constitutive equation for elastoplastic materials is as follows:
(7)dεij=12GdSij+dλSij+1−2υ3EdσkkδijwhereE is the modulus of elasticity, Pa; G is the shear modulus, Pa; υ is Poisson’s ratio; δij is the unit tensor; μ,λ is the lame constant; Θ is the volumetric strain, where Θ=εx+εy+εZ; SIJ is the partial stress tensor; K is the bulk modulus; dλ is the proportional factor greater than zero related to the loading history; eij is the total strain in the plastic stage; and Θ is the total strain in the plastic stage; σm is called the average stress, σm=σx+σy+σz/3For 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 incrementdλ and a uniform condition fσij=0. In the plastic region, the constitutive equation is nonlinear.


## 3.2. Elastic Thin Plate Theory

Thin plate refers to the plate whose thicknessh 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,γzx=γxz=0, εz=0(2)
The normal stress perpendicular to the middle plane of the thin plate is relatively small and neglected,σz=0(3)
There is no in-plane expansion and shear deformation in the thin plate, i.e.,ux,y,0=vx,y,0=0.


### 3.2.1. Constitutive Relation

(1) Stress-Deformation Relationship. According to the hypothesis, the relationship between the strain component and the midplane deflection w 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 w as shown in equation (9):
(8)εx=−z∂2w∂x2,εy=−z∂2w∂y2,εxy=−2z∂2w∂x∂y,(9)σx=−Ez1−υ2∂2w∂x2+∂2w∂y2,σy=−Ez1−υ2∂2w∂y2+∂2w∂x2,τxy=−Ez1+υ∂2w∂x∂y.(2) Internal Force-Deformation Relationship. The relationship between the internal force and deformation can be expressed as
(10)Mx=−D∂2w∂x2+υ∂2w∂y2,My=−D∂2w∂y2+υ∂2w∂x2,Mxy=Myx=−D1−υ∂2w∂x∂y,Qx=∫h/2h/2τxzdz=−D∂∂x∇2w,Qy=∫h/2h/2τyzdz=−D∂∂y∇2w,whereD is the bending rigidity of thin plate, expressed as D=Eh3/121−υ; Mx and My are bending moment per unit length of the thin plate cross-section; Mxy are Myx the torque per unit length of the thin plate cross-section, respectively; and Qx and Qy 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 qx,y 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
(11)Qx,y=qx,y−k⋅wx,y,wherek is the coefficient of elastic foundation.Therefore, the differential equation of the roof bending surface can be expressed as(12)∂4w∂x4+2∂4w∂x2∂y2+∂4w∂y4=qx,y−k⋅wx,yD.


## 3.2.1. Constitutive Relation

(1) Stress-Deformation Relationship. According to the hypothesis, the relationship between the strain component and the midplane deflection w 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 w as shown in equation (9):
(8)εx=−z∂2w∂x2,εy=−z∂2w∂y2,εxy=−2z∂2w∂x∂y,(9)σx=−Ez1−υ2∂2w∂x2+∂2w∂y2,σy=−Ez1−υ2∂2w∂y2+∂2w∂x2,τxy=−Ez1+υ∂2w∂x∂y.(2) Internal Force-Deformation Relationship. The relationship between the internal force and deformation can be expressed as
(10)Mx=−D∂2w∂x2+υ∂2w∂y2,My=−D∂2w∂y2+υ∂2w∂x2,Mxy=Myx=−D1−υ∂2w∂x∂y,Qx=∫h/2h/2τxzdz=−D∂∂x∇2w,Qy=∫h/2h/2τyzdz=−D∂∂y∇2w,whereD is the bending rigidity of thin plate, expressed as D=Eh3/121−υ; Mx and My are bending moment per unit length of the thin plate cross-section; Mxy are Myx the torque per unit length of the thin plate cross-section, respectively; and Qx and Qy 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 qx,y 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
(11)Qx,y=qx,y−k⋅wx,y,wherek is the coefficient of elastic foundation.Therefore, the differential equation of the roof bending surface can be expressed as(12)∂4w∂x4+2∂4w∂x2∂y2+∂4w∂y4=qx,y−k⋅wx,yD.


## 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 atx=0 and x=a, y=0 and y=b of the model are fixed, and the deflection and rotation angle are both 0, so the boundary conditions can be written as
(13)wx=0=wx=a=0,∂w∂xx=0=∂w∂xx=a=0,wy=0=wy=b=0,∂w∂yy=0=∂w∂yy=b=0.


### 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:(14)wx,y=∑m=1∞∑n=1∞Amnsinmπxasinnπxb,wheremandn are any positive integers, Amn is the undetermined coefficient, aandb 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:(15)σxx,y,z=12Mxh3z,σyx,y,z=12Myh3z,σxx,y,z=12Mxyh3z.


## 3.3.1. Boundary Conditions

The boundary conditions atx=0 and x=a, y=0 and y=b of the model are fixed, and the deflection and rotation angle are both 0, so the boundary conditions can be written as
(13)wx=0=wx=a=0,∂w∂xx=0=∂w∂xx=a=0,wy=0=wy=b=0,∂w∂yy=0=∂w∂yy=b=0.


## 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:(14)wx,y=∑m=1∞∑n=1∞Amnsinmπxasinnπxb,wheremandn are any positive integers, Amn is the undetermined coefficient, aandb 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:(15)σxx,y,z=12Mxh3z,σyx,y,z=12Myh3z,σxx,y,z=12Mxyh3z.


## 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.Table 1
Calculation parameter value.
NumberParameter name (unit)Parameter value1Average density of rock (kg/m3)28002Elastic modulus (GPa)303Foundation coefficient (MPa/m)404Panel length (m)1265Panel width (m)606Roof thickness (m)107Buried depth (m)4508Poisson’s ratio0.25The distribution of roof settlement under this parameter is shown in Figure2.Figure 2
Distribution of roof settlement under basic parameters.In Figure2, the X direction represents the advancing direction of the stope, and the Y direction represents the width of the panel area. It can be seen from the figure that the maximum settlement occurs at x=63 m and y=30 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.Figure 3
Curve of the maximum settlement of roof with the variation of foundation coefficient.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 Figure4.Figure 4
Influence of different mining depths on roof settlement.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 theqx,y 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 Figure5.Figure 5
Influence of different aspect ratios of the panel on roof settlement.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.Figure 6
Influence of different roof thicknesses on roof settlement.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.Figure 7
Influence of different elastic moduli of the roof rock on roof settlement.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.


## 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.Figure 3
Curve of the maximum settlement of roof with the variation of foundation coefficient.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 Figure4.Figure 4
Influence of different mining depths on roof settlement.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 theqx,y 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 Figure5.Figure 5
Influence of different aspect ratios of the panel on roof settlement.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.Figure 6
Influence of different roof thicknesses on roof settlement.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.Figure 7
Influence of different elastic moduli of the roof rock on roof settlement.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 z=0 is 0, so the stress changes at different z will be analyzed below.


### 5.1. Stress atz=h/2

It can be seen from Table1 that the thickness of the thin plate is 10 m, so when z=5 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.Figure 8
Stress diagram atz=h/2.
(a)
Stress inx direction(b)
Stress iny direction(c)
Shear stress inxy directionThe stress inx and y directions and shear stress inxy 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 x direction is 60 MPa, while the maximum stress in the y direction is 138 MPa, because the model takes 126 m in the x direction and 60 m in the y 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.Figure8(c) shows the shear stress in the xy 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 atz=−h/2

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 atz=−h/2, and the distribution of stress and shear stress is shown in Figure 9.Figure 9
Stress diagram atz=−h/2.
(a)
Stress inx direction(b)
Stress iny direction(c)
Shear stress inxy directionThe three figures in Figure9 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 x and y 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.


## 5.1. Stress atz=h/2

It can be seen from Table1 that the thickness of the thin plate is 10 m, so when z=5 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.Figure 8
Stress diagram atz=h/2.
(a)
Stress inx direction(b)
Stress iny direction(c)
Shear stress inxy directionThe stress inx and y directions and shear stress inxy 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 x direction is 60 MPa, while the maximum stress in the y direction is 138 MPa, because the model takes 126 m in the x direction and 60 m in the y 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.Figure8(c) shows the shear stress in the xy 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 atz=−h/2

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 atz=−h/2, and the distribution of stress and shear stress is shown in Figure 9.Figure 9
Stress diagram atz=−h/2.
(a)
Stress inx direction(b)
Stress iny direction(c)
Shear stress inxy directionThe three figures in Figure9 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 x and y 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 Figure10.Figure 10
Maximum roof settlement displacement curve under different mining depths.


#### 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.Table 2
Data comparison between the mathematical model calculation and numerical simulation results.
Influence factorParameterCalculated value of the model (m)Numerical simulation value (m)Difference rate (%)Mining depth350 m0.1800.1884.4450 m0.2330.26112.0550 m0.2850.33417.2660 m0.3370.45635.3Aspect ratio3 : 10.1000.12727.02 : 10.2330.26112.11.5 : 10.3100.289-6.81 : 10.3700.317-14.3Roof thickness5 m0.3780.296-21.78 m0.2960.280-5.410 m0.2330.26112.012 m0.1770.24035.6Elastic modulus15 GPa0.2980.3052.330 GPa0.2330.26112.045 GPa0.1920.1962.160 GPa0.1600.145-9.4Note:differencerate=simulationvalue−calculatedvalue/calculatedvalue∗100%.


#### 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 Table2.


#### 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 Figure11.Figure 11
Maximum settlement displacement curve of roof under different elastic moduli of rock.


### 6.2. Comparative Analysis of Mathematical Model Calculation and Numerical Simulation Results

It can be seen from Figure11 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 Table2 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.


## 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 Figure10.Figure 10
Maximum roof settlement displacement curve under different mining depths.


### 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.Table 2
Data comparison between the mathematical model calculation and numerical simulation results.
Influence factorParameterCalculated value of the model (m)Numerical simulation value (m)Difference rate (%)Mining depth350 m0.1800.1884.4450 m0.2330.26112.0550 m0.2850.33417.2660 m0.3370.45635.3Aspect ratio3 : 10.1000.12727.02 : 10.2330.26112.11.5 : 10.3100.289-6.81 : 10.3700.317-14.3Roof thickness5 m0.3780.296-21.78 m0.2960.280-5.410 m0.2330.26112.012 m0.1770.24035.6Elastic modulus15 GPa0.2980.3052.330 GPa0.2330.26112.045 GPa0.1920.1962.160 GPa0.1600.145-9.4Note:differencerate=simulationvalue−calculatedvalue/calculatedvalue∗100%.


### 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 Table2.


### 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 Figure11.Figure 11
Maximum settlement displacement curve of roof under different elastic moduli of rock.


## 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 Figure10.Figure 10
Maximum roof settlement displacement curve under different mining depths.


## 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.Table 2
Data comparison between the mathematical model calculation and numerical simulation results.
Influence factorParameterCalculated value of the model (m)Numerical simulation value (m)Difference rate (%)Mining depth350 m0.1800.1884.4450 m0.2330.26112.0550 m0.2850.33417.2660 m0.3370.45635.3Aspect ratio3 : 10.1000.12727.02 : 10.2330.26112.11.5 : 10.3100.289-6.81 : 10.3700.317-14.3Roof thickness5 m0.3780.296-21.78 m0.2960.280-5.410 m0.2330.26112.012 m0.1770.24035.6Elastic modulus15 GPa0.2980.3052.330 GPa0.2330.26112.045 GPa0.1920.1962.160 GPa0.1600.145-9.4Note:differencerate=simulationvalue−calculatedvalue/calculatedvalue∗100%.


## 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 Table2.


## 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 Figure11.Figure 11
Maximum settlement displacement curve of roof under different elastic moduli of rock.


## 6.2. Comparative Analysis of Mathematical Model Calculation and Numerical Simulation Results

It can be seen from Figure11 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 Table2 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.



---

*Source: 5929116-2021-08-06.xml*