# Nonlinear Elastic Foundation Plate Model Study on Overlying Strata Movement in Working Face with Solid Backfilling Mining

Authors: Baifu An; Nailu Li; Qiaomei Yi; Dong Zhang; Hailong Wang
Journal: Shock and Vibration
 (2021)

---

## Abstract

Although solid backfilling materials are featured with obvious nonlinear stress-strain properties, for a long time, they have been usually simplified as linear elastic materials for approximate calculation in mechanical analysis, so it is difficult to accurately reflect their deformation process. Based on test results of solid backfilling materials’ compaction characteristics, this paper provides a solution method to generate their elastic foundation coefficient. One multiparameter elastic foundation has been used to reflect stress-strain characteristics of solid backfilling material. In addition, the paper establishes a thin plate on a nonlinear elastic foundation model by adopting semianalytical and seminumerical method and obtains the relational expression between roof deflection, roof stress, and backfilling material’s compressive deformation. In combination with geological conditions in a specific mine, the paper probes into what influence both backfilling material’s particle size and the initial compaction force that the backfilling material bears could exert on roof subsidence and stress. Finally, the proposed model has been verified with measured data from industrial tests.


---

## Body

## 1. Introduction

China’s economic development is highly dependent on coal resources. As a matter of fact, coal accounts for over 56.8% of China’s energy consumption system and the demand is still growing. Especially since the 21st century, there has been a surge in demand for coal, which has grown from 1.1 billion tons in 2001 to 3.9 billion tons in 2020. Due to enormous demand, the magnitude of coal mining has been multiplied and recoverable reserves in numerous mines (especially in eastern China) have plummeted as well. Faced with this situation, many coal mines have to set their eyes on the excavation of the coal resources under buildings, water bodies, and railways. As a matter of fact, a great amount of valuable coal resources is being buried under buildings, water bodies, and railways; for instance, for large state-owned coal mines alone, the figure hits over 13.0 billion tons. Hence, how to effectively emancipate this coal resource has become one of the most important research topics in promoting sustainable development of mining area [1–3]. Meanwhile, gangue and other solid wastes are piling up on ground surface in various mining areas, which account for over 4.5 billion tons right now and are still rapidly increasing at an annual rate of 150 million tons. Their accumulation has taken up large areas of farmland and brought about many problems such as spontaneous combustion, explosion of coal gangue piles, and water resources pollution.Aimed at emancipating the coal resources and solving the problem of gangue accumulation on the ground surface, solid backfilling technology has been developed by backfilling solid waste like gangue into mined-out areas so as to effectively control strata movement and reduce surface deformation and ensure the safety of ground buildings [4–8]. Involved in the self-organization process of strata movement, the backfilling materials under roof pressure would undergo vertical compression and deformation [9–12]. Nevertheless, in the process of compression, the backfilling body would play an increasingly supportive role for the roof and hence share part of the overlying strata load. In this sense, strata movements are jointly controlled by backfilling material and hydraulic support as well as coal body which altogether constitute a new stope support system [13–16].In order to study the behavior of overlying strata in the practice of solid backfilling mining (SBM), scholars have established various relevant models, including but not limited to equivalent mining height model, key strata mechanics analysis model of combinatorial elastic base, linear elastic foundation plate mechanics model, roof deformation and energy change in surrounding rock calculation model, and simplified and coupling model for overlying rock and backfill mechanics [17–25]. Every single model has reflected overlying strata behavior characteristics in solid backfilling mining from a certain perspective and provides favorable theoretical support for stope pressure and surface subsidence control. Nevertheless, the strain processes of backfilling material in the above models have been all simplified as linear elastic process and hence failed to showcase nonlinear stress-strain characteristics of backfilling material in mined-out area. Based on test results analysis concerning backfilling materials’ mechanical properties, this literature establishes a nonlinear elastic foundation plate model using seminumerical and semianalytical method for the purpose of exploring characteristic of strata behavior in SBM. This study has also been verified in combination with industrial test.


## 2. Main Factors Influencing Overlying Strata Subsidence

In the case of SBM, a key factor in stope pressure control is whether the mined-out area could be filled to the maximum extent or not. When the mined-out area is fully filled, the possible strata subsidence would be reduced and the backfilling material would play a greater supportive role for the roof. To be specific, a goaf’s backfilling degree is largely determined by roof subsidence prior to backfilling, distance from backfilling body to the roof, and compression deformation of backfilling material.


### 2.1. Roof Subsidence prior to Backfilling

Coal mining upsets original rock stress balance and causes stress redistribution, which would inevitably result in overburden strata deformation and movement. Therefore, the exposed roof subsidence is sure to occur. In SBM,hx, the roof subsidence prior to backfilling, refers to the height that the roof subsides before packing materials are filled into the goaf, as shown in Figure 1. Unfortunately, given the current technology, it is still not feasible to get the prior roof sinking restored to the original position, so the only option so far is to rely on an increased hydraulic support’s force so as to limit prior roof subsidence as much as possible.Figure 1
Overlying strata subsidence in solid backfilling and mining.


### 2.2. Distance from the Backfilling Body to the Roof

Since granular backfilling materials have certain flowability, they slide easily once they are filled into the mined-out area, thus resulting in an interspace between solid backfilling material piling up in the gob and the roof. The vertical distance of the interspace is defined ashq, as shown in Figure 2. Obviously, higher hq indicates worse control effect of backfilling material for roof subsidence. Therefore, backfilling effect could be improved through optimizing the parameters of backfilling equipment so that it could adapt to different mining heights. In addition, strengthening backfilling operation management also contributes to lower hq.Figure 2
Distance from the backfilling body to the roof.


### 2.3. Compressive Deformation of Backfilling Material

The compressive resistance of backfilling material tends to differ markedly from one source to another. In actual production, considering cost control and backfilling material available in mining areas, there are only rather limited choices. Nevertheless, a couple of methods can still be applied to enhance their compressibility, such as modifying backfilling material gradation and using mixed material.As shown in Figure3, when backfilling materials are initially fed into the mined-out area, they still have a low degree of compaction and a high porosity. As overlying strata subside, they are further compacted, bringing about compression deformation. Hence, hk, the compressive deformation of backfilling material, is the difference between hc, the initial height of backfilling body, and hm, its final height after it is fully compacted by complete roof subsidence.Figure 3
Compressive deformation of backfilling body.(1)hk=σ2E2−σ1E1hc.Extensive research results demonstrated that the value ofhk is jointly determined by the following factors: σ1, the initial compressive stress that backfilling material bears; σ2, the ultimate compressive stress that backfilling bears; and the respective elastic moduli E1 and E2, as well as hc. The corresponding relationship amongst these factors is shown as follows:As shown in formula (1), we find that both higher σ1 and E2 and lower σ2 and E1 cause the decrease of compression deformation of backfilling material, amongst which E1 and E2 are decided by backfilling materials’ inherent nature, whereas σ2 represents the original rock stress, so it is difficult to change these three factors. In this case, the only way to decrease compressive deformation is to improve the initial compressive stress of the backfilling materials.Thanks to significantly elevated design level and machinery manufacturing level, hydraulic support and rammer compactor being used in SBM nowadays have made it possible to basically eliminate prior roof subsidence and the interspace from backfilling body to the roof. In this case, how to effectively decrease the compression deformation of backfilling material has become a determinant factor, which therefore would be our focus of research hereunder.


## 2.1. Roof Subsidence prior to Backfilling

Coal mining upsets original rock stress balance and causes stress redistribution, which would inevitably result in overburden strata deformation and movement. Therefore, the exposed roof subsidence is sure to occur. In SBM,hx, the roof subsidence prior to backfilling, refers to the height that the roof subsides before packing materials are filled into the goaf, as shown in Figure 1. Unfortunately, given the current technology, it is still not feasible to get the prior roof sinking restored to the original position, so the only option so far is to rely on an increased hydraulic support’s force so as to limit prior roof subsidence as much as possible.Figure 1
Overlying strata subsidence in solid backfilling and mining.


## 2.2. Distance from the Backfilling Body to the Roof

Since granular backfilling materials have certain flowability, they slide easily once they are filled into the mined-out area, thus resulting in an interspace between solid backfilling material piling up in the gob and the roof. The vertical distance of the interspace is defined ashq, as shown in Figure 2. Obviously, higher hq indicates worse control effect of backfilling material for roof subsidence. Therefore, backfilling effect could be improved through optimizing the parameters of backfilling equipment so that it could adapt to different mining heights. In addition, strengthening backfilling operation management also contributes to lower hq.Figure 2
Distance from the backfilling body to the roof.


## 2.3. Compressive Deformation of Backfilling Material

The compressive resistance of backfilling material tends to differ markedly from one source to another. In actual production, considering cost control and backfilling material available in mining areas, there are only rather limited choices. Nevertheless, a couple of methods can still be applied to enhance their compressibility, such as modifying backfilling material gradation and using mixed material.As shown in Figure3, when backfilling materials are initially fed into the mined-out area, they still have a low degree of compaction and a high porosity. As overlying strata subside, they are further compacted, bringing about compression deformation. Hence, hk, the compressive deformation of backfilling material, is the difference between hc, the initial height of backfilling body, and hm, its final height after it is fully compacted by complete roof subsidence.Figure 3
Compressive deformation of backfilling body.(1)hk=σ2E2−σ1E1hc.Extensive research results demonstrated that the value ofhk is jointly determined by the following factors: σ1, the initial compressive stress that backfilling material bears; σ2, the ultimate compressive stress that backfilling bears; and the respective elastic moduli E1 and E2, as well as hc. The corresponding relationship amongst these factors is shown as follows:As shown in formula (1), we find that both higher σ1 and E2 and lower σ2 and E1 cause the decrease of compression deformation of backfilling material, amongst which E1 and E2 are decided by backfilling materials’ inherent nature, whereas σ2 represents the original rock stress, so it is difficult to change these three factors. In this case, the only way to decrease compressive deformation is to improve the initial compressive stress of the backfilling materials.Thanks to significantly elevated design level and machinery manufacturing level, hydraulic support and rammer compactor being used in SBM nowadays have made it possible to basically eliminate prior roof subsidence and the interspace from backfilling body to the roof. In this case, how to effectively decrease the compression deformation of backfilling material has become a determinant factor, which therefore would be our focus of research hereunder.


## 3. Elastic Foundation Coefficient Calculations of Solid Backfilling Materials

### 3.1. Stress-Strain Relations of Backfilling Material under Confined Uniaxial Compression Test

Confined uniaxial compression tests were conducted to study the stress-strain characteristics of different backfilling materials, using a chamber. The maximum compressive strength is 10.0 MPa with a loading rate of 2 kN/s. However, while the value of stress is greater than 10.0 MPa, the change rules of stress and strain can be found rapidly and accurately. Before testing, materials were placed in a compaction steel cylinder (Figure4) and the initial height h of material was recorded. During the compaction, axial stress and the corresponding axial displacement were monitored and collected by using a data acquisition system.Figure 4
YAS-5000 material test system and compaction steel cylinder: YAS-5000 electric servo hydraulic material test system (a); compaction steel cylinder (b); compactor (c).
(a)(b)(c)


### 3.2. The Computing Method of Elastic Foundation Coefficient

Elastic foundation coefficient of backfilling material in the gob could be obtained by the following steps:(a)
Get specimen deformation modulus according to the following formula:(2)E0=σε,whereE0 refers to deformation modulus, σ represents axial compressive strength, and ε is the corresponding axial strain.(b)
Get specimen’s elastic modulus:(3)E=dσdε,wheredσ represents the axial compressive strength differential, while dε represents the axial strain differential.(c)
According toE=fε, namely, the relation expression between elastic modulus and strain, we could deduce the relationship between elastic foundation coefficient parameters and strain, which reflects the compressive capacity of backfilling body. Thus, the elastic foundation coefficient of backfilling materials in steel cylinder can be obtained according to the following formula:(4)Kt=σhε=Eh,whereKt is elastic foundation coefficient; σ is loading stress; h is instantaneous height during compression; ε is instantaneous deformation.(d)
Formula (5) indicates the relation between elastic foundation coefficient and the strain of the backfilling materials specimen.(5)Kt=fεh=fε1−εH,whereKt refers to the elastic foundation coefficient of materials in steel cylinder; f(ε) is function relation about the strain; H is the original height of backfilling material.(e)
To further study sample materials’ compressive capacity in the compression process, it is defined that elastic foundation coefficient parameterK’ = KH; then its relation with strain can be expressed by the following formula:(6)K′=fε1−ε,whereK represents the elastic foundation coefficient parameter.(f)
The elastic foundation coefficient in the stope could be derived from the following formula:(7)Kc=K′hc,whereKc is the elastic foundation coefficient in the stope and hc is the backfilling height in stope.


## 3.1. Stress-Strain Relations of Backfilling Material under Confined Uniaxial Compression Test

Confined uniaxial compression tests were conducted to study the stress-strain characteristics of different backfilling materials, using a chamber. The maximum compressive strength is 10.0 MPa with a loading rate of 2 kN/s. However, while the value of stress is greater than 10.0 MPa, the change rules of stress and strain can be found rapidly and accurately. Before testing, materials were placed in a compaction steel cylinder (Figure4) and the initial height h of material was recorded. During the compaction, axial stress and the corresponding axial displacement were monitored and collected by using a data acquisition system.Figure 4
YAS-5000 material test system and compaction steel cylinder: YAS-5000 electric servo hydraulic material test system (a); compaction steel cylinder (b); compactor (c).
(a)(b)(c)


## 3.2. The Computing Method of Elastic Foundation Coefficient

Elastic foundation coefficient of backfilling material in the gob could be obtained by the following steps:(a)
Get specimen deformation modulus according to the following formula:(2)E0=σε,whereE0 refers to deformation modulus, σ represents axial compressive strength, and ε is the corresponding axial strain.(b)
Get specimen’s elastic modulus:(3)E=dσdε,wheredσ represents the axial compressive strength differential, while dε represents the axial strain differential.(c)
According toE=fε, namely, the relation expression between elastic modulus and strain, we could deduce the relationship between elastic foundation coefficient parameters and strain, which reflects the compressive capacity of backfilling body. Thus, the elastic foundation coefficient of backfilling materials in steel cylinder can be obtained according to the following formula:(4)Kt=σhε=Eh,whereKt is elastic foundation coefficient; σ is loading stress; h is instantaneous height during compression; ε is instantaneous deformation.(d)
Formula (5) indicates the relation between elastic foundation coefficient and the strain of the backfilling materials specimen.(5)Kt=fεh=fε1−εH,whereKt refers to the elastic foundation coefficient of materials in steel cylinder; f(ε) is function relation about the strain; H is the original height of backfilling material.(e)
To further study sample materials’ compressive capacity in the compression process, it is defined that elastic foundation coefficient parameterK’ = KH; then its relation with strain can be expressed by the following formula:(6)K′=fε1−ε,whereK represents the elastic foundation coefficient parameter.(f)
The elastic foundation coefficient in the stope could be derived from the following formula:(7)Kc=K′hc,whereKc is the elastic foundation coefficient in the stope and hc is the backfilling height in stope.


## 4. Modelling of a Thin Plate on a Nonlinear Elastic Foundation

The bending of roof is the root cause for coal pillar and backfilling materials deformation, since the loading of overlying strata is transferred to coal pillars and backfilling materials through roof. Strata behaviors in stope are determined by the coupled action of backfilling materials, coal pillar, and roof. Therefore, roof bending deformation has been the research object of the literature, which would conduct further analysis on the characteristics concerning the pillar stress, backfilling materials deformation, and stress distribution according to the coordination relations mentioned above.As shown in Figure5, the mechanical model illustrating the interaction amongst the backfilling body, coal body, hydraulic support, and the roof is shown in Figure 4. Considering the fact that the thickness of the roof is far less than the other two dimensions (i.e., length and width), it can be simplified as a thin plate that is fixed on four sides XX. Specifically, the upper surface undertakes a uniformly distributed vertical load of q that is perpendicular to the transverse. Coal pillars are treated as homogeneous elastomer XX and their support reaction force for the roof is proportional to their deformation; that is, pc=kcwc, where kc is the elastic foundation coefficient of coal body and wc is the coal deformation. As far as the backfilling body is concerned, it is nonuniform elastomer (will be discussed later below) and, hence, there exists nonlinear relationship between its support reaction force for the roof and its deformation. In addition, Pg=k1+k2w+k3w2+k4w3+k5w4+⋯, where k1,k2,k3,k4,k5 are elastic foundation coefficients of the backfilling body.Figure 5
Mechanical model of the main roof.


### 4.1. Mechanical Model of the Main Roof

Force bearing analysis has been conducted on the roof whose length isD and finite element method has been used to recover the coal pillars below. The roof size is A+2B×C+D+E+F and the backfilling body size is A×D. In general, B, C, and D should be not less than 45 m. Furthermore, the middle face of coal seam roof has been divided into 12 areas and we did discretization using eight different sizes of quadrilateral element, thus producing a discrete roof system. As shown in Figure 6, the roof is partitioned by quadrilateral elements of six different dimensions, including a total of m=n1+2n3×n2+n4+n6+1 quadrilateral elements and the according node number is m=n1+2n3+1×n2+n4+n6+2.Figure 6
Roof mesh generation above the workface.


### 4.2. Quadrilateral Thin Plate Element Analysis

One element was chosen randomly from those above-mentioned for analysis. This method applies to the other 7 kinds of elements as well.


#### 4.2.1. Quadrilateral Thin Plate Element Displacement Mode

Under the only effect of bending moment, the roof would have a displacement ofω solely in the direction of z-axis whose rotation angles with ε-axis and y-axis are, respectively, θε and θy. In addition, θx=∂ω/∂y and θy=−z∂ω/∂x. Since all elements’ deflection and slope at the nodes are of continuity, the deflection and its first-order partial derivatives in the direction of ε-axis and y-axis have to be designated as nodal displacement (called generalized displacement). Take a second kind of thin plate element and establish a local coordinate system. Figure 7 is an illustration of relevant node displacement. In this element, there are three displacement components in each node, namely, the according node freedoms with ω being the displacement and θε and θy being the rotation angles. Each element has 4 nodes with a total of 12 nodal displacement components, thus producing an element freedom degree of 12.Figure 7
Thin plate rectangular element and its nodal displacement diagram.At the node ofi, the matching displacement component and nodal force component can be represented as(8)δi=ωiθxiθyi,Fi=WiMθxiMθyi.Nodal displacement of the quadrilateral thin plate element is(9)δe=ωi,θxi,θyi,ωj,θxj,θyj,ωk,θxk,θyk,ωl,θxl,θylT.Generally speaking, the deflectionω and the matching nodal force W are positive along the forward axis direction. Likewise, for θε and θy and the matching moments Mθx and Mθy, the vectors obtained in accordance with right-hand rule are positive along the axis direction.Select a polynomial containing 12 parameters as the displacement mode for plate element, namely, the interpolation polynomials concerning its displacement and rotation angle, as shown in the following formulas:(10)ω=a1+a2x+a3y+a4xy+a5x2+a6y2+a7x2y+a8xy2+a9x3+a10y3+a11x3y+a12xy3,(11)θx=∂ω∂y=a3+a4x+2a6y+a7x2+2a8xy+3a10y2+a11x3+3a12xy2,(12)θy=−z∂ω∂x=−a2+a4y+2a5x+2a7xy+a8y2+3a9x2+3a11x2y+a12y3.When the above equations are expressed in matrix form, the displacement mode is(13)ωθxθy=1xyxyx2y2x2yxy2x3y3x3yxy3001x02yx22xy03y2x33xy20−10−y−2x0−2xy−y2−3x30−3x2y−y3a1⋮a12.When we plug the coordinates and corresponding displacement components of the four nodes into formula (13) above, we would get(14)ωiθxiθyiωjθxiθyjωkθxkθykωlθxlθyl=1xiyixiyixi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyl2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3a1a2a3a4a5a6a7a8a9a10a11a12.By solving the equation above, we get a coefficient vector ofa1…a12T. Plug it into the displacement mode equation and we get an expression using a shape function matrix of N.(15)ω=Nδe,where the shape function matrix of N is as follows:(16)N=N1N2N3N4N5N6N7N8N9N10N11N12=1xyxyx2y2x2yxy2x3y3x3yxy31xiyixiylxi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyk2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3.


#### 4.2.2. Stiffness Matrix of Quadrilateral Thin Plate Element

Since there are 12 degrees of freedom in the 4-node plate element, the stiffness matrix is a square matrix of12×12. The geometric equation of bending plate is(17)ε=εxεyγxy=Bδe,where B is a strain matrix of 3×12.(18)B=−Z∂2∂x2∂2∂y22∂2∂x∂yN.The strain energy of the plate element is(19)U=∫V12εTσdV=∫V12εTDεdV.Using the principle of minimum potential energy, we get(20)∂U∂δe=Kδe=∫BTDBz2dVδe.Therefore, the solving equation for the element stiffness matrix is(21)k=∫BTDBz2dV=∫∫h/2h/2BTDBz2dzdA,where D is the elastic matrix for materials:(22)D=E1−μ21μ0μ10001−μ2.Because the thickness of the plate element is uniform, we get(23)∫−h/2h/2z2dz=h312.Based on the above equation, we get stiffness matrix(24)k=∫BTDBz2dV=h3121−μ2∫−4.54.5dx∫−4.54.5BTDBdy.By combining all above-mentioned equations, the element stiffness matrix was solved using MATLAB programming. However, due to the complexity of solution results, they are not specifically listed hereunder. In addition, the stiffness matrix of elements of the other four sizes was solved using the same methodology, which will not be presented to keep the paper at a reasonable length.


### 4.3. Element Stiffness Matrix Assembly

Having solved stiffness matrix of the above six elements, we need to solve the overall stiffness matrix. On the entire roof, there arem elements and n nodes. Hence, the dimension of overall stiffness matrix is 3n1+2n3+1×3n2+n4+n6+2. Furthermore, element node freedom degrees are introduced to correspond to transition matrix G that has been expanded to be a structure element node. Suppose that there are a total of n structure nodes and the overall node numbers corresponding to quadrilateral thin plate elements are m, r, s, and t; then the transformation matrix for the element’s nodal freedom degrees is as follows:(25)That is, in the transformation matrix ofG, the subblocks where the overall serial numbers (m, r, s, and t) are located are set as third-order unit matrix, and all the others are zero. Using the method of transformation matrix, according to node numbering of each unit, the corresponding transformation matrix is G1,G2,G3,Gm. Since each transformation matrix is of 12×3n, the total stiffness matrix is(26)K=∑i=1mGiTkiGi,i=1,2,3,…,m.The total stiffness matrix is of3n×3n, yet, due to its huge dimension, the solution results are not elaborated upon herein.


### 4.4. Statics Solution Equation

When the transverse distributed load that the element exerts on the middle surface perpendicularly isq, the unit load vector could be calculated using the following formula:(27)Pe=∫∫ΩNTqdxdy.In this plate element, the transverse distributed load is as follows:q1=q, q2=q−kcwc, and q3=q−k1+k2w+k3w2+k4w3+k5w4+⋯. When q1=q, since q is a constant, the specific outcomes for Pe are presented in formula (27), where the unit length and width are assumed, respectively, as 2a and 2b.(28)Pe=4qab14b12−a1214b12a1214−b12a1214−b12−a12T.Besides, with formula (21), we get q3=q−k1+k2w+k3w2+k4w3+k5w4+⋯ and q2=q−kcw. After obtaining node load of the unit, the overall roof node load array can also be obtained by transition matrix(29)P=∑i=1mGiTPe,i=1,2,3,…,m.The roof is dotted with a total ofn nodes; therefore, the overall nodal force matrix P is of 3n×1bit, while the overall node displacement array δ is of 3n×1bit. But, due to their huge dimension, their expressions are not being demonstrated herein. Therefore, the total statics solving equation is(30)P3n×1=K3n×3nδ3n×1.


### 4.5. Introduction of Boundary Conditions and Equation Solutions

Since the boundaries of the roof are assumed to be fixed, displacement and rotation at the boundaries are 0. On basis of these conditions, the overall stiffness matrix can be modified as described below. The number of nodes on the boundary isp=2n1+2n2+4n3+2n4+2n6−2. If the three rows and three columns at each corresponding node in the total stiffness matrix are removed, then the new stiffness matrix would be K3n−p×3n−p′ and the new array of node displacement would be δ3n×1′. Now let us bring in force boundary conditions and the nodal force at the other nodes is Pi; then we remove a number of p nodes off the boundary and get new nodal force of P3n−p×1′. Consequently, we get a new statics balance equation as follows:(31)P3n−p×1′=K3n−p×13n−p′δ3n−p×1′,the specific solution equation of which is(32)δ3n−p×1′=K3n−p×3n−p′−1P3n−p×1′.Programming calculation for the above process is accomplished using Femap and Nastran. Specifically, Newton iterative method is applied to solve the nonlinear equations herein because this method is characteristic of simple programming and quick convergence rate; what is more, its precision can generally meet the requirements of this calculation.Firstly, we transform nonlinear equations in formula (30) into formula (32) to solve ω1, θx1, θy1, …, ω(3n-p), ωx(3n-p), θy(3n-p). Assume a certain spot as z and we could obtain the linear terms and formula (33) according to Taylor expansion of multivariate function at the point of Z(k).(33)f1ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f2ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f3ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,⋯,f3n−p−1ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f3n−pω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0.(34)f1kω1k,θx1k,θy1k,…ω3n−pk,θx3n−pk,θy3n−pkf2kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf3kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk⋯fn−1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkfnkω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk+f′zkω1k+1−ωkθx1k+1−θx2kθy1k+1−θy2k⋯θx3n−p−1k+1−θxn−p−1kθy3n−pk+1−θy3n−pk=0,where(35)f′z=∂f1∂ω1⋯∂f1∂θy3n−p⋯∂f3n−p∂ω1⋯∂f3n−p∂θy3n−p.We get the iterative formula as follows:(36)ω1k+1θx1k+1⋯θ3n−p−1xk+1θ3n−pyk+1=ω1kθx1k⋯θ3n−p−1xkθ3n−pyk−f′zk−1⋅f1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf2kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf3kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk⋯fn−1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkfnkω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk.During iteration process, the function equation and partial derivatives matrix function are programmed in advance using FORTRAN language. When it comes to calculation process, the above formula would be applied constantly to do iterative computation. The results obtained are node displacement, which are not going to be listed herein due to their complexity. Then regenerate the general node displacement arrayδ3n×1 and, according to the initial static balance equation, we get nodal force at each node. Then we would acquire strain according to the geometric equations and stress according to the physics equation.We getMx, My, σx, andσy in corresponding model by solving formula (36) as follows:(37)M=MxMyMxy=∫−h/2h/2σz2dz=−h312D∂2w∂2x2∂2w∂2y22∂2w∂x∂y=−Eh3121−μ21μ0μ10001−μ/2∂2w∂2x2∂2w∂2y22∂2w∂x∂y,σ=σxσy=12M⋅zh2=12Mx⋅zh212My⋅zh2.


## 4.1. Mechanical Model of the Main Roof

Force bearing analysis has been conducted on the roof whose length isD and finite element method has been used to recover the coal pillars below. The roof size is A+2B×C+D+E+F and the backfilling body size is A×D. In general, B, C, and D should be not less than 45 m. Furthermore, the middle face of coal seam roof has been divided into 12 areas and we did discretization using eight different sizes of quadrilateral element, thus producing a discrete roof system. As shown in Figure 6, the roof is partitioned by quadrilateral elements of six different dimensions, including a total of m=n1+2n3×n2+n4+n6+1 quadrilateral elements and the according node number is m=n1+2n3+1×n2+n4+n6+2.Figure 6
Roof mesh generation above the workface.


## 4.2. Quadrilateral Thin Plate Element Analysis

One element was chosen randomly from those above-mentioned for analysis. This method applies to the other 7 kinds of elements as well.


### 4.2.1. Quadrilateral Thin Plate Element Displacement Mode

Under the only effect of bending moment, the roof would have a displacement ofω solely in the direction of z-axis whose rotation angles with ε-axis and y-axis are, respectively, θε and θy. In addition, θx=∂ω/∂y and θy=−z∂ω/∂x. Since all elements’ deflection and slope at the nodes are of continuity, the deflection and its first-order partial derivatives in the direction of ε-axis and y-axis have to be designated as nodal displacement (called generalized displacement). Take a second kind of thin plate element and establish a local coordinate system. Figure 7 is an illustration of relevant node displacement. In this element, there are three displacement components in each node, namely, the according node freedoms with ω being the displacement and θε and θy being the rotation angles. Each element has 4 nodes with a total of 12 nodal displacement components, thus producing an element freedom degree of 12.Figure 7
Thin plate rectangular element and its nodal displacement diagram.At the node ofi, the matching displacement component and nodal force component can be represented as(8)δi=ωiθxiθyi,Fi=WiMθxiMθyi.Nodal displacement of the quadrilateral thin plate element is(9)δe=ωi,θxi,θyi,ωj,θxj,θyj,ωk,θxk,θyk,ωl,θxl,θylT.Generally speaking, the deflectionω and the matching nodal force W are positive along the forward axis direction. Likewise, for θε and θy and the matching moments Mθx and Mθy, the vectors obtained in accordance with right-hand rule are positive along the axis direction.Select a polynomial containing 12 parameters as the displacement mode for plate element, namely, the interpolation polynomials concerning its displacement and rotation angle, as shown in the following formulas:(10)ω=a1+a2x+a3y+a4xy+a5x2+a6y2+a7x2y+a8xy2+a9x3+a10y3+a11x3y+a12xy3,(11)θx=∂ω∂y=a3+a4x+2a6y+a7x2+2a8xy+3a10y2+a11x3+3a12xy2,(12)θy=−z∂ω∂x=−a2+a4y+2a5x+2a7xy+a8y2+3a9x2+3a11x2y+a12y3.When the above equations are expressed in matrix form, the displacement mode is(13)ωθxθy=1xyxyx2y2x2yxy2x3y3x3yxy3001x02yx22xy03y2x33xy20−10−y−2x0−2xy−y2−3x30−3x2y−y3a1⋮a12.When we plug the coordinates and corresponding displacement components of the four nodes into formula (13) above, we would get(14)ωiθxiθyiωjθxiθyjωkθxkθykωlθxlθyl=1xiyixiyixi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyl2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3a1a2a3a4a5a6a7a8a9a10a11a12.By solving the equation above, we get a coefficient vector ofa1…a12T. Plug it into the displacement mode equation and we get an expression using a shape function matrix of N.(15)ω=Nδe,where the shape function matrix of N is as follows:(16)N=N1N2N3N4N5N6N7N8N9N10N11N12=1xyxyx2y2x2yxy2x3y3x3yxy31xiyixiylxi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyk2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3.


### 4.2.2. Stiffness Matrix of Quadrilateral Thin Plate Element

Since there are 12 degrees of freedom in the 4-node plate element, the stiffness matrix is a square matrix of12×12. The geometric equation of bending plate is(17)ε=εxεyγxy=Bδe,where B is a strain matrix of 3×12.(18)B=−Z∂2∂x2∂2∂y22∂2∂x∂yN.The strain energy of the plate element is(19)U=∫V12εTσdV=∫V12εTDεdV.Using the principle of minimum potential energy, we get(20)∂U∂δe=Kδe=∫BTDBz2dVδe.Therefore, the solving equation for the element stiffness matrix is(21)k=∫BTDBz2dV=∫∫h/2h/2BTDBz2dzdA,where D is the elastic matrix for materials:(22)D=E1−μ21μ0μ10001−μ2.Because the thickness of the plate element is uniform, we get(23)∫−h/2h/2z2dz=h312.Based on the above equation, we get stiffness matrix(24)k=∫BTDBz2dV=h3121−μ2∫−4.54.5dx∫−4.54.5BTDBdy.By combining all above-mentioned equations, the element stiffness matrix was solved using MATLAB programming. However, due to the complexity of solution results, they are not specifically listed hereunder. In addition, the stiffness matrix of elements of the other four sizes was solved using the same methodology, which will not be presented to keep the paper at a reasonable length.


## 4.2.1. Quadrilateral Thin Plate Element Displacement Mode

Under the only effect of bending moment, the roof would have a displacement ofω solely in the direction of z-axis whose rotation angles with ε-axis and y-axis are, respectively, θε and θy. In addition, θx=∂ω/∂y and θy=−z∂ω/∂x. Since all elements’ deflection and slope at the nodes are of continuity, the deflection and its first-order partial derivatives in the direction of ε-axis and y-axis have to be designated as nodal displacement (called generalized displacement). Take a second kind of thin plate element and establish a local coordinate system. Figure 7 is an illustration of relevant node displacement. In this element, there are three displacement components in each node, namely, the according node freedoms with ω being the displacement and θε and θy being the rotation angles. Each element has 4 nodes with a total of 12 nodal displacement components, thus producing an element freedom degree of 12.Figure 7
Thin plate rectangular element and its nodal displacement diagram.At the node ofi, the matching displacement component and nodal force component can be represented as(8)δi=ωiθxiθyi,Fi=WiMθxiMθyi.Nodal displacement of the quadrilateral thin plate element is(9)δe=ωi,θxi,θyi,ωj,θxj,θyj,ωk,θxk,θyk,ωl,θxl,θylT.Generally speaking, the deflectionω and the matching nodal force W are positive along the forward axis direction. Likewise, for θε and θy and the matching moments Mθx and Mθy, the vectors obtained in accordance with right-hand rule are positive along the axis direction.Select a polynomial containing 12 parameters as the displacement mode for plate element, namely, the interpolation polynomials concerning its displacement and rotation angle, as shown in the following formulas:(10)ω=a1+a2x+a3y+a4xy+a5x2+a6y2+a7x2y+a8xy2+a9x3+a10y3+a11x3y+a12xy3,(11)θx=∂ω∂y=a3+a4x+2a6y+a7x2+2a8xy+3a10y2+a11x3+3a12xy2,(12)θy=−z∂ω∂x=−a2+a4y+2a5x+2a7xy+a8y2+3a9x2+3a11x2y+a12y3.When the above equations are expressed in matrix form, the displacement mode is(13)ωθxθy=1xyxyx2y2x2yxy2x3y3x3yxy3001x02yx22xy03y2x33xy20−10−y−2x0−2xy−y2−3x30−3x2y−y3a1⋮a12.When we plug the coordinates and corresponding displacement components of the four nodes into formula (13) above, we would get(14)ωiθxiθyiωjθxiθyjωkθxkθykωlθxlθyl=1xiyixiyixi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyl2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3a1a2a3a4a5a6a7a8a9a10a11a12.By solving the equation above, we get a coefficient vector ofa1…a12T. Plug it into the displacement mode equation and we get an expression using a shape function matrix of N.(15)ω=Nδe,where the shape function matrix of N is as follows:(16)N=N1N2N3N4N5N6N7N8N9N10N11N12=1xyxyx2y2x2yxy2x3y3x3yxy31xiyixiylxi2yi2xi2yixiyi2xi3yi3xi3yixiyi3001xi02yixi22xiyi03yi2xi33xiyi20−10−yi−2xi0−2xiyi−yi2−3xi20−3xi2yi−yi31xjyjxjyjxj2yj2xj2yjxjyj2xj3yj3xj3yjxjyj3001xj02yjxj22xjyj03yj2xj33xjyj20−10−yj−2xj0−2xjyj−yj2−3xj20−3xj2yj−yj31xkykxkykxk2yk2xk2ykxkyk2xk3yk3xk3ykxkyk3001xk02ykxk22xkyk03yk2xk33xkyk20−10−yk−2xk0−2xkyk−yk2−3xk20−3xk2yk−yk31xlylxlylxl2yl2xl2ylxlyk2xl3yl3xl3ylxlyl3001xl02ylxl22xlyl03yl2xl33xlyl20−10−yl−2xl0−2xlyl−yl2−3xl20−3xl2yl−yl3.


## 4.2.2. Stiffness Matrix of Quadrilateral Thin Plate Element

Since there are 12 degrees of freedom in the 4-node plate element, the stiffness matrix is a square matrix of12×12. The geometric equation of bending plate is(17)ε=εxεyγxy=Bδe,where B is a strain matrix of 3×12.(18)B=−Z∂2∂x2∂2∂y22∂2∂x∂yN.The strain energy of the plate element is(19)U=∫V12εTσdV=∫V12εTDεdV.Using the principle of minimum potential energy, we get(20)∂U∂δe=Kδe=∫BTDBz2dVδe.Therefore, the solving equation for the element stiffness matrix is(21)k=∫BTDBz2dV=∫∫h/2h/2BTDBz2dzdA,where D is the elastic matrix for materials:(22)D=E1−μ21μ0μ10001−μ2.Because the thickness of the plate element is uniform, we get(23)∫−h/2h/2z2dz=h312.Based on the above equation, we get stiffness matrix(24)k=∫BTDBz2dV=h3121−μ2∫−4.54.5dx∫−4.54.5BTDBdy.By combining all above-mentioned equations, the element stiffness matrix was solved using MATLAB programming. However, due to the complexity of solution results, they are not specifically listed hereunder. In addition, the stiffness matrix of elements of the other four sizes was solved using the same methodology, which will not be presented to keep the paper at a reasonable length.


## 4.3. Element Stiffness Matrix Assembly

Having solved stiffness matrix of the above six elements, we need to solve the overall stiffness matrix. On the entire roof, there arem elements and n nodes. Hence, the dimension of overall stiffness matrix is 3n1+2n3+1×3n2+n4+n6+2. Furthermore, element node freedom degrees are introduced to correspond to transition matrix G that has been expanded to be a structure element node. Suppose that there are a total of n structure nodes and the overall node numbers corresponding to quadrilateral thin plate elements are m, r, s, and t; then the transformation matrix for the element’s nodal freedom degrees is as follows:(25)That is, in the transformation matrix ofG, the subblocks where the overall serial numbers (m, r, s, and t) are located are set as third-order unit matrix, and all the others are zero. Using the method of transformation matrix, according to node numbering of each unit, the corresponding transformation matrix is G1,G2,G3,Gm. Since each transformation matrix is of 12×3n, the total stiffness matrix is(26)K=∑i=1mGiTkiGi,i=1,2,3,…,m.The total stiffness matrix is of3n×3n, yet, due to its huge dimension, the solution results are not elaborated upon herein.


## 4.4. Statics Solution Equation

When the transverse distributed load that the element exerts on the middle surface perpendicularly isq, the unit load vector could be calculated using the following formula:(27)Pe=∫∫ΩNTqdxdy.In this plate element, the transverse distributed load is as follows:q1=q, q2=q−kcwc, and q3=q−k1+k2w+k3w2+k4w3+k5w4+⋯. When q1=q, since q is a constant, the specific outcomes for Pe are presented in formula (27), where the unit length and width are assumed, respectively, as 2a and 2b.(28)Pe=4qab14b12−a1214b12a1214−b12a1214−b12−a12T.Besides, with formula (21), we get q3=q−k1+k2w+k3w2+k4w3+k5w4+⋯ and q2=q−kcw. After obtaining node load of the unit, the overall roof node load array can also be obtained by transition matrix(29)P=∑i=1mGiTPe,i=1,2,3,…,m.The roof is dotted with a total ofn nodes; therefore, the overall nodal force matrix P is of 3n×1bit, while the overall node displacement array δ is of 3n×1bit. But, due to their huge dimension, their expressions are not being demonstrated herein. Therefore, the total statics solving equation is(30)P3n×1=K3n×3nδ3n×1.


## 4.5. Introduction of Boundary Conditions and Equation Solutions

Since the boundaries of the roof are assumed to be fixed, displacement and rotation at the boundaries are 0. On basis of these conditions, the overall stiffness matrix can be modified as described below. The number of nodes on the boundary isp=2n1+2n2+4n3+2n4+2n6−2. If the three rows and three columns at each corresponding node in the total stiffness matrix are removed, then the new stiffness matrix would be K3n−p×3n−p′ and the new array of node displacement would be δ3n×1′. Now let us bring in force boundary conditions and the nodal force at the other nodes is Pi; then we remove a number of p nodes off the boundary and get new nodal force of P3n−p×1′. Consequently, we get a new statics balance equation as follows:(31)P3n−p×1′=K3n−p×13n−p′δ3n−p×1′,the specific solution equation of which is(32)δ3n−p×1′=K3n−p×3n−p′−1P3n−p×1′.Programming calculation for the above process is accomplished using Femap and Nastran. Specifically, Newton iterative method is applied to solve the nonlinear equations herein because this method is characteristic of simple programming and quick convergence rate; what is more, its precision can generally meet the requirements of this calculation.Firstly, we transform nonlinear equations in formula (30) into formula (32) to solve ω1, θx1, θy1, …, ω(3n-p), ωx(3n-p), θy(3n-p). Assume a certain spot as z and we could obtain the linear terms and formula (33) according to Taylor expansion of multivariate function at the point of Z(k).(33)f1ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f2ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f3ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,⋯,f3n−p−1ω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0,f3n−pω1,θx1,θy1,…,ω3n−p,θx3n−p,θy3n−p=0.(34)f1kω1k,θx1k,θy1k,…ω3n−pk,θx3n−pk,θy3n−pkf2kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf3kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk⋯fn−1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkfnkω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk+f′zkω1k+1−ωkθx1k+1−θx2kθy1k+1−θy2k⋯θx3n−p−1k+1−θxn−p−1kθy3n−pk+1−θy3n−pk=0,where(35)f′z=∂f1∂ω1⋯∂f1∂θy3n−p⋯∂f3n−p∂ω1⋯∂f3n−p∂θy3n−p.We get the iterative formula as follows:(36)ω1k+1θx1k+1⋯θ3n−p−1xk+1θ3n−pyk+1=ω1kθx1k⋯θ3n−p−1xkθ3n−pyk−f′zk−1⋅f1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf2kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkf3kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk⋯fn−1kω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pkfnkω1k,θx1k,θy1k,…,ω3n−pk,θx3n−pk,θy3n−pk.During iteration process, the function equation and partial derivatives matrix function are programmed in advance using FORTRAN language. When it comes to calculation process, the above formula would be applied constantly to do iterative computation. The results obtained are node displacement, which are not going to be listed herein due to their complexity. Then regenerate the general node displacement arrayδ3n×1 and, according to the initial static balance equation, we get nodal force at each node. Then we would acquire strain according to the geometric equations and stress according to the physics equation.We getMx, My, σx, andσy in corresponding model by solving formula (36) as follows:(37)M=MxMyMxy=∫−h/2h/2σz2dz=−h312D∂2w∂2x2∂2w∂2y22∂2w∂x∂y=−Eh3121−μ21μ0μ10001−μ/2∂2w∂2x2∂2w∂2y22∂2w∂x∂y,σ=σxσy=12M⋅zh2=12Mx⋅zh212My⋅zh2.


## 5. Engineering Example Analysis

### 5.1. Project Profile

The field trial working face is the first working face in No. 13080 backfilling and mining area of No. 12 Pingdingshan coal mine. The working face is 100 meters long with a continuous advancing length of 350 meters. The coal seam is 3.3 meters thick and its dip angle is 0∼4°. The main roof, 31.5 meters thick, is composed of fine-grained sandstone with an elasticity modulus of 28.0 GPa and Poisson’s ratio of 0.25; and the uniformly distributed load on the overlying strata is calculated to be 4.0 MPa. Description concerning the roof and floor of coal seam is shown in Figure8. The working face width, namely, the support length, is 8.0 m. Roof-contacted backfilling could be accomplished using No. ZZC8800/20/38 compactor (the interspace left for backfilling is 0). hq is 0 in this case.Figure 8
Coal seam roof descriptions in the mining area.


### 5.2. Elastic Foundation Coefficient of Backfilling Materials

Coal gangue, which is the only viable material in this coal mine, was chosen as the backfilling material in the field trial due to its easy availability and low cost. Compaction test was carried out on coal gangue samples with different particle size range and their respective stress-strain curves were plotted in Figure9. The elastic foundation coefficient corresponding to each sample is calculated and shown in Table 1. As depth of this working face is 350 m, the max stress is set to 9.0 MPa.Figure 9
Stress-strain curve of backfilling materials.Table 1
Elastic foundation coefficient of the backfilling materials in different gradation.
Particle size distribution (mm)Elastic foundation coefficientHomogeneous ratio27296.0ε5 − 12160.0ε4 + 8524.6ε3 − 1961.5ε2 + 168.0ε0∼5 mm13060.0ε5 − 17588.0ε4 + 11657.0ε3 − 3228.6ε2 + 241.5ε5∼10 mm16170.0ε5 − 5391.0ε4 + 4187.0ε3 − 1183.2ε2 + 120.7ε10∼15 mm6365.9ε5 − 2460.2ε4 + 1970.8ε3 − 575.4ε2 + 76.8ε15∼20 mm5801.7ε5 − 2255.1ε4 + 1665.0ε3 − 504.5ε2+69.8ε


### 5.3. Influencing Factors Analysis

#### 5.3.1. Influence of Particle Size Distribution

The particle size of backfilling material can evidently influence this change of roof subsidence and roof tensile stress, as shown in Figure10. The results show that the sample with homogeneous ratio produces the minimum roof tensile stress and roof subsidence volume, which indicates the best roof controlling ability.Figure 10
Roof deformation and stress under different particle size condition: roof tensile stress (a); roof subsidence volume (b).
(a)(b)


#### 5.3.2. Influence of Initial Compaction Force

Considering the aforementioned experimental results, backfilling material of homogeneous ratio was chosen to analyze the effect of initial compaction force on roof subsidence and its tensile stress. As illustrated in Figure11, backfilling material is a nonlinear elastic material. Specifically, backfilling material deforms in a strain-hardening manner during compression. Figure 11 illustrates the according roof tensile stress and roof subsidence after the initial compaction force. It is obvious that as compaction force enhances, both the roof tensile stress and its subsidence decrease but at a decreasing rate. Since tensile strength of roof strata is 7.6 MPa, when the initial compaction force exceeds 0.8 MPa, roof fracture does not occur. In industrial practice, a larger pressure means a greater cost per request, so the initial compaction force is limited under 2.0 MPa, and hence, this value is set to 2.0 MPa.Figure 11
Roof deformation and stress under different initial compaction condition: roof tensile stress (a); roof subsidence (b).
(a)(b)


### 5.4. Field Measurement

As shown in Figure12, given equipment manufacturing capacity and roof control requirements, a 2.0 MPa of initial compaction force was used. To monitor roof subsidence in backfilled area, roof subsidence meters were installed in the gob. Figure 13 displays the roof subsidence curves based on data collected from the measuring point which was located 30 meters away from the open-off cut position. When working face was advanced to 100 m, the subsidence measured by measuring points at 30 m and 50 m position was, respectively, 0.329 m and 0.363 m, which is close to the theoretical calculation results of 0.310 m and 0.350 m with errors of 6.12% and 8.86%, respectively. Therefore, it can be concluded that the proposed theoretical model is reliable and could be used for roof convergence prediction in working face with solid backfilling mining.Figure 12
Site measuring equipment and method of roof subsidence volume.Figure 13
Roof subsidence measured curve in the goaf.


## 5.1. Project Profile

The field trial working face is the first working face in No. 13080 backfilling and mining area of No. 12 Pingdingshan coal mine. The working face is 100 meters long with a continuous advancing length of 350 meters. The coal seam is 3.3 meters thick and its dip angle is 0∼4°. The main roof, 31.5 meters thick, is composed of fine-grained sandstone with an elasticity modulus of 28.0 GPa and Poisson’s ratio of 0.25; and the uniformly distributed load on the overlying strata is calculated to be 4.0 MPa. Description concerning the roof and floor of coal seam is shown in Figure8. The working face width, namely, the support length, is 8.0 m. Roof-contacted backfilling could be accomplished using No. ZZC8800/20/38 compactor (the interspace left for backfilling is 0). hq is 0 in this case.Figure 8
Coal seam roof descriptions in the mining area.


## 5.2. Elastic Foundation Coefficient of Backfilling Materials

Coal gangue, which is the only viable material in this coal mine, was chosen as the backfilling material in the field trial due to its easy availability and low cost. Compaction test was carried out on coal gangue samples with different particle size range and their respective stress-strain curves were plotted in Figure9. The elastic foundation coefficient corresponding to each sample is calculated and shown in Table 1. As depth of this working face is 350 m, the max stress is set to 9.0 MPa.Figure 9
Stress-strain curve of backfilling materials.Table 1
Elastic foundation coefficient of the backfilling materials in different gradation.
Particle size distribution (mm)Elastic foundation coefficientHomogeneous ratio27296.0ε5 − 12160.0ε4 + 8524.6ε3 − 1961.5ε2 + 168.0ε0∼5 mm13060.0ε5 − 17588.0ε4 + 11657.0ε3 − 3228.6ε2 + 241.5ε5∼10 mm16170.0ε5 − 5391.0ε4 + 4187.0ε3 − 1183.2ε2 + 120.7ε10∼15 mm6365.9ε5 − 2460.2ε4 + 1970.8ε3 − 575.4ε2 + 76.8ε15∼20 mm5801.7ε5 − 2255.1ε4 + 1665.0ε3 − 504.5ε2+69.8ε


## 5.3. Influencing Factors Analysis

### 5.3.1. Influence of Particle Size Distribution

The particle size of backfilling material can evidently influence this change of roof subsidence and roof tensile stress, as shown in Figure10. The results show that the sample with homogeneous ratio produces the minimum roof tensile stress and roof subsidence volume, which indicates the best roof controlling ability.Figure 10
Roof deformation and stress under different particle size condition: roof tensile stress (a); roof subsidence volume (b).
(a)(b)


### 5.3.2. Influence of Initial Compaction Force

Considering the aforementioned experimental results, backfilling material of homogeneous ratio was chosen to analyze the effect of initial compaction force on roof subsidence and its tensile stress. As illustrated in Figure11, backfilling material is a nonlinear elastic material. Specifically, backfilling material deforms in a strain-hardening manner during compression. Figure 11 illustrates the according roof tensile stress and roof subsidence after the initial compaction force. It is obvious that as compaction force enhances, both the roof tensile stress and its subsidence decrease but at a decreasing rate. Since tensile strength of roof strata is 7.6 MPa, when the initial compaction force exceeds 0.8 MPa, roof fracture does not occur. In industrial practice, a larger pressure means a greater cost per request, so the initial compaction force is limited under 2.0 MPa, and hence, this value is set to 2.0 MPa.Figure 11
Roof deformation and stress under different initial compaction condition: roof tensile stress (a); roof subsidence (b).
(a)(b)


## 5.3.1. Influence of Particle Size Distribution

The particle size of backfilling material can evidently influence this change of roof subsidence and roof tensile stress, as shown in Figure10. The results show that the sample with homogeneous ratio produces the minimum roof tensile stress and roof subsidence volume, which indicates the best roof controlling ability.Figure 10
Roof deformation and stress under different particle size condition: roof tensile stress (a); roof subsidence volume (b).
(a)(b)


## 5.3.2. Influence of Initial Compaction Force

Considering the aforementioned experimental results, backfilling material of homogeneous ratio was chosen to analyze the effect of initial compaction force on roof subsidence and its tensile stress. As illustrated in Figure11, backfilling material is a nonlinear elastic material. Specifically, backfilling material deforms in a strain-hardening manner during compression. Figure 11 illustrates the according roof tensile stress and roof subsidence after the initial compaction force. It is obvious that as compaction force enhances, both the roof tensile stress and its subsidence decrease but at a decreasing rate. Since tensile strength of roof strata is 7.6 MPa, when the initial compaction force exceeds 0.8 MPa, roof fracture does not occur. In industrial practice, a larger pressure means a greater cost per request, so the initial compaction force is limited under 2.0 MPa, and hence, this value is set to 2.0 MPa.Figure 11
Roof deformation and stress under different initial compaction condition: roof tensile stress (a); roof subsidence (b).
(a)(b)


## 5.4. Field Measurement

As shown in Figure12, given equipment manufacturing capacity and roof control requirements, a 2.0 MPa of initial compaction force was used. To monitor roof subsidence in backfilled area, roof subsidence meters were installed in the gob. Figure 13 displays the roof subsidence curves based on data collected from the measuring point which was located 30 meters away from the open-off cut position. When working face was advanced to 100 m, the subsidence measured by measuring points at 30 m and 50 m position was, respectively, 0.329 m and 0.363 m, which is close to the theoretical calculation results of 0.310 m and 0.350 m with errors of 6.12% and 8.86%, respectively. Therefore, it can be concluded that the proposed theoretical model is reliable and could be used for roof convergence prediction in working face with solid backfilling mining.Figure 12
Site measuring equipment and method of roof subsidence volume.Figure 13
Roof subsidence measured curve in the goaf.


## 6. Conclusions

(1)
Confined compaction test on solid backfilling materials with various particle size distribution is conducted and multiple-parameter elastic foundation coefficient for each sample is deduced.(2)
Based on coordinating relations analysis between coal body compression, backfilling materials’ deformation, and roof subsidence, a multiparameter nonlinear elastic foundation plate model using semianalytical and seminumerical method is developed. Industrial filed measurement demonstrates that the roof subsidence values from the proposed theoretical model are very close to the measured values, with an error of merely 3.24%.(3)
In combination with specific mining geological conditions, the influence of particle size distribution of material and initial compaction force on overlying strata subsidence and stress is investigated and results show the following: (i) gangues with uniform ratio have the best roof controlling effect; (ii) roof subsidence and its tensile stress decrease as initial compaction force increases; (iii) when the initial compaction has increased significantly to a certain extent, the inhibition effect that the backfilling body has for roof convergence would be weakened.



---

*Source: 4829761-2021-09-20.xml*