Stability Analysis of Paste Filling Roof by Cut and Fill Mining

: Ensuring the stability of paste false rooves is an important issue in the study of the process of paste ﬁlling and slicing mining. Here, a mechanical model of a paste false roof is created to analyze its stability in the process of lower slicing mining in order to determine the minimum slicing thickness of the false roof. We use FLAC3D to simulate and analyze the inﬂuence of changes in paste false roof thickness on the stability of the roof. The quantitative functional relationship between the thickness and the subsidence of a false roof, and the optimal thickness of the artiﬁcial paste roof, is ﬁnally obtained by the development law of the plastic zone in the lower slicing face. The results show that when the thickness of the paste false roof is 3.2 m, the roof can maintain its self-stabilization state and ensure the normal mining of lower layers. Because the same thickness of the upper and lower layers is beneﬁcial for mining replacement and equipment selection in different layered working faces, the optimal thickness of a paste false roof is determined to be 3.2 m.


Introduction
Paste filling mining can effectively control surface deformation and overlying rock movement and has become one of the most important technical means to recover coal resources under buildings (structures) [1]. Paste false roofs need to have sufficient strength and a stable bearing capacity. Mastering the instability law of paste filling and layered mining is a prerequisite for safe and efficient mining under pressure [2]. A series of studies have been carried out.
Some researchers established a mechanical model of the basic roof key blocks for paste filling mining and calculated the thickness of the cracked roof [3][4][5]. For example, the mechanical model of the roof rock beam for backfill mining and its differential equations were established. The thin plate theory was used to analyze the stress conditions of the false roof under the mine stope [6,7]. The "embedded beam" mechanical model for the overfilled roof of the roadway was established. The finite-length beam model was used to mechanically analyze the "filling body-direct roof" integral support structure [8][9][10].
At present, the existing research mainly analyzes the stability of the false roof as it relates to the backfill, regarding aspects such as the amount of underfilling, the compression rate of the filling, the filling rate, the mining method, and the overburden load [11,12]. There are few reports on the thickness, strength, and stability of the false roof layer during the layered filling and recovery operation under the paste filling and layered mining [13,14]. Because of the limited space, this study will focus on the influencing factor of the thickness of the paste top layer, and use the method of numerical simulation to study the effect of different thicknesses of the false top layer on its stability [15,16].
Under some mining geological conditions, to effectively control surface deformation, recover the "three-under" compressed coal resources, and digest gangue to protect the ecological environment, paste filling mining is the best solution [17]. At present, the paste filling method has been used to implement paste filling in multiple mines, such as the Taiping Coal Mine of Jining Mining Group, the Daizhuang Coal Mine of Shandong Energy Zi Mining Group, the Xima Coal Mine (Shenyang Province, China)., and the Xiaotun Coal Mine of Jizhong Energy Fengfeng Group. Paste filling mining the full height at one time means extracting the full height of the coal seam at one time and then implementing the filling [18]. The paste filling layer mining method means filling the upper layer first, and the upper layer filling body is used as the false roof for the lower layer. The practice of backfill mining shows that layered mining with paste filling has a better ground deformation control effect than one-time mining with paste filling, and the early strength requirement of the filling is lower than that of the layered mining, so the former has stronger industrial applicability than the latter.
This study takes the paste filling mining project at the E1302 face of the Gaohe Coal Mine of Shanxi Lu'an Group as the background. According to the mining geological conditions of the Gaohe Coal Mine, the paste filling layered coal mining method is selected to recover coal resources. Theoretical analysis combined with numerical simulation is applied to study the relationship between the thickness of the paste false roof and the stability of the false roof.

Project Overview
The E1302 paste filling face is located in the East first panel of the Gaohe Mine, Lu'an, with a length of 230 m and an advanced length of 430 m. The main mining seam is the No. 3 coal seam, and the average coal thickness is 6.4 m. The paste filling is used for two-layer down mining. The coal hardness coefficient is 0.7, the average inclination angle of the coal seam is 8 • , and the roof of the No. 3 coal seam is mudstone, sandy mudstone, siltstone, and partly sandstone. The floor is black mudstone, sandy mudstone, and dark-gray siltstone. The lithological characteristics of the roof and floor rocks are shown in Table 1.
Formula: h 1 -the thickness of the direct roof, m. L 1 -suspended ceiling distance, m. γ 1 -Volume force, MN/m 3 . If the suspended ceiling distance is equal to the controlled ceiling distance, then: (2) Basic top load Q 2 According to the composite beam theory, the high-strength rock layer and the upper, softer rock layer form a composite beam structure, and the load formed by the upper, softer rock layer and the lower key layer can be solved using composite beam theory.
Considering the influence of n layers of overlying rock on the bottom layer, replace (q 1 ) x with (q n ) 1 , and then: Using this formula to calculate the load formed on the bottom rock beam, it must be calculated gradually from bottom to top, that is, the sequence is (q 1 ) 1 , (q 2 ) 1 , ..., (q n ) 1 ; when the calculation is (q i + 1 ) 1 < (q i ), then the calculation stops at one o'clock. It shows that the i + 1 layer does not generate additional load on the bottom layer. At the time (q n ) 1 = (q i ) 1 , it can be judged that the n + 1th layer does not affect the load of the first layer, that is, the strength of the nth layer on top of the rock. It is higher and can support the overburdened load without load on the basic roof. After calculating (q 6 ) 1 < (q 5 ) 1 , the sixth layer of siltstone above the basic roof has the characteristics of large thickness and high strength, and the overlying load will not be transmitted downwards. Therefore, the load on the basic roof is Q 2 = (q 5 ) 1 = 0.443 MPa.
(3) Calculation of self-load of paste false top The load of the paste false top is considered to be safe, and it is calculated according to the larger thickness of 5.2 m: In summary, the false roof bearing load q is calculated as:

Theoretical Analysis of False Roof Instability Thickness
There will be no collapse zone on the roof above the mining face with the paste filling, but only a fracture zone and curved subsidence zone [19]. The practice has shown that the roof control effect of paste filling mining is generally affected by factors such as the amount of roof and floor moving before filling, the amount of underfilling, the compression rate, and the filling rate of the filling body. With the continuous improvement of the filling process and the development of key filling equipment, paste filling technology has become more mature. The amount of top and bottom plate moving before filling, the amount of underfilling, and the compression rate of the filling body can be controlled within a small range, while still meeting the requirements [20,21]. The filling body fills the mined-out area and controls the roof in time, and therefore effectively controls the roof subsidence and the filling requirements of the surface collapse. At present, the commonly used mechanical models of the instability mechanism of the paste false top are the thin "plate" theoretical model and the simply supported "beam" theoretical model [22]. Both theoretical models assume that the paste false top is continuous and homogeneous, and that it was a linear elastic body before yield failure. Since the length of the filling working face in the inclined direction is much longer than the span in the strike direction, the roof in the middle of the working face can be regarded as a beam of countless unit widths inserted into the front and rear rock masses, so it will be filled. The artificial top is considered to be a simply supported beam for research purposes [22]. The load of the simply supported beam and above is transmitted from the beam to the front and rear of the stope before it breaks. At this time, the overburden load is regarded as the uniform load q, the false roof thickness h, and the rock beam length L, as shown in Figure 1. The mechanical model analyzes its force [23]. From elastic mechanics, we know that the bending stress in the beam σ x is caused by the bending moment, the extrusion stress σ y is caused by the overlying load q, and τ xy is caused by the shear stress. Assuming that σ y is only a function of y, we have: The Airy stress function is: Bringing the formula into the Airy stress function Formula (7), we obtain: Solve the stress function from the compatibility equation, put the Equation (8) into the compatibility equation, and simplify it to obtain: Omit the primary term and constant term, and further simplify to obtain the stress function: Φ = x 2 2 Ay 3 + By 2 + Cy + D + x Ey 3 + Fy 2 + Gy − A 10 y 5 − B 6 y 4 + Hy 3 + Ky 2 Solve the stress component from the stress function, and put the Formula (10) into the Airy stress function (7) to obtain the stress component: Considering that the span of the beam is much larger than the thickness of the beam, the stress boundary conditions on the main boundary can be satisfied, and the integral stress boundary of the Saint-Venant principle is used to replace the minor part of the boundary so that the boundary conditions are approximately satisfied: Therefore, first consider the main boundary conditions, combine the formula with the formula, and solve: Next, consider the left and right boundary conditions. Since the model is arranged symmetrically, one side can be considered. On the right side of the beam, there is no horizontal plane force, that is, σ x = 0, so: Combining Equations (15)- (17) with conditional Equations (18) and (19), the stress component is solved as: It can be seen from the stress expressions (20)-(22) that in a long beam whose length is far greater than its thickness, the bending stress σ x is of the same order as q(L/h) 2 , which is the main stress; the shear stress τ xy is the same as q(L/h) 2 . The same order of q(L/h) is the secondary stress; the extrusion stress σ y is the same order of q, and this is secondary stress, which is not considered. Therefore, the failure of the false roof is mainly caused by the bending moment and shearing force.
Breaking of the false roof is caused by the bending moment: The maximum bending stress on the beam reaches its maximum value at y = h/2, so: Omit the high-order infinitesimal terms, then: If the tensile strength of the false top is F T , the critical conditions for the bending instability and breaking of the paste false top are considered as follows: Then there are: Breaking of the false roof is caused by shearing force: If the maximum bending stress on the beam reaches its maximum value on the neutral axis at y = 0, omitting the high-order infinitesimal terms, then: If the tensile strength of the false roof is F S , the critical condition of shear instability and fracture of the paste false roof is considered as follows: Then there are: According to the formula, the minimum thickness of false roof h min is obtained because it also satisfies: This is the case considering that x = 0.5 L is the dangerous section of the false roof's simply supported beam. According to the paste material mechanics experiment, the 28 d tensile strength of the paste false roof is 0.5 MPa, and the shear strength is 1.5 MPa. The false roof span L is 8 m, and the false roof bears the load of 0.754 MPa, and Formulas (30) and (31) are used to obtain: Therefore, the comprehensive analysis shows that the thickness of the false roof layer should be greater than 3.01 m.

Numerical Simulation Modeling of Layered Mining
The length of the E1302 working face is 230 m, the advance length is 430 m, the average thickness is 6.4 m, and the buried depth of the three coal seams is about 420 m. According to the actual filling and mining process of the working face, the geometric model of the stope structure of paste filling and layered mining is established as shown in Figure 2. By referring to the geometric model and combining the physical and mechanical parameters of the rock formation in Table 2, we can use FLAC 3D [24,25] to establish a numerical simulation model as shown in Figure 3. Horizontal displacement constraints are imposed on the north and south boundaries of the model, and the bottom boundary of the model is fixed in the vertical direction. The top rock layer is buried at about 310 m. Therefore, a vertical downward uniform load of 7.75 MPa is applied to the top to simulate the overburden load. To eliminate the influence of the boundary of the model, a boundary coal pillar of 60 m is left on both sides of the working face along the advancing direction. The size of the model is 550 m × 300 m × 160 m in the order of length × width × height. In the process of analysis and calculation, the weight of each rock layer is considered in the vertical direction of the model, and the Mohr-Coulomb strength criterion is used to solve the constitutive relationship [26,27].

Numerical Simulation Schemes of Different False Roof Thicknesses
According to the actual mining and filling operation of the E1302 working face of Gaohe filling mining, the span of the roof control area of the working face is about 8 m. The upper stratified filling body is considered by the unidirectional force state, and the later strength requirement of the filling body is determined to be 5 MPa according to the Bieniawski paste late strength formula [28,29]. In the simulation, the controlled variable method is adopted to control the layer thickness as a single variable. According to the actual mining and filling operation, the compressive strength of the filling body false top is set to 5 MPa, and the collapse distance of the lower layer working face is 8 m. Other simulation variables are assumed to be consistent. The sum of the height of the false top layer and the lower layer working face of the model is the total thickness of the coal seam, which is 6.4 m. Based on the experience in the field, by changing the thickness of the false roof (2.2 m, 3.2 m, 4.2 m, 5.2 m), the distribution of the false roof plastic zone and the amount of roof subsidence have been analyzed. Previously, scholars have used the method of finite element numerical simulation to analyze the relationship between the displacement of the middle part of the false roof section and the instability and fracture of the false roof under the action of the weight and the upper rock formation according to the simulation results of the displacement field [10]. Some scholars also use FLAC to establish a numerical simulation model and analyze the critical conditions for roof failure and the instability of the working face under different filling steps by observing the plastic zone.
In this paper, the FLAC3D 5.0 numerical simulation method is used to study the roof of the working face through the development of the plastic zone, and the influence of the false roof thickness on the stability of the false roof is studied [30]. In the 3D geometric model, the blue survey lines and pink survey points are arranged at the bottom of the upper paste filling roof, and 19 survey points are set at equal intervals with an interval of 1 m. The survey line layout section is shown in Figure 4.  According to the different thickness conditions in the scheme, a simulation is carried out, and the displacement cloud diagram of the lower layer working face is shown in Figure 5. It can be seen from Figure 5 that: when the thickness of the false roof is 2.2 m, the maximum sinking amount of the roof is 81 mm; when the thickness of the false roof is 3.2 m, the maximum sinking amount of the roof is 75 mm; and when the thickness of the false roof is 4.2 m, the sinking amount of the roof is at its maximum. When the false roof is 72 mm and the thickness of the false roof is 5.2 m, the maximum subsidence of the roof is 58 mm. It can be roughly judged that the maximum subsidence of the false roof is negatively correlated with the thickness of the false roof.
One can draw the false roof subsidence curve of each measurement point according to the cloud map displacement. As shown in Figure 6, the false roof subsidence area can be divided into three areas, A, B, and C, according to the distribution law of the subsidence curve; that is, the false roof subsidence boundary zone (zone A), the false roof subsidence transition zone (zone B), and the false roof subsidence central zone (zone C). Select different sections in each area of A, B, and C in Figure 6, and one can draw the false roof in sections. The relationship between the amount of subsidence and the thickness of the false roof is shown in Figure 7.
To further quantitatively analyze the functional relationship between the false roof thickness and the false roof subsidence in the central area of false roof subsidence (zone C), the distance between the working face and the origin of the survey line is x, and when 52 m ≤ x ≤ 57 m, C is obtained the area function relationship is Formula (32), the degree of fit R2 = 0.989, and the fitted curve is shown in Figure 8.
From Figure 8 and Formula (32), it can be seen that in the central area of false roof subsidence (zone C), the thickness of the false roof is negatively correlated with the amount of false roof subsidence, and the amount of subsidence decreases with the increase of false roof thickness. The slope of the curve first gradually decreases and then increases, which means that when the false roof thickness is between 2.2 and 3.2 m, the subsidence speed of the false roof decreases with the increase in thickness. When the false roof thickness is between 3.2 and 4.2 m, the subsidence value of the false roof changes little. The sinking speed of the false roof increases with the increase of thickness after 4.2 m.
Similarly, quantitatively analyze areas B and C, and one can obtain the piecewise function of the correlation between the false roof thickness and the false roof subsidence:

Analysis of the Development Law of Plastic Zone under Different False Top Thickness
The plastic zone distribution of the lower layer working face with different filling body false roof thickness is shown in Figure 9. It can be seen from the figure that the failure rule of the plastic zone around the working face is the same. The false roof subsidence center area mainly has a tensile failure, and the subsidence boundary area has a shear failure. In the lower part of the subsidence transition area, the working face bottom plate is subjected to shear stress and tension at the same time. Plastic failure occurs under the influence of stress, and it is in a state of multi-directional stress failure.
When the thickness of the false roof is 2.2 m, tensile failure occurs mainly in the sinking center area of the false roof, while shear failure occurs in the sinking boundary area. When the false roof thickness is 3.2 m, there is a partial tensile failure within 1 m of the subsidence center area, and shear failure occurs in the subsidence boundary area. When the thickness of the false roof is 3.2~4.2 m, the deep part of the plastic zone does not increase significantly, but rather shear failure occurs in the sinking boundary area, and the plastic failure range does not increase significantly with the increase of the false roof thickness in the vertical and horizontal directions. When the thickness of the false roof is 4.2~5.2 m, the distribution area of the plastic zone around the working face is significantly reduced, and only less shear deformation occurs in the subsidence boundary area, and the false roof is in a stable state.
In summary, when the thickness of the false roof is less than 3.2 m, the extent of the plastic failure zone is significantly reduced as the thickness of the false roof increases; when the thickness of the false roof is 3.2-4.2 m, the thickness of the false roof increases, and in the plastic failure zone there is no significant expansion of the range; when the false roof thickness is greater than 4.2 m, the plastic failure zone will be further reduced with the increase of the false roof thickness. It should be pointed out that though the study is based on theoretical and numerical investigation, the results are relatively limited due to the selection of parameters. In the future, a field study will be conducted and then compared to the numerical study, which will be a guideline in this study.

Conclusions
(1) Based on the geological conditions of the Gaohe Mine, a mechanical model of the paste false top was established. By analyzing the two limit states of the paste false top in tensile failure and shear failure, it was concluded that the thickness of the upper layer should be at least 3.01 m for stability.
(2) The false roof subsidence area is divided into three areas: false roof subsidence boundary area, false roof subsidence transition area, and false roof subsidence central area, and the false roof thickness in the central area of false roof subsidence is obtained by the fitting curve. The cubic function relationship of the subsidence and the corresponding change trend of the false roof subsidence with the increase of thickness is analyzed through the slope of the curve.
(3) By analyzing the distribution of the plastic zone in the lower layered working face, it can be seen that when the false roof thickness is 2.2-3.2 m, the plastic failure range decreases with the increase of the false roof thickness; when the false roof thickness is 3.2-4.2 m, the plastic failure zone does not expand significantly with the increase of the false roof thickness; when the false roof thickness is 4.2-5.2 m, the plastic failure zone is further reduced with the increase of the false roof thickness. Since the total thickness of the top layer and the bottom layer is 6.4 m, the thickness of the top layer is too large or the thickness of the top layer is too small, which will increase the difficulty of mining, increase the requirements for support and equipment selection, and is not conducive to the work of upper and lower layers. Therefore, it is finally determined that the thickness of the false top layer of the paste filling is 3.2 m.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.