Theoretical Analysis on Stress and Deformation of Overburden Key Stratum in Solid Filling Coal Mining Based on the Multilayer Winkler Foundation Beam Model

Solid backfill coal mining (SBCM) is a green mining technology which can effectively alleviate the environmental problems induced by traditional coal mining techniques, such as surface subsidence, water resources loss, coal gangue occupation, and pollution. In this study, a multilayer Winkler foundation beam model for the overburden key strata is proposed, and the model with two key strata is solved. The subsidence, rotating angle, inner force, and stress of the overburden key strata are systematically analyzed under various backfill elastic modulus, mining height, and soft layer thickness. The results show that the subsidence of the key strata exhibit “basin”-shape curves, and the backfill elastic modulus, mining height, and the thickness of the soft strata have significant influences on the subsidence of the key strata. The shear stress, horizontal stress, and vertical stress of key stratum can be effectively reduced by increasing the backfill elastic modulus. The increase of mining height has little influence on the stress of key stratum that close to the coal seam (key stratum #1), but has a significant effect on the stress of key stratum that above the soft layers (key stratum #2). On the contrary, the effect of increasing soft layer thickness on the stress of key stratum is opposite to that of increasing mining height. In addition, the shear failure of key stratum #1 at mining boundary and the tensile failures on both sides of mining boundary should be preferentially considered in SBCM engineering design. Due to the low shear stress level of key stratum #2, the tensile failure on both sides of the mining boundary should be mainly considered.


Introduction
Solid backfill coal mining (SBCM) is a coal mining technology in which underground coal resources are replaced by backfill [1,2]. The backfill normally consists of industrial wastes, such as coal gangue and fly ash [3][4][5]. In this respect, SBCM can effectively alleviate the environmental problems caused by gangue discharge in coal mining, including occupying the ground and spontaneous combustion [6][7][8][9]. Besides, SBCM can improve the mining rate efficiently, especially for mining the coal resources under the surface water bodies, buildings, and railways (roads) [1]. In addition, compared with the traditional underground coal mining technol-ogy, the backfill can support the overburden strata in SBCM, which can effectively reduce a series of problems caused by overburden strata deformation and failure, such as surface subsidence [10,11] and water resources loss [12][13][14][15][16]. Therefore, it is of great significance to study the deformation and stress distribution of overburden strata for SBCM engineering.
The research on SBCM mainly focuses on the backfill material and the overburden control. At present, the mechanical properties, components, and particle sizes of backfill materials have relatively systematic research results [17][18][19][20][21]. However, to the best of our knowledge, there are few theoretical and numerical modeling studies on the overburden strata control.
For the numerical simulation studies, Zha et al. studied the characteristics of overburden deformation and movement by the numerical simulation [22]; Huang et al. analyzed the effect of compaction rates of backfill on controlling movement of overburden and surface subsidence by UDEC [23]; Li et al. proposed a backfill material model that considers the coefficient of horizontal pressure and numerically simulated the overburden movement deformation characteristics by Flac 3D [24]. However, the above works only focus on the characteristics of overburden movement and deformation without considering the stress distribution and failure of the overburden key strata. Zhang et al. proposed a negative exponential function model of backfill materials and simulated the deformation and stress distribution of overburden of SBCM by the software of Abaqus [25]. However, this work did not pay attention to the stress distribution of overburden key strata.
For the theory studies, Miao et al. proposed an equivalent mining height theory, which assumes that the actual mining height of SBCM can be equivalent to the deformation of backfill [2]. Based on this theory, the Winkler foundation beam model of the upper roof was established, and the analytical solution of the subsidence and internal force of the old roof were obtained. Similar to the work of Miao et al., Chen et al. proposed a beam model for elastic foundation of roof and systematically analyzed the influences of mining depth and backfill foundation coefficient on roof subsidence [26]. In addition, Li et al. simplified the basic roof into a thin plate, proposed a thin plate model for the overburden key strata of SBCM, and systematically analyzed the influence of compaction rate of backfill on deformation of overburden key strata [27]. The above theoretical models only consider the key stratum close to the coal seam, and the load on the key stratum is simplified to uniformly distribute; however, in practical SBCM engineering, there are usually several key strata, and their deformations are not necessarily uniform subsidence, resulting in complicated load distribution on the key stratum. Consequently, the previous theory models cannot capture the stress and deformation of the key strata appropriately. In addition, the previous studies mainly focus on the deformation of overburden, but it is necessary to grasp the fracture cases of the key stratum such as the key stratum for water retention in overburden. Therefore, the stress distribution of key stratum, in overburden, should also be systematically studied.
In this study, a multilayer beam model of SBCM based on Winkler foundation theory is established, and the case with 2 key strata is solved. The subsidence, internal forces, stresses distribution, and failure characteristics of key strata are systematically analyzed under various elastic modulus of backfill, mining heights, and thicknesses of soft layers.

Mechanical Mole of Overburden Strata in SBCM
In longwall mining, the overburden strata movement and deformation along tilt direction of the working face are the same, except both ends of working face. Therefore, the overburden strata mechanical model in the central stope can be simplified to a plane strain problem, shown as Figure 1.
Using the filled mining technology to replace the coal resources with filled materials which has adequate filling rate and compaction, the movement of overburden only induces slight bending deformation. Consequently, the overburden under solid filling mining only has continuous deformation, rather than discontinuous deformations such as the fracture of rock stratum and the development of fracture zone above the stope. That is, the theory of continuum mechanics can be employed to study the overburden deformation and land subsidence induced by solid filling mining.
To study the deformation characteristics of the overburden in longwall mining, the overburden structure can be simplified to a laminated beam on elastic foundation, shown in Figure 2(a). Moreover, based on the symmetry, only half of the model is considered, shown in Figure 2(b). The boundary condition of the model is that the horizontal displacement of left side is fixed, and both the horizontal and vertical displacement of the infinity at right side are fixed. For the case that the depth of coal seam is relatively shallow, it is more probably to have only one key stratum in the overburden. Because of this, the simplified single-layer beam on elastic foundation model, which is easily to be solved through the basic beam theory, could be used to study the overburden deformation. However, the above single beam model is inapplicable for the cases that the overburden has more than one key stratum. Therefore, a new elastic foundation beam model including more than one beams should be developed. Based on the key stratum theory, the overburden can be classified as key stratum and soft layer. In the model shown in Figure 2, the soft layer can be simplified to elastic foundation because of the small deformation in SBCM. Consequently, the model in Figure 2 can be further simplified to a multilayer Winkler foundation beam model, shown in Figure 3. In this model, all of the overburden key strata are simplified to beams, and the backfill, coal, and soft layer are simplified to elastic foundation. K 1 -K n denote the elastic foundation coefficients of soft layer 1 to n, K 11 and K 12 represent the elastic foundation coefficients of backfill and coal, respectively.  Figure 4. In this model, the two key strata are simplified to beams, while the backfill, coal, and soft layer are simplified to elastic foundation. K 11 , K 12 , and K 2 are the elastic foundation coefficients of backfill, coal, and soft layer, respectively. q is the weight of the rock and topsoil layers. In the actual geological conditions, the soft layer normally consists of several soft strata, resulting in that K 2 is an equivalent parameter, which could be calculated as follows:

Model Solution
where k i and h i are the foundation coefficient and thickness of layer i in the soft layers, respectively.

Geofluids
According to the multilayer beam model on Winkler foundation, the foundation coefficient K has a significant influence on overburden deformation, so it is necessary to solve the Winkler elastic foundation coefficient. According to the Winkler hypothesis, the elastic foundation coefficient of rock stratum can be calculated as follows: where E is the elastic modulus of the rock strata.
In Winkler foundation beam theory, the deflection at any point on the foundation surface is proportional to the pressure on the unit area of the point.
where PðxÞ is the support force of the foundation to the beam and wðxÞ is the deflection equation of the beam. The relation between wðxÞ, load, and PðxÞ is

Geofluids
Here, qðxÞ and EI are the load concentration and flexural rigidity of the beam, respectively. According to Eq. (4), the governing equation of the double-layer Winkler foundation beam model can be obtained as follows: where ω ij ðxÞ is the deflection equation of key stratum i, subscript j = 1 represents the part of the key stratum above the stope, and j = 2 represents the part of the key stratum above the unmined coal seam. q 1 , q 2 , and q soft are the weight of key stratum #1, key stratum #2, and soft layer, respectively, which determined by equationsq 1 = ρ 1 gh 1 , q 2 = ρ 2 gh 2 , and q soft = ρ soft gh soft , in which ρ and h are the densities and thicknesses of the rock strata. By adding Eqs. (5) and (6), it can be obtained that where q pl = q + q 1 + q 2 + q soft . The fourth derivative of both sides of Eq. (5) with respect to x can be obtained as follows: Substitute Eq. (10) into Eq. (9).
The above equation is an 8-order linear inhomogeneous differential equation with constant coefficients, of which the solution is the general solution of secondary differential equation plus a particular solution of the inhomogeneous equation. For the secondary differential equation of Eq. (11), its characteristic equation is as follows: The above equation can be converted to Substituting A and B in Eq. (11) into Eq. (13), it is easy to get M > 0; then, It is obvious that ffiffiffiffi ffi M p − B/2A < 0, the solution of the characteristic (Eq. (13)) can be obtained as follows: Key stratum #1 (BL) Boundary line of the stope Axis of symmetry: middle of stope

Geofluids
The real and imaginary part of Eq. (15) are denoted as real i and imag i , respectively, and combining with the Euler formula, the solution of Eq. (11) is as follows: where q pl /k 11 is a particular solution of Eq. (11). Based on the relation between ω 11 and ω 21 in Eq. (5), ω 21 can be expressed as follows: Equation (16) and Eq. (17) are the deflection equations of the part of key strata #1 and #2 above the goaf. In the same way, ω 12 and ω 22 can also be obtained. ð18Þ Equations (16)- (19) are the analytical solutions of the double-layer Winkler foundation beam model for solid backfill coal mining, but there are still 16 integral constants to be determined. Firstly, according to the boundary condition, the key strata at the infinitely far from stope are not affected by mining, and its bending moment and rotating angle are zero.
Combined with Eqs. (15), (18), and (19), it can be easily obtained that C 9 = C 12 = C 13 = C 16 = 0. In addition, the rotating angle and shearing force are 0 when x = 0, and the deflection, rotating angle, bending moment, and shearing force are continuous when x = L 1 ; then, : where θ, M, and Q are the rotating angle, bending moment, and shearing force of the key strata, respectively, and the relations of them with deflection are as follows: Based on the above equations, Eq. (21) can be used to obtain the linear system of equations with respect to the remaining 12 integral constants. Furthermore, the expressions for ω 1 ðxÞ and ω 2 ðxÞ can be obtained. The details of expressions are not listed here since they are too long.

n-Layer Winkler Foundation Beam Model.
For the multilayer Winkler foundation beam model that with n key strata in SBCM, the governing equation can be easily obtained according to Eqs. (5)- (8).
5 Geofluids By solving the double-layer Winkler foundation beam model, it is concluded that the key process is to get the solution of 8-order differential equation. For a three-layer beam, the order of the differential equation will increase to 12. The increase in number of layers will increase the order of the equations, resulting that is difficult to get the analytical expression of the key strata deformation of overburden. In order to solve the above equations, the corresponding program should be compiled with the help of mathematical calculation software and the numerical method should be applied.

Result and Discussion
According to the solution of the double-layer Winkler foundation beam model, the deformation, stress and inner force of the key strata are systematically analyzed under various elastic modulus E fill , mining height h, soft layer thickness h s , the calculation parameters, and schemes are shown in Tables 1 and 2. Based on the mechanical properties of backfill [3,6,20], the E fill is determined as 80-160 MPa. Figures 5 and 6 show the curves of key strata subsidence w and rotating angle θ with the distance to the middle of the stope (x) under various elastic modulus of backfill, mining heights, and thicknesses of soft layers. It can be seen that the subsidence curves of overburden key strata under all cases exhibit the "basin" shaped, which remains almost unchanged in a large range in the middle of stope, but decreases sharply near the BL (boundary line of the stope, Figure 4). This characteristic reflects that its rotating angle firstly rises and then falls with the increase of distance x from the middle of stope.

Deformation of the Key Strata.
With the increase of the elastic modulus of backfill, the peak value of subsidence of key stratum #1 decreases significantly. The range of subsidence curve "basin" gradually increases, while the horizontal scope of rotating angle that dramatically rises decreases, indicating that the role of filling material on supporting key stratum #1 is gradual significant. However, with the increase of the elastic modulus of the backfill, the increase of the "basin" range of the key stratum #2 is not significant, and its rotation-angle starts to increase at x = 110 m under all of the E fill . It is indicated that the elastic modulus of the backfill only has a great influence on the bending deformation of the key stratum #2, while slight influence on the location of bending.
With an increase in mining height, the subsidence of the key strata increase gradually. While the "basin" area decreases gradually, and the range of rotating angle increases gradually, which indicates that the greater the mining height is, the more significant the bending deformation of the overburden key stratum is. With the increase of the thickness of the soft layers, both the subsidence and "basin" range of key stratum #1 gradually rise, and the range of the rotating angle also increases. However, the above characteristic of key stratum #2 is opposite to that of key stratum #1. In addition, by comparing the subsidence curves of key strata #1 and #2, it can be seen that the subsidence of the two strata is almost the same in a larger range in the middle of the stope, indicating that the subsidence in this range is caused by the compressive deformation of backfill and the vertical downward movement only occurs in the 2 key strata and the soft layers. Figure 7 shows the bending moment curves of key strata with the distance to the middle of the stope under various elastic modulus of backfill, mining height, and thickness of soft layers, where points A and B represent the location of which bending moment is 0 in key strata #1 and #2, respectively. As can be seen that the bending moment curve shapes of the two key strata are analogical, the bending moment is almost 0 in the location of the key strata that far from the BL while the   The E fill and h s have a significant effect on the bending moment, which is relatively slightly influenced by h. Figure 8 shows the shear force curves of key strata with the distance to the middle of the stope under various elastic modulus of backfill, mining height, and thickness of soft layers. It can be seen that the shear force distributions of the two key strata are obvious different. For the key stratum #1, the shear force peak locates at BL, where the shear force has a sudden change, caused by the concentration force at   Geofluids the location of x = 200 due to greater stiffness difference between backfill and coal. For the key stratum #2, the shear force peak is located to the left of BL, which is significantly lower than that of key stratum #1. Moreover, with the increase of the elastic modulus of the backfill, the shear force peak of the two key strata decreases. When E fill increases from 60 MPa to 140 MPa, the shear force peak of the key stratum #1 decreases by 38.5% and #2 decreases by 52.4%. Therefore, it is effective to prevent the shear fracture of the key strata by increasing the elastic modulus of the backfill. In addition, with the increase of the thickness of the soft layer, the shear force of key stratum #2 decreases gradually, but the shear peak of key stratum #1 increases significantly. Consequently, it should be fully considered in engineering design that the thickness of the soft layer is significant on the shear failure of key stratum #1. Besides, mining height has little influence on the shear force of key stratum #1, and the shear peak of key stratum #2 changes moderately.

Stress of Overburden Key Strata.
For rectangular beam, bending moment and maximum normal stress, shear force, and maximum shear stress have the following relations: where W and A are the section modulus in bending and section area, respectively. According to the above equations, it is easy to obtain the absolute value of maximum horizontal stress and shear stress distribution of the two key strata, shown in Figures 9 and 10. It can be seen that distributions of the horizontal stress in the two key strata are similar, and the maximum values are located on the right of BL. Therefore, in the engineering design of solid backfill mining, the key strata tensile fracture at right of BL should be preferentially considered, because of the lower tensile strength of the rock material. For the maximum shear stress, the position of key stratum #1 is located on BL, while the position of key stratum #2 is located on the left of BL. Consequently, the shear strength of these two positions should attract more attentions in the design of solid backfill mining. In addition, with the increase of the elastic modulus of the backfill, the maximum horizontal stress and shear stress of the two key strata decrease significantly, indicating that improving the compaction of the backfill can effectively prevent the fracture of the key strata. However, with the increase of mining height, the maximum horizontal stress and shear stress of key stratum #1 change slightly, but the ones of key strata #2 increase significantly, and this characteristic denotes that special attention should be paid to the fracture of key stratum #2 in large mining height backfill engineering. Besides, with the increase of the thickness of the soft layer, the peak value of the horizontal stress and the shear stress of the key stratum #1 increase significantly, while the peak value of key stratum #2 slightly decreases. Therefore, it should attract more attentions to the fracture of the key stratum #1 with larger thickness of the soft layer in the engineering geological condition.

Geofluids
The vertical stress distribution is also an important aspect to analyze the damage of underground coal mining. For the Winkler foundation beam, the upper surface is subjected to vertical downward distributed load, while the lower surface is subjected to vertical upward distributed load due to the support of the foundation. Therefore, the foundation beam is generally under compressive stress in the vertical direction. In order to analyze the vertical stress distribution characteristics of the beam, a microelement dx in the beam is taken into account. The microelement is cut by a cross section that is perpendicular to y direction, and the part of the microelement above the cross section is taken as a research object, of which the force analysis diagram is shown in Figure 11. According to the equilibrium conditions, the following equation can be obtained.
Based on the distribution characteristics of shear stress in the beam section in bending, it can be obtained that Combining Eqs. (25) and (26), the vertical stress can be expressed as follows: According to Eqs. (5), (7), and (22), the vertical stress can be expressed withpðxÞ, qðxÞ, and y.
Equation (28)   11 Geofluids stresses of key strata #1 and #2 under various E fill , h, and h s can be easily obtained, of which the distribution is shown in Figures 12 and 13.
For key stratum #1, it can be seen that the maximum and minimum values of vertical stresses are located on both sides of BL, and the vertical stress of key stratum #1 changes sharply at this location because of the sudden change in shear force.
For key stratum #2, differing from key stratum #1, there is no mutation value. The maximum vertical stress is located on the right of BL. In addition, the vertical stress peak of key stra-tum #1 is significantly greater than that of key stratum #2. With the increase of the elastic modulus of the backfill, the vertical stress peaks, located on the bottom, of the two key strata decrease significantly. With the increase in mining height, the peak of vertical stress in key stratum #1 almost unchanged, while the one of key stratum #2 gradually increases. With the increase of thickness of soft layer, the vertical stress of both key strata increases. Within the scope of the model, the increase in the load on the top of key strata #1 results in increasing the support force from soft layer to key stratum #2.

Strength Analysis of Overburden Key Strata
Generally, rock material is typical brittle material, of which failure forms include tensile failure and shear failure. Tensile failure is generally determined by the maximum tensile stress criterion, while shear failure is determined by the More-Coulomb model. According to the horizontal stress and shear stress distribution of the key strata, their maximum values are distributed in different locations. Besides, stresses of different locations on the same cross section are also different on the basis of the stress distribution of rectangular section beam under transverse bending. Consequently, it is necessary to study the stress state in different locations of the cross section which has the maximum tensile stress and shear stress for further analyzing the strength of the key strata. Figure 14 shows the section positions of the maximum tensile stress and shear stress of the two key strata. Where points a, b, and c represent the top, middle, and bottom of the beam at the same section, respectively. According to the stress distribution characteristics of the key strata, the horizontal stress, shear stress characteristics, and stress states of each section are shown in Figure 15. It can be seen that the most dangerous point in the S1-S1 section of the key stratum #1 is point b, because the horizontal tensile stress is low at point c of S1-S1 section. The horizontal stress of point b in section S1-S1 is 0, which is mainly affected by shear stress and vertical stress. According to the More-Coulomb criterion, the τ max should meet the following requirements:  where c is the cohesive force of the key stratum material, n is the safety factor, and ½c is the allowable cohesive force. In addition, the dangerous location of section T1-T1 is at the point c, and its strength should meet the maximum tensile stress criterion.
where σ t is the tensile strength of the key stratum material and ½σ t is the allowable tensile stress. For key stratum #2, its strength condition is the same with that of key stratum #1 according to its stress distribution and state.
To summarize, the strength condition of key strata in solid backfill mining is as follows:

Conclusion
(1) A multilayer Winkler foundation beam model of the overburden strata in SBCM is established, and its analytical solution is obtained. According to the theoretical analysis results, the subsidence of key strata in overburden shows a "basin"-shape curve, and all of the backfill elastic modulus, mining height, and soft layer thickness have significant influences on the subsidence of key strata (2) The maximum shear stress of the key stratum close to the coal seam is located on BL, while the one of the key stratum above the soft strata is located on the cross section that horizontal position of which is on the right of BL. The maximum horizontal stresses of both key strata are located to the cross section that horizontal position of which is on the right of BL (3) The shear stress, horizontal stress, and vertical stress of both the two key strata can be effectively reduced by increasing the elastic modulus of backfill. The increase of mining height has negligible influence on the shear and vertical stresses of key stratum #1, but has a significant influence on the all of the stresses of key stratum #2. The increase of the soft layer thickness has remarkable influence on the horizontal and shear stress of key stratum #1, but has slight influence on the stress of key stratum #2 (4) In the engineering design of SBCM, the shear failure of key stratum #1 at BL and the tensile failures on both sides of BL should be preferentially considered. Due to the low shear stress level of key stratum #2, the tensile failure on both sides of BL should be mainly considered

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.