# Method to Calculate Working Surface Abutment Pressure Based on Key Strata Theory

Authors: Hongkai Han; Jialin Xu; Xiaozhen Wang; Jianlin Xie; Yantuan Xing
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Overburden key strata (KS) have a significant influence on abutment pressure distribution. However, current calculation methods for working surface abutment pressure do not consider the influence of the overburden KS. This study uses KS theory to analyze the overburden load transferred to coal-rock masses on both sides of the stope through fractured blocks in different layers of the KS in the fissure zone and KS in different layers of the curve subsidence zone. Using Winkler’s elastic foundation beam theory, we consider the fissure zone KS on the coal mass side and the curve subsidence zone KS as many elastic foundation beams interact with each other. A method to calculate the abutment pressure of the coal mass and the goaf was then established, considering the influence of the overburden KS. The abutment pressure distribution of working surface 207 after mining was then calculated using our method, based on mining conditions present in the Tingnan coal mine. The calculated results were verified using measurements from borehole stress meters and microseismic monitoring systems, as well as numerical simulations. In addition, the calculation results were used to determine a reasonable position for the stopping line and remaining width of the roadway’s protection coal pillar in working surface 207. The results of this study can be used to calculate the abutment pressure distribution of the working surface under a variety of overburden KS conditions. The results can also provide guidance for forecasting and preventing mine dynamic hazards, controlling the surrounding rock in mining roadways, determining reasonable widths for protection coal pillars, and designing the layout of mining roadways.


---

## Body

## 1. Introduction

Primary rock stress is normally in equilibrium before coal seams are excavated. During the course of mining, the overlying strata are shifted and stress is redistributed to the surrounding rock, resulting in abutment pressure. Failure to account for abutment pressure distribution leads to coal mine hazards, such as rockbursts, coal and gas outbursts, and roadway deformation and instability [1–9]. Forecasting and preventing these hazards, controlling the surrounding rock in mining roadways, determining a reasonable width for protection coal pillars, and designing the layout of mining roadways are all based on the distribution of abutment pressure. Empirical data from previous studies [7, 10] present a clear correlation between abutment pressure and shifts in overlying strata, and the formation of abutment pressure is closely related to overburden movement caused by coal seam mining. The overlying strata consist of dozens of rock layers, each with a different role in this movement. In China, the key strata (KS) theory is widely used in the study of rock strata movements [4, 7, 9, 11–17]. The theory states that when multiple rock layers are present in the overburden, the rock layer that controls all or part of the rock mass movements is called the KS [4, 7, 9, 11–17]. Determining the KS is mainly based on rock deformation and fracturing characteristics [7, 11–13]. Previous studies [7, 11–13] have proposed a method for determining the overlying KS. The movement of the KS results in the overall movement of all or a considerable part of the overburden. The change in the stress field caused by the movement of rock after mining is mainly controlled by the movement of the KS [7, 11–13].Previous studies have found that abutment pressure distribution is closely related to the overburden KS (number, thickness, strength, and location of KS) [4, 7, 9, 11, 12, 14, 18]. For example, calculations using the Fast Lagrangian Analysis of Continua (FLAC) have shown the distribution of abutment pressure under 100 m thick igneous KS during the advancement of a working face, thereby discovering that the advancement distance over which the abutment pressure becomes stable is far greater under the 100 m thick igneous KS than under normal geological conditions [18]. The distribution of abutment pressure under a 140 m thick igneous sill KS has also been estimated using the Universal Distinct Element Code (UDEC), which revealed that abutment pressure width is much larger under the 140 m thick igneous sill than when no igneous sill is present [4]. Xie and Xu [14] used a UDEC simulation and found that different KS locations have a significant influence on the abutment pressure distribution, and that when the KS is located further away from the coal seam, the peak value of the abutment pressure is lower while the abutment pressure width is larger. However, the current calculation methods used for determining working surface abutment pressure do not consider the influence of the overburden KS. Abutment pressure in a goaf has been shown to be positively proportional to the distance to the coal wall, which suggests that abutment pressure will recover to the initial rock stress when the distance is equal to 30% of the seam mining depth [19]. However this is not constant under different overburden KS conditions. Shi [20] considered the main roof above the coal seam as an elastic foundation beam and established a calculation method of the abutment pressure caused by the main roof; however, this method did not consider the impact of the KS above the main roof on abutment pressure. Zhu et al. [9] reported that 593 m thick alluvial stratum in the overburden had a significant influence on the distribution of abutment pressure and established an estimation method for abutment pressure under those conditions. Their method, however, is not suitable for conditions where the overburden has a thick and hard KS [9]. Theoretical estimation formulas for correlating the abutment pressure to the seam burial depth, mining width, and mining height have also been derived [21–25], but these calculation methods do not consider the influence of overburden KS on abutment pressure. For example, equation d0 = 5.13⋅h was developed by Peng and Chiang [21], where d0 is the abutment pressure width and h is the cover depth. In situ measurements of the abutment pressure from two mines in the Western United States, in conjunction with statistical analyses of the abutment pressure reported in a number of studies, have demonstrated [26] that the actual value varies greatly from results calculated using Peng and Chiang’s method [21]. Abutment pressure has been suggested to be closely associated with the composition of overburden strata [26], but an abutment pressure calculation method based on such data has not been established. Furthermore, many studies have extensively calculated the distribution of abutment pressure using a limit equilibrium [7, 27–29] which does not account for relationships between the mining width, overburden KS, and abutment pressure. Existing calculation methods also neglect the influence of the overburden KS on abutment pressure, but the characteristics of the overburden KS inevitably impact the size of the overlying load that shifts toward either side of the coal-rock mass in a coal mine, thereby affecting abutment pressure distribution. Therefore, in this study, we attempt to establish a method to calculate abutment pressure, taking into consideration the influence of the overburden KS.This study uses KS theory [7, 11] and Winkler’s theory of elastic foundation beams [12, 15, 16, 20, 30] to establish a method to calculate abutment pressure that considers the influence of the overburden KS. The abutment pressure distribution of working surface 207 after mining was then calculated using our method, based on mining conditions present in the Tingnan coal mine. These results were verified by measurements from borehole stress meters and microseismic monitoring systems, as well as by numerical simulation results. In addition, the calculation results were used to determine a reasonable position for the stopping line and the remaining width of the roadway’s protection coal pillar for working surface 207.


## 2. Case Analysis of Abutment Pressure Distribution Affected by Overburden KS

Neglecting the influence of overburden KS on the distribution of abutment pressure often leads to catastrophic accidents, as described in the two following cases.


### 2.1. Deformation and Damage to Roadways due to Abutment Pressure

The overburden of the thick sandstone KS on abutment pressure distribution was not considered during the operation of the Tingnan coal mine, and as a result, the width of the protective coal pillars in the main roadway was insufficient. As a result, the main roadway 200 m from the stopping line underwent significant deformation. The Tingnan coal mine is located in the Changwu County, Xianyang, Shaanxi Province, and it has an average mining height of 7.0 m and a coal seam dip angle of 0° in the second panel. Figure1 shows the layout of the working area in the second panel. There is an ultrathick (222.14 m), hard sandstone KS, which is a moderately water-rich aquifer with an average uniaxial compressive strength of 42.2 MPa. On the basis of the integrated stratigraphic column in the Tingnan coal mine [30], the location of the overburden KS is shown in Figure 2. The recovery periods of working faces 204, 205, and 206 were Nov. 2011–Nov. 2012, Jun. 2013–Nov. 2014, and Jan. 2015–Apr. 2016, respectively. Currently, the panel is recovering working face 207. The mine design left a 200 m-wide protection coal pillar between the panel working surfaces stopping line and the main roadway, outside of the experience width of the abutment pressure [7, 21]. There was significant deformation and damage to the main roadway caused by the working surfaces after extraction, as shown in Figure 3.Figure 1
Layout sketch map of the working area in the second panel.Figure 2
Location of the overburden KS in the Tingnan coal mine.Figure 3
Photos of the main roadway damage.
(a)
(b)


### 2.2. Coal and Gas Outburst Caused by Abutment Pressure

The Haizi coal mine is located in Suixi County, Huaibei City, Anhui Province, China. The average coal seam thickness in the II102 mining area is 2.5 m, the average depth is 623 m, and a 140 m thick igneous KS is located approximately 170 m from the coal seam. After the extraction of working surfaces II1022 and II1024, a coal and gas outburst occurred 180 m from the goaf while a transport roadway was being excavated in working surface II1026, as shown in Figure4. Again, this outburst occurred outside of the experience width of the abutment pressure [7, 21] because the influence of the overburden of thick igneous KS on the abutment pressure distribution was neglected [4].Figure 4
Layout sketch map of the working area in the Haizi coal mine.The deformation and damage to the main roadway in the Tingnan coal mine and the coal and gas outburst in the Haizi coal mine [4] were closely related to the distribution of abutment pressure. Previous studies have determined empirical values [7, 21] regarding the abutment pressure width. Qian et al. [7] determined that the abutment pressure width is generally 15–40 m. According to Peng and Chiang’s [21] calculation method, the abutment pressure widths of the Tingnan and Haizi coal mines were 118 m and 128 m, respectively. According to these results, the roadways 200 m from the stopping line in the Tingnan coal mine and 180 m from the goaf in the Haizi coal mine should have been outside of the width of abutment pressure and under primary rock stress. This does not explain the deformation and damage to the roadway in the Tingnan coal mine nor the coal and gas outburst in the Haizi coal mine, both of which are due to neglecting the presence of thick, hard sandstone KS and igneous rock KS in the overburdens of the Tingnan and Haizi coal mines, respectively. The KS causes the actual width of abutment pressure to be significantly larger than the experience width [7, 21], resulting in the above phenomena. Therefore, the significant impact of overburden KS on abutment pressure distribution must not be ignored, and a calculation method for abutment pressure that accounts for this influence should be established.


## 2.1. Deformation and Damage to Roadways due to Abutment Pressure

The overburden of the thick sandstone KS on abutment pressure distribution was not considered during the operation of the Tingnan coal mine, and as a result, the width of the protective coal pillars in the main roadway was insufficient. As a result, the main roadway 200 m from the stopping line underwent significant deformation. The Tingnan coal mine is located in the Changwu County, Xianyang, Shaanxi Province, and it has an average mining height of 7.0 m and a coal seam dip angle of 0° in the second panel. Figure1 shows the layout of the working area in the second panel. There is an ultrathick (222.14 m), hard sandstone KS, which is a moderately water-rich aquifer with an average uniaxial compressive strength of 42.2 MPa. On the basis of the integrated stratigraphic column in the Tingnan coal mine [30], the location of the overburden KS is shown in Figure 2. The recovery periods of working faces 204, 205, and 206 were Nov. 2011–Nov. 2012, Jun. 2013–Nov. 2014, and Jan. 2015–Apr. 2016, respectively. Currently, the panel is recovering working face 207. The mine design left a 200 m-wide protection coal pillar between the panel working surfaces stopping line and the main roadway, outside of the experience width of the abutment pressure [7, 21]. There was significant deformation and damage to the main roadway caused by the working surfaces after extraction, as shown in Figure 3.Figure 1
Layout sketch map of the working area in the second panel.Figure 2
Location of the overburden KS in the Tingnan coal mine.Figure 3
Photos of the main roadway damage.
(a)
(b)


## 2.2. Coal and Gas Outburst Caused by Abutment Pressure

The Haizi coal mine is located in Suixi County, Huaibei City, Anhui Province, China. The average coal seam thickness in the II102 mining area is 2.5 m, the average depth is 623 m, and a 140 m thick igneous KS is located approximately 170 m from the coal seam. After the extraction of working surfaces II1022 and II1024, a coal and gas outburst occurred 180 m from the goaf while a transport roadway was being excavated in working surface II1026, as shown in Figure4. Again, this outburst occurred outside of the experience width of the abutment pressure [7, 21] because the influence of the overburden of thick igneous KS on the abutment pressure distribution was neglected [4].Figure 4
Layout sketch map of the working area in the Haizi coal mine.The deformation and damage to the main roadway in the Tingnan coal mine and the coal and gas outburst in the Haizi coal mine [4] were closely related to the distribution of abutment pressure. Previous studies have determined empirical values [7, 21] regarding the abutment pressure width. Qian et al. [7] determined that the abutment pressure width is generally 15–40 m. According to Peng and Chiang’s [21] calculation method, the abutment pressure widths of the Tingnan and Haizi coal mines were 118 m and 128 m, respectively. According to these results, the roadways 200 m from the stopping line in the Tingnan coal mine and 180 m from the goaf in the Haizi coal mine should have been outside of the width of abutment pressure and under primary rock stress. This does not explain the deformation and damage to the roadway in the Tingnan coal mine nor the coal and gas outburst in the Haizi coal mine, both of which are due to neglecting the presence of thick, hard sandstone KS and igneous rock KS in the overburdens of the Tingnan and Haizi coal mines, respectively. The KS causes the actual width of abutment pressure to be significantly larger than the experience width [7, 21], resulting in the above phenomena. Therefore, the significant impact of overburden KS on abutment pressure distribution must not be ignored, and a calculation method for abutment pressure that accounts for this influence should be established.


## 3. Calculation Method for Abutment Pressure

### 3.1. Mechanical Model for the Calculation of Abutment Pressure

In Chinese coal mines, the most commonly used mining method is longwall mining in which the stoping space created by the stoping of one or more working faces is rectangular in shape and the mining area contains two main sections [30, 31]: the along-dip main section I and the along-strike main section II, as shown in Figure 1. As the working face is advanced along main section II, often for thousands of meters, this section is deemed to be under full subsidence [30, 31]. A mechanical model for calculating abutment pressure that considers the overburden KS was established along main section I, as shown in Figure 5. After a coal seam is mined, a breaking movement and a flexural subsidence are expected in the overburden strata of the mining site, which is composed of the caving zone, fissure zone, and curve subsidence zone [11, 12, 17, 21, 30]. The KS, which is in the curve subsidence zone, is subject to a curve subsidence rather than a breaking movement [11, 12, 17, 30]. The height of the fissure zone can be estimated by applying Xu et al.’s prediction method [17] based on the KS position. Since the caving zone rock mass does not exert a horizontal force onto the coal-rock mass around the stope [7, 11, 12], it does not transfer its load to the coal-rock masses on either side of the stope. KS theory demonstrates that [7, 9, 11, 12, 15, 16] after the KS is fractured, the load is shifted to both sides of the stope through the fractured blocks of the KS. When the KS is not fractured, it transfers its load to both sides of the stope. Thus, part of the load of the stope’s overburden will be exerted onto the coal-rock masses on both sides of the stope through the fractured blocks in different KS of the fissure zone (e.g., the fractured blocks A1, A2, …, Am−1, as shown in Figure 5) as well as the KS in different layers of the curve subsidence zone. The other part of the load is transferred to the goaf. We assumed that the solid coal-rock masses and the rock masses in the curve subsidence zone, caving zone, and fissure zone satisfy Winkler’s elastic foundation assumption [15, 16, 20, 30, 32, 33]. Using Winkler’s elastic foundation beam theory, we consider the fissure zone KS on the coal mass side and curve subsidence zone KS as many elastic foundation beams interact with each other. For any two adjacent KS, the load of the upper KS in the adjacent KS is transmitted to the lower KS by a nonuniform load through the elastic foundation consisting of the lower KS and its rock load, while the weight of the lower KS and its load acts on the lower KS in a uniform manner. The load is transferred layer by layer through the KS and ultimately forms abutment pressure in the coal mass and in the goaf.Figure 5
Mechanical model for abutment pressure calculation.In Figure5, the lowermost KS in the fissure zone is KS1 and the uppermost KS in that zone is KSm−1. The lowermost KS in the curve subsidence zone is KSm, and the uppermost KS in that zone is KSn. The weight of different locations of KS and the rock loads they control is represented by q1, q2, …, qm−1, qm, …, qn−1, qn. For example, qm refers to the weight of the interlayer rock between KSm+1 and KSm and that of KSm and is measured in MPa. The weight of the rock layer in the caving zone, measured in MPa, is represented by q0. In adjacent KS, the foundation coefficients of the elastic foundation consisting of the lower KS and its rock load are represented by k1, k2, …, km−1, km, …, kn−1, kn. For example, km refers to the foundation coefficient of the foundation formed by the interlayer rock between KSm and KSm−1 as well as KSm−1 and is measured in N/m3. The foundation coefficient of the KSm on the goaf side is kcm and is measured in N/m3. The elastic modulus of different layers of KS, E1, E2, …, Em−1, Em, …, En−1, En, is measured in GPa. The moment of inertia of different layers of KS, I1, I2, …, Im−1, Im, …, In−1, In, is measured in m4. The thickness of different layers of KS, H1, H2, …, Hm−1, Hm, …, Hn−1, Hn, are measured in m. The mining width is L and is measured in m.As the KS in Figure5 are symmetric about x = L/2, only half of the KS were used for our analysis. The stress analysis of the fissure zone KS on the coal mass side and curve subsidence zone KS is as follows.


#### 3.1.1. Stress Analysis of the Fissure Zone KS on the Coal Mass Side

Using KSm−1 in Figure 5 as an example, the corresponding stress analysis is shown in Figure 6. In this study, Q1(x), Q2(x), …, Qm−1(x), Qm(x), …, Qn−1(x), Qn(x) and Qcm(x), …, Qcn−1(x), Qcn(x) represent the shearing forces experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. For example, in Figure 6, Qm−1(0) is the shearing force borne by the side of the coal mass KSm−1 at x = 0. The shearing force of the different layers of KS in the fissure zone at x = 0 is approximately equal to half of the product of the lengths of the fractured KS blocks and their load in the goaf [9]. The lengths of the KS fractured blocks A1, A2, …, Am−1 are represented by l1, l2, …, lm−1, respectively, and can be calculated using the KS discriminant method [7, 11, 12]. If the number of KS layers in the fissure zone is small, the lengths of the KS fractured blocks can be determined by a measured pressure curve analysis of the hydraulic support of the mine’s working surface [34]. Hence, Qm−1(0) in Figure 6 is(1)Qm−10=lm−1qm−12.Figure 6
Stress analysis of the fissure zone KS on the coal mass side.y
1(x), y2(x), …, ym−1(x), ym(x), …, yn−1(x), yn(x) and ycm(x), …, ycn−1(x), ycn(x) are the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf, respectively, and are measured in m. For example, ym−1(x) is the deflection curve equation of KSm−1 on the coal mass side. p1(x), p2(x), …, pm−1(x), pm(x), …, pn−1(x), pn(x) and pcm(x), …, pcn−1(x), pcn(x) represent the force on the foundation exerted by the different layers of KS on the side of the coal mass and the side of the goaf, which is also the force exerted by KS in different locations on the adjacent KS below them. This is measured in MPa. As shown in Figure 6, pm(x) is the force of KSm on the foundation and also the load of KSm acting on KSm−1. Hence, pm−1(x) and pm(x) in Figure 6 are calculated as follows:(2)pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x.


#### 3.1.2. Stress Analysis of the Curve Subsidence Zone KS

Using KSn−1 in Figure 5 as an example, the stress analysis is shown in Figures 7(a) and 7(b). M1(x), M2(x), …, Mm−1(x), Mm(x), …, Mn−1(x), Mn(x), and Mcm(x), …, Mcn−1(x), Mcn(x) represent the bending moments borne by different layers of KS on the side of the coal mass and the side of the goaf, respectively. As shown in Figure 7, Mn−1(x) and Mcn−1(x) are the bending moments of KSn−1 on the side of the coal mass and side of the goaf, respectively. Qn−1(0) and Qcn−1(0) are the shearing forces of KSn−1 on the side of the coal mass and the side of the goaf at x = 0, respectively.Figure 7
Stress analysis of the curve subsidence zone KS: (a) side of the coal mass; (b) side of the goaf.
(a)
(b)p
n−1(x) and pn(x) in Figure 7(a) are as follows:(3)pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x.p
cn−1(x) and pcn(x) in Figure 7(b) are as follows:(4)pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x.The relationship between the elastic foundation beam’s deflection, foundation pressure, and the load it bears is as follows [15, 16, 20, 30, 32, 33]:(5)EId4yxdx4+px=qx,where E is the elastic modulus of the beam, I is its moment of inertia, y(x) is its deflection equation, p(x) is the foundation pressure, and q(x) is the load that the beam bears.According to the established mechanical model, the following groups of deflection curve differential equations can be formed:(6)E1I1d4y1xdx4+bk1y1x=bq1+bk2y2x−y1x,E2I2d4y2xdx4+bk2y2x−y1x=bq2+bk3y3x−y2x,⋮Em−1Im−1d4ym−1xdx4+bkm−1ym−1x−ym−2x=bqm−1+bkmymx−ym−1x,EmImd4ymxdx4+bkmymx−ym−1x=bqm+bkm+1ym+1x−ymx,⋮En−1In−1d4yn−1xdx4+bkn−1yn−1x−yn−2x=bqn−1+bknynx−yn−1x,EnInd4ynxdx4+bknynx−yn−1x=bqn,x≤0,(7)EmImd4ycmxdx4+bkcmycmx=bqm+bkm+1ycm+1x−ycmx,⋮En−1In−1d4ycn−1xdx4+bkn−1ycn−1x−ycn−2x=bqn−1+bknycnx−ycn−1x,EnInd4ycnxdx4+bknycnx−ycn−1x=bqn,x≥0,where b(m) is the width of the beam, considered to be equal to unity, i.e., b = 1.Equation group (6) is the deflection curve differential equation for different KS on the side of the coal mass. Equation group (7) is the deflection curve differential equation of different KS on the side of the goaf. Solving the equations in equation group (6) gives(8)and4ny1xdx4n+an−1d4n−1y1xdx4n−1+⋯+a2d8y1xdx8+a1d4y1xdx4+an+1y1x=an+2.As the expressionsa1, a2, …, an, an+1, an+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; k1, k2, …, kn−1, kn; and q1, q2, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (6).The characteristic equation corresponding to equation (8) is(9)anr4n+an−1r4n−1+⋯+a2r8+a1r4+an+1=0.The root of equation (9) is r1,2,3,4 = ±t1 ± t1 · i, r5,6,7,8 = ±t2 ± t2 · i,  …,  r4(n−1)−3,4(n−1)−2,4(n−1)−1,4(n−1) = ±tn−1 ± tn−1 · i, and r4n−3,4n−2,4n−1,4n = ±tn ± tn · i. As the expressions t1, t2, …, tn−1, tn are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into equation (9).A special solution for equation (8) is(10)y∗=q1+q2+⋯+qn−1+qnk1.Then, the general solution for equation (8) is(11)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A3e−x⋅t1cosx⋅t1+A4e−x⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+A7e−x⋅t2cosx⋅t2+A8e−x⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−1−1e−x⋅tn−1cosx⋅tn−1+A4n−1e−x⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+A4n−1e−x⋅tncosx⋅tn+A4ne−x⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1.Asx ⟶ −∞ and y1(x) ⟶ (q1 + q2 + ⋯ + qn−1 + qn)/k1. Thus, it is only satisfied when A3, A4, A7, A8, …, A4(n−1)−1, A4(n−1), A4n−1, A4n all have a value of 0, and the deflection curve equation y1(x) is obtained. Substituting this into equation group (6), the deflection curve equations of the different layers of KS on the side of the coal mass are shown in equation (12).As the expressionsB2,1, B2,2, …, B2,n−1, B2,n; …; Bn−1,1, Bn−1,2, …, Bn−1,n−1, Bn−1,n; Bn,1, Bn,2, …, Bn,n−1, Bn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into group (6). The deflection curve equation group of the different layers of KS on the side of the coal mass has 2n unknowns, A1, A2, A5, A6, …, A4(n−1)−3, A4(n−1)−2, A4n−3, A4n−2:(12)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1,y2x=B2,1A1ex⋅t1cosx⋅t1+B2,1A2ex⋅t1sinx⋅t1+B2,2A5ex⋅t2cosx⋅t2+B2,2A6ex⋅t2sinx⋅t2+⋯+B2,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+B2,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+B2,nA4n−3ex⋅tncosx⋅tn+B2,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2,⋮ym−1x=Bm−1,1A1ex⋅t1cosx⋅t1+Bm−1,1A2ex⋅t1sinx⋅t1+Bm−1,2A5ex⋅t2cosx⋅t2+Bm−1,2A6ex⋅t2sinx⋅t2+⋯+Bm−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm−1,nA4n−3ex⋅tncosx⋅tn+Bm−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1,ymx=Bm,1A1ex⋅t1cosx⋅t1+Bm,1A2ex⋅t1sinx⋅t1+Bm,2A5ex⋅t2cosx⋅t2+Bm,2A6ex⋅t2sinx⋅t2+⋯+Bm,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm,nA4n−3ex⋅tncosx⋅tn+Bm,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1+qm+⋯+qn−1+qnkm,⋮yn−1x=Bn−1,1A1ex⋅t1cosx⋅t1+Bn−1,1A2ex⋅t1sinx⋅t1+Bn−1,2A5ex⋅t2cosx⋅t2+Bn−1,2A6ex⋅t2sinx⋅t2+⋯+Bn−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn−1,nA4n−3ex⋅tncosx⋅tn+Bn−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1,ynx=Bn,1A1ex⋅t1cosx⋅t1+Bn,1A2ex⋅t1sinx⋅t1+Bn,2A5ex⋅t2cosx⋅t2+Bn,2A6ex⋅t2sinx⋅t2+⋯+Bn,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn,nA4n−3ex⋅tncosx⋅tn+Bn,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1+qnkn.Solving the equations in equation group (7) gives(13)bnd4n−m+1ycmxdx4n−m+1+bn−1d4n−mycmxdx4n−m+⋯+bm+1d8ycmxdx8+bmd4ycmxdx4+bn+1ycmx=bn+2.As the expressionsbm, bm+1, …, bn, bn+1, bn+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; kcm, km+1, …, kn−1, kn; and qm, qm+1, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (7).The characteristic equation corresponding to equation (13) is(14)bnr4n−m+1+bn−1r4n−m+⋯+bm+1r8+bmr4+bn+1=0.The root of equation (14) is r4m−3,4m−2,4m−1,4m = ±sm ± sm · i, …, r4(n−1)−3, 4(n−1)−2, 4(n−1)−1, 4(n−1) = ±sn−1 ± sn−1 · i, r4n−3,4n−2,4n−1,4n = ±sn ± sn · i. As the expressions sm, …, sn−1, sn are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into equation (14).A special solution for equation (13) is(15)y∗=qm+⋯+qn−1+qnkcm.Then, the general solution for equation (13) is(16)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,(17)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,⋮ycn−1x=Dn−1,mC4m−3ex⋅smcosx⋅sm+Dn−1,mC4m−2ex⋅smsinx⋅sm+Dn−1,mC4m−1e−x⋅sm·cosx⋅sm+Dn−1,mC4me−x⋅smsinx⋅sm+⋯+Dn−1,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn−1,n−1C4n−1−2ex⋅sn−1sinx⋅sn−1+Dn−1,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn−1,n−1·C4n−1e−x⋅sn−1sinx⋅sn−1+Dn−1,nC4n−3ex⋅sncosx⋅sn+Dn−1,nC4n−2ex⋅snsinx⋅sn+Dn−1,nC4n−1e−x⋅sncosx⋅sn+Dn−1,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1,ycnx=Dn,mC4m−3ex⋅smcosx⋅sm+Dn,mC4m−2ex⋅smsinx⋅sm+Dn,mC4m−1e−x⋅smcosx⋅sm+Dn,mC4me−x⋅smsinx⋅sm+⋯+Dn,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1−2·ex⋅sn−1sinx⋅sn−1+Dn,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1e−x⋅sn−1sinx⋅sn−1+Dn,nC4n−3ex⋅sncosx⋅sn+Dn,nC4n−2ex⋅snsinx⋅sn+Dn,nC4n−1e−x⋅sncosx⋅sn+Dn,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1+qnkn.Substituting equation (16) into equation group (7), the deflection curve equations of the different layers of KS on the side of the goaf are shown in equation (17). As the expressions Dm+1,m, …, Dm+1,n−1, Dm+1,n; …; Dn−1,m, …, Dn−1,n−1, Dn−1,n; Dn,m, …, D n,n−1, Dn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into group (7). The deflection curve equation group of the different layers of KS on the side of the goaf has 4(n−m+1) unknowns, C4m−3, C4m−2, C4m−1, C4m; …; C4(n−1)−3, C4(n−1)−2, C4(n−1)−1, C4(n−1);C4n−3, C4n−2, C4n−1, C4n.In the paper,θ1(x), θ2(x), …, θm−1(x), θm(x), …, θn−1(x), θn(x), and θcm(x), …, θcn−1(x), θcn(x) represent the rotation angles experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. From the relationship among the rotation angle, bending moment, and shearing force of the KS with its deflection curve equation in equations (12) and (17) [15, 30, 32], by substituting the deflection curve equations of the different layers of KS into equation (18), we obtain θ1(x), …, θn−1(x), θn(x); M1(x), …, Mn−1(x), Mn(x); Q1(x), …, Qn−1(x), Qn(x); θcm(x), …, θcn−1(x), θcn(x); Mcm(x), …, Mcn−1(x), Mcn(x); and Qcm(x), …, Qcn−1(x), Qcn(x):(18)θx=dyxdx,Mx=−EIdy2xdx2,Qx=−EIdy3xdx3,(19)ym0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qm−1+qm+⋯+qn−1+qnkm−1−qm+⋯+qn−1+qnkm=ycm0,⋮yn−10−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1=ycn−10−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1,yn0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1−qnkn=ycn0−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1−qnkn,θm0=θcm0,…,θn−10=θcn−10,θn0=θcn0,Mm0=Mcm0,…,Mn−10=Mcn−10,Mn0=Mcn0,Qm0=Qcm0,…,Qn−10=Qcn−10,Qn0=Qcn0,Q10=l1q12,Q20=l2q22,…,Qm−10=lm−1qm−12;M10=0,M20=0,…,Mm−10=0,θcmL2=0,…,θcn−1L2=0,θcnL2=0,QcmL2=0,…,Qcn−1L2=0,QcnL2=0.There are a total of 4(n − m + 1) + 2n unknowns in equations (12) and (17), which are the deflection curve equations of the KS on the side of the coal mass and on the side of the goaf, respectively. The different curve subsidence zone KS on the side of the coal mass and the side of the goaf are continuous at x = 0. At x = 0, the KS on both sides have equal subsidence, rotation angle, bending moment, and shearing force. In the fissure zone, different KS on the side of the coal mass have a bending moment of 0 at x = 0, and the shearing force at x = 0 can be determined using equation (6). The model is symmetric about x = L/2; thus, the shearing force and rotation angle of the KS at x = L/2 is 0. Considering that in equations (12) and (17) (q1 + q2 + ⋯ + qn−1 + qn)/k1, (q2 + ⋯ + qn−1 + qn)/k2, …, (qm−1 + qm + ⋯ + qn−1 + qn)/km−1, (qm + ⋯ + qn−1 + qn)/km, …, (qn−1 + qn)/kn−1, qn/kn are caused by stratal compression before the coal seam is mined, the continuity and boundary conditions of the different KS are shown in equation (19).There are a total of 2n + 4(n − m + 1) equations that can be used to obtain the unknowns in equations (12) and (17) and obtain the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf. In particular, these unknown values will no longer vary with the increase of L after the land surface has fully subsided. The foundation pressure equations of the different layers of KS on the side of the coal mass and the side of the goaf are(20)p1x=k1y1x,p2x=k2y2x−y1x,⋮pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x,⋮pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x,x≤0,pcmx=kcmycmx,⋮pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x,x≥0.


### 3.2. Method for Calculating Abutment Pressure

#### 3.2.1. Calculation Method for Coal Mass Abutment Pressure

According to Section3.1, the force exerted by the different layers of KS on the coal mass is p1(x). Considering the weight of the rock mass on the side of the caving zone, q0, the coal mass stress is(21)σ1x=p1x+q0.When the coal lateral pressure in the proximity of the mining boundary is zero, i.e., when the area is under uniaxial compression, as the plastic failure [7, 29, 35] moves inward, the coal will gradually become triaxially stressed, thereby enhancing its compressive strength [7]. If the coal strength on the coal mass side near the mining boundary is smaller than the abutment pressure value in the corresponding location, which can be calculated using (21), a plastic failure will occur in this part of the coal. In this case, the equation for the abutment pressure in the plastic failure region (equation (22)), can be obtained from Ref. [7]:(22)σpx=N0e−2fx/M1+sinφ/1−sinφ,x≤0.The symbolσp(x) represents the abutment pressure of the coal plastic region, M is the mining height of the seam, φ is the internal angle of friction between the seam and the roof/floor slate stratum, f is the friction factor between the stratum surfaces, and N0 is the supporting capacity of the coal wall.As (22) is equal to (21), the width |x1| of the coal plastic region can be estimated, along with the peak abutment pressure σ(x1) on the interface between the coal plastic region and the elastic region. The abutment pressure in the plastic region and in the elastic region can be calculated using (21) and (22), respectively. Hence, the coal mass abutment pressure is(23)σx=σ1x,x≤x1,σpx,x1≤x≤0.Ifσ(x1) = K(q0 + q1 + q2 + ⋯ + qn−1 + qn), then K is the stress concentration coefficient. If σ1(xd) = q0 + q1 + q2 + ⋯ + qn−1 + qn, then |xd| is the abutment pressure width of the side of the coal mass.


#### 3.2.2. Calculation Method for Goaf Abutment Pressure

According to Section3.1, the stress distribution of the curve subsidence zone KS in the goaf is pcm(x). By superimposing it on the stress distribution caused by the fissure zone and the caving zone fractured rock masses in the goaf, we can obtain the distribution of abutment pressure in the goaf. It is assumed that the KS fractured blocks in the fissure zone far away from the mining boundary transfer the entire weight to the goaf [9], and the KS fractured blocks at the mining boundary A1, …, Am−1 transfer half of their weight to the goaf [9, 23]. The load distribution that they form in the goaf increases linearly [19]. Thus, the stress distribution formed by the fissure zone KS in the goaf is as follows:(24)σfx=q1l1+q2l2+⋯+qm−1lm−1x,0<x≤l1,q2l2+⋯+qm−1lm−1x+q1,l1≤x≤l2,⋮⋮qm−1lm−1x+q1+q2+⋯+qm−2,lm−2≤x≤lm−1,q1+q2+⋯+qm−2+qm−1,x≥lm−1.The distribution of the abutment pressure in the goaf can be calculated as follows:(25)σcx=pcmx+σfx+q0,x>0.


### 3.3. Determining Model Parameters

E
1, E2, …, Em−1, Em, …, En−1, En and H1, H2, …, Hm−1, Hm, …, Hn−1, Hn as well as L and M can be determined using the physical mechanical parameter tests on drill core samples and the specific mining conditions.


#### 3.3.1. Determination ofI1, I2, …, Im−1, Im, …, In−1, In

I
m is used as an example to illustrate the determination of the moment of inertia for different layers of KS. Im is derived using the following equation:(26)Im=112bHm3.


#### 3.3.2. Determination ofq1, q2, …, qm−1, qm, …, qn−1, qn

q
m−1 is used as an example to illustrate the respective load for different KS and the load layer weight they control. If the adjacent KS interlayer rock layers are named z1, z2, …, zi−1, zi, with a total of i layers, as shown in Figure 8, then qm−1 is(27)qm−1=γm−1Hm−1+∑i=1iγziHzi.Figure 8
Diagram of adjacent KS.In Figure8, γz1, γz2, …, γzi−1, γzi represent the density of the rock mass z1, z2, …, zi−1, zi, respectively; Ez1, Ez2, …, Ezi−1, Ezi represent the elastic modulus of the rock mass z1, z2, …, zi−1, zi, respectively; Hz1, Hz2, …, Hzi−1, Hzi represent the thickness of the rock mass z1, z2, …, zi−1, zi, respectively. In this study, γ1, γ2, …, γm−1, γm, …, γn−1, γn represent the density of the different layers of KS.


#### 3.3.3. Determination ofk1, k2, …, km−1, km, …, kn−1, kn

k
m is used as an example to illustrate the calculation method of k1, k2, …, km−1, km, …, kn−1, kn. The foundation coefficient [15, 16, 20, 30] formed by the total of i layers and KSm−1 shown in Figure 8 is(28)km=1Hm−1/Em−1+∑i=1iHzi/Ezi.In particular, the calculation ofk1 should consider the strata from KS1 to the thick hard rock formation in the coal seam floor [30].


#### 3.3.4. Determination ofkcm

Because the model assumes that the fractured rock in the caving zone and fissure zone satisfies Winkler’s elastic foundation hypothesis [15, 16, 20, 30, 32, 33], we can calculate the foundation coefficient kcm of KSm on the goaf side if we know the subsidence value of KSm and the corresponding load bearing of the foundation [30]. When the land surface is under full subsidence [31], the level of subsidence in the curve subsidence zone is consistent [25], and the subsidence of KSm reaches the maximum subsidence [31], δM. At this point, the load of the foundation that corresponds to the maximum subsidence of the side of the goaf, KSm, returns to primary rock stress [7, 11] qm + qm+1 + ⋯ + qn−1 + qn. Hence, the foundation coefficient of the KSm on the side of the goaf, kcm, is(29)kcm=qm+qm+1+…+qn−1+qnδM,where δ represents the subsidence coefficient [30, 31], which can be found in the literatures [30, 31], and M is the mining thickness of the seam.


## 3.1. Mechanical Model for the Calculation of Abutment Pressure

In Chinese coal mines, the most commonly used mining method is longwall mining in which the stoping space created by the stoping of one or more working faces is rectangular in shape and the mining area contains two main sections [30, 31]: the along-dip main section I and the along-strike main section II, as shown in Figure 1. As the working face is advanced along main section II, often for thousands of meters, this section is deemed to be under full subsidence [30, 31]. A mechanical model for calculating abutment pressure that considers the overburden KS was established along main section I, as shown in Figure 5. After a coal seam is mined, a breaking movement and a flexural subsidence are expected in the overburden strata of the mining site, which is composed of the caving zone, fissure zone, and curve subsidence zone [11, 12, 17, 21, 30]. The KS, which is in the curve subsidence zone, is subject to a curve subsidence rather than a breaking movement [11, 12, 17, 30]. The height of the fissure zone can be estimated by applying Xu et al.’s prediction method [17] based on the KS position. Since the caving zone rock mass does not exert a horizontal force onto the coal-rock mass around the stope [7, 11, 12], it does not transfer its load to the coal-rock masses on either side of the stope. KS theory demonstrates that [7, 9, 11, 12, 15, 16] after the KS is fractured, the load is shifted to both sides of the stope through the fractured blocks of the KS. When the KS is not fractured, it transfers its load to both sides of the stope. Thus, part of the load of the stope’s overburden will be exerted onto the coal-rock masses on both sides of the stope through the fractured blocks in different KS of the fissure zone (e.g., the fractured blocks A1, A2, …, Am−1, as shown in Figure 5) as well as the KS in different layers of the curve subsidence zone. The other part of the load is transferred to the goaf. We assumed that the solid coal-rock masses and the rock masses in the curve subsidence zone, caving zone, and fissure zone satisfy Winkler’s elastic foundation assumption [15, 16, 20, 30, 32, 33]. Using Winkler’s elastic foundation beam theory, we consider the fissure zone KS on the coal mass side and curve subsidence zone KS as many elastic foundation beams interact with each other. For any two adjacent KS, the load of the upper KS in the adjacent KS is transmitted to the lower KS by a nonuniform load through the elastic foundation consisting of the lower KS and its rock load, while the weight of the lower KS and its load acts on the lower KS in a uniform manner. The load is transferred layer by layer through the KS and ultimately forms abutment pressure in the coal mass and in the goaf.Figure 5
Mechanical model for abutment pressure calculation.In Figure5, the lowermost KS in the fissure zone is KS1 and the uppermost KS in that zone is KSm−1. The lowermost KS in the curve subsidence zone is KSm, and the uppermost KS in that zone is KSn. The weight of different locations of KS and the rock loads they control is represented by q1, q2, …, qm−1, qm, …, qn−1, qn. For example, qm refers to the weight of the interlayer rock between KSm+1 and KSm and that of KSm and is measured in MPa. The weight of the rock layer in the caving zone, measured in MPa, is represented by q0. In adjacent KS, the foundation coefficients of the elastic foundation consisting of the lower KS and its rock load are represented by k1, k2, …, km−1, km, …, kn−1, kn. For example, km refers to the foundation coefficient of the foundation formed by the interlayer rock between KSm and KSm−1 as well as KSm−1 and is measured in N/m3. The foundation coefficient of the KSm on the goaf side is kcm and is measured in N/m3. The elastic modulus of different layers of KS, E1, E2, …, Em−1, Em, …, En−1, En, is measured in GPa. The moment of inertia of different layers of KS, I1, I2, …, Im−1, Im, …, In−1, In, is measured in m4. The thickness of different layers of KS, H1, H2, …, Hm−1, Hm, …, Hn−1, Hn, are measured in m. The mining width is L and is measured in m.As the KS in Figure5 are symmetric about x = L/2, only half of the KS were used for our analysis. The stress analysis of the fissure zone KS on the coal mass side and curve subsidence zone KS is as follows.


### 3.1.1. Stress Analysis of the Fissure Zone KS on the Coal Mass Side

Using KSm−1 in Figure 5 as an example, the corresponding stress analysis is shown in Figure 6. In this study, Q1(x), Q2(x), …, Qm−1(x), Qm(x), …, Qn−1(x), Qn(x) and Qcm(x), …, Qcn−1(x), Qcn(x) represent the shearing forces experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. For example, in Figure 6, Qm−1(0) is the shearing force borne by the side of the coal mass KSm−1 at x = 0. The shearing force of the different layers of KS in the fissure zone at x = 0 is approximately equal to half of the product of the lengths of the fractured KS blocks and their load in the goaf [9]. The lengths of the KS fractured blocks A1, A2, …, Am−1 are represented by l1, l2, …, lm−1, respectively, and can be calculated using the KS discriminant method [7, 11, 12]. If the number of KS layers in the fissure zone is small, the lengths of the KS fractured blocks can be determined by a measured pressure curve analysis of the hydraulic support of the mine’s working surface [34]. Hence, Qm−1(0) in Figure 6 is(1)Qm−10=lm−1qm−12.Figure 6
Stress analysis of the fissure zone KS on the coal mass side.y
1(x), y2(x), …, ym−1(x), ym(x), …, yn−1(x), yn(x) and ycm(x), …, ycn−1(x), ycn(x) are the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf, respectively, and are measured in m. For example, ym−1(x) is the deflection curve equation of KSm−1 on the coal mass side. p1(x), p2(x), …, pm−1(x), pm(x), …, pn−1(x), pn(x) and pcm(x), …, pcn−1(x), pcn(x) represent the force on the foundation exerted by the different layers of KS on the side of the coal mass and the side of the goaf, which is also the force exerted by KS in different locations on the adjacent KS below them. This is measured in MPa. As shown in Figure 6, pm(x) is the force of KSm on the foundation and also the load of KSm acting on KSm−1. Hence, pm−1(x) and pm(x) in Figure 6 are calculated as follows:(2)pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x.


### 3.1.2. Stress Analysis of the Curve Subsidence Zone KS

Using KSn−1 in Figure 5 as an example, the stress analysis is shown in Figures 7(a) and 7(b). M1(x), M2(x), …, Mm−1(x), Mm(x), …, Mn−1(x), Mn(x), and Mcm(x), …, Mcn−1(x), Mcn(x) represent the bending moments borne by different layers of KS on the side of the coal mass and the side of the goaf, respectively. As shown in Figure 7, Mn−1(x) and Mcn−1(x) are the bending moments of KSn−1 on the side of the coal mass and side of the goaf, respectively. Qn−1(0) and Qcn−1(0) are the shearing forces of KSn−1 on the side of the coal mass and the side of the goaf at x = 0, respectively.Figure 7
Stress analysis of the curve subsidence zone KS: (a) side of the coal mass; (b) side of the goaf.
(a)
(b)p
n−1(x) and pn(x) in Figure 7(a) are as follows:(3)pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x.p
cn−1(x) and pcn(x) in Figure 7(b) are as follows:(4)pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x.The relationship between the elastic foundation beam’s deflection, foundation pressure, and the load it bears is as follows [15, 16, 20, 30, 32, 33]:(5)EId4yxdx4+px=qx,where E is the elastic modulus of the beam, I is its moment of inertia, y(x) is its deflection equation, p(x) is the foundation pressure, and q(x) is the load that the beam bears.According to the established mechanical model, the following groups of deflection curve differential equations can be formed:(6)E1I1d4y1xdx4+bk1y1x=bq1+bk2y2x−y1x,E2I2d4y2xdx4+bk2y2x−y1x=bq2+bk3y3x−y2x,⋮Em−1Im−1d4ym−1xdx4+bkm−1ym−1x−ym−2x=bqm−1+bkmymx−ym−1x,EmImd4ymxdx4+bkmymx−ym−1x=bqm+bkm+1ym+1x−ymx,⋮En−1In−1d4yn−1xdx4+bkn−1yn−1x−yn−2x=bqn−1+bknynx−yn−1x,EnInd4ynxdx4+bknynx−yn−1x=bqn,x≤0,(7)EmImd4ycmxdx4+bkcmycmx=bqm+bkm+1ycm+1x−ycmx,⋮En−1In−1d4ycn−1xdx4+bkn−1ycn−1x−ycn−2x=bqn−1+bknycnx−ycn−1x,EnInd4ycnxdx4+bknycnx−ycn−1x=bqn,x≥0,where b(m) is the width of the beam, considered to be equal to unity, i.e., b = 1.Equation group (6) is the deflection curve differential equation for different KS on the side of the coal mass. Equation group (7) is the deflection curve differential equation of different KS on the side of the goaf. Solving the equations in equation group (6) gives(8)and4ny1xdx4n+an−1d4n−1y1xdx4n−1+⋯+a2d8y1xdx8+a1d4y1xdx4+an+1y1x=an+2.As the expressionsa1, a2, …, an, an+1, an+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; k1, k2, …, kn−1, kn; and q1, q2, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (6).The characteristic equation corresponding to equation (8) is(9)anr4n+an−1r4n−1+⋯+a2r8+a1r4+an+1=0.The root of equation (9) is r1,2,3,4 = ±t1 ± t1 · i, r5,6,7,8 = ±t2 ± t2 · i,  …,  r4(n−1)−3,4(n−1)−2,4(n−1)−1,4(n−1) = ±tn−1 ± tn−1 · i, and r4n−3,4n−2,4n−1,4n = ±tn ± tn · i. As the expressions t1, t2, …, tn−1, tn are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into equation (9).A special solution for equation (8) is(10)y∗=q1+q2+⋯+qn−1+qnk1.Then, the general solution for equation (8) is(11)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A3e−x⋅t1cosx⋅t1+A4e−x⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+A7e−x⋅t2cosx⋅t2+A8e−x⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−1−1e−x⋅tn−1cosx⋅tn−1+A4n−1e−x⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+A4n−1e−x⋅tncosx⋅tn+A4ne−x⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1.Asx ⟶ −∞ and y1(x) ⟶ (q1 + q2 + ⋯ + qn−1 + qn)/k1. Thus, it is only satisfied when A3, A4, A7, A8, …, A4(n−1)−1, A4(n−1), A4n−1, A4n all have a value of 0, and the deflection curve equation y1(x) is obtained. Substituting this into equation group (6), the deflection curve equations of the different layers of KS on the side of the coal mass are shown in equation (12).As the expressionsB2,1, B2,2, …, B2,n−1, B2,n; …; Bn−1,1, Bn−1,2, …, Bn−1,n−1, Bn−1,n; Bn,1, Bn,2, …, Bn,n−1, Bn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into group (6). The deflection curve equation group of the different layers of KS on the side of the coal mass has 2n unknowns, A1, A2, A5, A6, …, A4(n−1)−3, A4(n−1)−2, A4n−3, A4n−2:(12)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1,y2x=B2,1A1ex⋅t1cosx⋅t1+B2,1A2ex⋅t1sinx⋅t1+B2,2A5ex⋅t2cosx⋅t2+B2,2A6ex⋅t2sinx⋅t2+⋯+B2,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+B2,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+B2,nA4n−3ex⋅tncosx⋅tn+B2,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2,⋮ym−1x=Bm−1,1A1ex⋅t1cosx⋅t1+Bm−1,1A2ex⋅t1sinx⋅t1+Bm−1,2A5ex⋅t2cosx⋅t2+Bm−1,2A6ex⋅t2sinx⋅t2+⋯+Bm−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm−1,nA4n−3ex⋅tncosx⋅tn+Bm−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1,ymx=Bm,1A1ex⋅t1cosx⋅t1+Bm,1A2ex⋅t1sinx⋅t1+Bm,2A5ex⋅t2cosx⋅t2+Bm,2A6ex⋅t2sinx⋅t2+⋯+Bm,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm,nA4n−3ex⋅tncosx⋅tn+Bm,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1+qm+⋯+qn−1+qnkm,⋮yn−1x=Bn−1,1A1ex⋅t1cosx⋅t1+Bn−1,1A2ex⋅t1sinx⋅t1+Bn−1,2A5ex⋅t2cosx⋅t2+Bn−1,2A6ex⋅t2sinx⋅t2+⋯+Bn−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn−1,nA4n−3ex⋅tncosx⋅tn+Bn−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1,ynx=Bn,1A1ex⋅t1cosx⋅t1+Bn,1A2ex⋅t1sinx⋅t1+Bn,2A5ex⋅t2cosx⋅t2+Bn,2A6ex⋅t2sinx⋅t2+⋯+Bn,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn,nA4n−3ex⋅tncosx⋅tn+Bn,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1+qnkn.Solving the equations in equation group (7) gives(13)bnd4n−m+1ycmxdx4n−m+1+bn−1d4n−mycmxdx4n−m+⋯+bm+1d8ycmxdx8+bmd4ycmxdx4+bn+1ycmx=bn+2.As the expressionsbm, bm+1, …, bn, bn+1, bn+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; kcm, km+1, …, kn−1, kn; and qm, qm+1, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (7).The characteristic equation corresponding to equation (13) is(14)bnr4n−m+1+bn−1r4n−m+⋯+bm+1r8+bmr4+bn+1=0.The root of equation (14) is r4m−3,4m−2,4m−1,4m = ±sm ± sm · i, …, r4(n−1)−3, 4(n−1)−2, 4(n−1)−1, 4(n−1) = ±sn−1 ± sn−1 · i, r4n−3,4n−2,4n−1,4n = ±sn ± sn · i. As the expressions sm, …, sn−1, sn are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into equation (14).A special solution for equation (13) is(15)y∗=qm+⋯+qn−1+qnkcm.Then, the general solution for equation (13) is(16)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,(17)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,⋮ycn−1x=Dn−1,mC4m−3ex⋅smcosx⋅sm+Dn−1,mC4m−2ex⋅smsinx⋅sm+Dn−1,mC4m−1e−x⋅sm·cosx⋅sm+Dn−1,mC4me−x⋅smsinx⋅sm+⋯+Dn−1,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn−1,n−1C4n−1−2ex⋅sn−1sinx⋅sn−1+Dn−1,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn−1,n−1·C4n−1e−x⋅sn−1sinx⋅sn−1+Dn−1,nC4n−3ex⋅sncosx⋅sn+Dn−1,nC4n−2ex⋅snsinx⋅sn+Dn−1,nC4n−1e−x⋅sncosx⋅sn+Dn−1,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1,ycnx=Dn,mC4m−3ex⋅smcosx⋅sm+Dn,mC4m−2ex⋅smsinx⋅sm+Dn,mC4m−1e−x⋅smcosx⋅sm+Dn,mC4me−x⋅smsinx⋅sm+⋯+Dn,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1−2·ex⋅sn−1sinx⋅sn−1+Dn,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1e−x⋅sn−1sinx⋅sn−1+Dn,nC4n−3ex⋅sncosx⋅sn+Dn,nC4n−2ex⋅snsinx⋅sn+Dn,nC4n−1e−x⋅sncosx⋅sn+Dn,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1+qnkn.Substituting equation (16) into equation group (7), the deflection curve equations of the different layers of KS on the side of the goaf are shown in equation (17). As the expressions Dm+1,m, …, Dm+1,n−1, Dm+1,n; …; Dn−1,m, …, Dn−1,n−1, Dn−1,n; Dn,m, …, D n,n−1, Dn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into group (7). The deflection curve equation group of the different layers of KS on the side of the goaf has 4(n−m+1) unknowns, C4m−3, C4m−2, C4m−1, C4m; …; C4(n−1)−3, C4(n−1)−2, C4(n−1)−1, C4(n−1);C4n−3, C4n−2, C4n−1, C4n.In the paper,θ1(x), θ2(x), …, θm−1(x), θm(x), …, θn−1(x), θn(x), and θcm(x), …, θcn−1(x), θcn(x) represent the rotation angles experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. From the relationship among the rotation angle, bending moment, and shearing force of the KS with its deflection curve equation in equations (12) and (17) [15, 30, 32], by substituting the deflection curve equations of the different layers of KS into equation (18), we obtain θ1(x), …, θn−1(x), θn(x); M1(x), …, Mn−1(x), Mn(x); Q1(x), …, Qn−1(x), Qn(x); θcm(x), …, θcn−1(x), θcn(x); Mcm(x), …, Mcn−1(x), Mcn(x); and Qcm(x), …, Qcn−1(x), Qcn(x):(18)θx=dyxdx,Mx=−EIdy2xdx2,Qx=−EIdy3xdx3,(19)ym0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qm−1+qm+⋯+qn−1+qnkm−1−qm+⋯+qn−1+qnkm=ycm0,⋮yn−10−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1=ycn−10−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1,yn0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1−qnkn=ycn0−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1−qnkn,θm0=θcm0,…,θn−10=θcn−10,θn0=θcn0,Mm0=Mcm0,…,Mn−10=Mcn−10,Mn0=Mcn0,Qm0=Qcm0,…,Qn−10=Qcn−10,Qn0=Qcn0,Q10=l1q12,Q20=l2q22,…,Qm−10=lm−1qm−12;M10=0,M20=0,…,Mm−10=0,θcmL2=0,…,θcn−1L2=0,θcnL2=0,QcmL2=0,…,Qcn−1L2=0,QcnL2=0.There are a total of 4(n − m + 1) + 2n unknowns in equations (12) and (17), which are the deflection curve equations of the KS on the side of the coal mass and on the side of the goaf, respectively. The different curve subsidence zone KS on the side of the coal mass and the side of the goaf are continuous at x = 0. At x = 0, the KS on both sides have equal subsidence, rotation angle, bending moment, and shearing force. In the fissure zone, different KS on the side of the coal mass have a bending moment of 0 at x = 0, and the shearing force at x = 0 can be determined using equation (6). The model is symmetric about x = L/2; thus, the shearing force and rotation angle of the KS at x = L/2 is 0. Considering that in equations (12) and (17) (q1 + q2 + ⋯ + qn−1 + qn)/k1, (q2 + ⋯ + qn−1 + qn)/k2, …, (qm−1 + qm + ⋯ + qn−1 + qn)/km−1, (qm + ⋯ + qn−1 + qn)/km, …, (qn−1 + qn)/kn−1, qn/kn are caused by stratal compression before the coal seam is mined, the continuity and boundary conditions of the different KS are shown in equation (19).There are a total of 2n + 4(n − m + 1) equations that can be used to obtain the unknowns in equations (12) and (17) and obtain the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf. In particular, these unknown values will no longer vary with the increase of L after the land surface has fully subsided. The foundation pressure equations of the different layers of KS on the side of the coal mass and the side of the goaf are(20)p1x=k1y1x,p2x=k2y2x−y1x,⋮pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x,⋮pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x,x≤0,pcmx=kcmycmx,⋮pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x,x≥0.


## 3.1.1. Stress Analysis of the Fissure Zone KS on the Coal Mass Side

Using KSm−1 in Figure 5 as an example, the corresponding stress analysis is shown in Figure 6. In this study, Q1(x), Q2(x), …, Qm−1(x), Qm(x), …, Qn−1(x), Qn(x) and Qcm(x), …, Qcn−1(x), Qcn(x) represent the shearing forces experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. For example, in Figure 6, Qm−1(0) is the shearing force borne by the side of the coal mass KSm−1 at x = 0. The shearing force of the different layers of KS in the fissure zone at x = 0 is approximately equal to half of the product of the lengths of the fractured KS blocks and their load in the goaf [9]. The lengths of the KS fractured blocks A1, A2, …, Am−1 are represented by l1, l2, …, lm−1, respectively, and can be calculated using the KS discriminant method [7, 11, 12]. If the number of KS layers in the fissure zone is small, the lengths of the KS fractured blocks can be determined by a measured pressure curve analysis of the hydraulic support of the mine’s working surface [34]. Hence, Qm−1(0) in Figure 6 is(1)Qm−10=lm−1qm−12.Figure 6
Stress analysis of the fissure zone KS on the coal mass side.y
1(x), y2(x), …, ym−1(x), ym(x), …, yn−1(x), yn(x) and ycm(x), …, ycn−1(x), ycn(x) are the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf, respectively, and are measured in m. For example, ym−1(x) is the deflection curve equation of KSm−1 on the coal mass side. p1(x), p2(x), …, pm−1(x), pm(x), …, pn−1(x), pn(x) and pcm(x), …, pcn−1(x), pcn(x) represent the force on the foundation exerted by the different layers of KS on the side of the coal mass and the side of the goaf, which is also the force exerted by KS in different locations on the adjacent KS below them. This is measured in MPa. As shown in Figure 6, pm(x) is the force of KSm on the foundation and also the load of KSm acting on KSm−1. Hence, pm−1(x) and pm(x) in Figure 6 are calculated as follows:(2)pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x.


## 3.1.2. Stress Analysis of the Curve Subsidence Zone KS

Using KSn−1 in Figure 5 as an example, the stress analysis is shown in Figures 7(a) and 7(b). M1(x), M2(x), …, Mm−1(x), Mm(x), …, Mn−1(x), Mn(x), and Mcm(x), …, Mcn−1(x), Mcn(x) represent the bending moments borne by different layers of KS on the side of the coal mass and the side of the goaf, respectively. As shown in Figure 7, Mn−1(x) and Mcn−1(x) are the bending moments of KSn−1 on the side of the coal mass and side of the goaf, respectively. Qn−1(0) and Qcn−1(0) are the shearing forces of KSn−1 on the side of the coal mass and the side of the goaf at x = 0, respectively.Figure 7
Stress analysis of the curve subsidence zone KS: (a) side of the coal mass; (b) side of the goaf.
(a)
(b)p
n−1(x) and pn(x) in Figure 7(a) are as follows:(3)pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x.p
cn−1(x) and pcn(x) in Figure 7(b) are as follows:(4)pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x.The relationship between the elastic foundation beam’s deflection, foundation pressure, and the load it bears is as follows [15, 16, 20, 30, 32, 33]:(5)EId4yxdx4+px=qx,where E is the elastic modulus of the beam, I is its moment of inertia, y(x) is its deflection equation, p(x) is the foundation pressure, and q(x) is the load that the beam bears.According to the established mechanical model, the following groups of deflection curve differential equations can be formed:(6)E1I1d4y1xdx4+bk1y1x=bq1+bk2y2x−y1x,E2I2d4y2xdx4+bk2y2x−y1x=bq2+bk3y3x−y2x,⋮Em−1Im−1d4ym−1xdx4+bkm−1ym−1x−ym−2x=bqm−1+bkmymx−ym−1x,EmImd4ymxdx4+bkmymx−ym−1x=bqm+bkm+1ym+1x−ymx,⋮En−1In−1d4yn−1xdx4+bkn−1yn−1x−yn−2x=bqn−1+bknynx−yn−1x,EnInd4ynxdx4+bknynx−yn−1x=bqn,x≤0,(7)EmImd4ycmxdx4+bkcmycmx=bqm+bkm+1ycm+1x−ycmx,⋮En−1In−1d4ycn−1xdx4+bkn−1ycn−1x−ycn−2x=bqn−1+bknycnx−ycn−1x,EnInd4ycnxdx4+bknycnx−ycn−1x=bqn,x≥0,where b(m) is the width of the beam, considered to be equal to unity, i.e., b = 1.Equation group (6) is the deflection curve differential equation for different KS on the side of the coal mass. Equation group (7) is the deflection curve differential equation of different KS on the side of the goaf. Solving the equations in equation group (6) gives(8)and4ny1xdx4n+an−1d4n−1y1xdx4n−1+⋯+a2d8y1xdx8+a1d4y1xdx4+an+1y1x=an+2.As the expressionsa1, a2, …, an, an+1, an+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; k1, k2, …, kn−1, kn; and q1, q2, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (6).The characteristic equation corresponding to equation (8) is(9)anr4n+an−1r4n−1+⋯+a2r8+a1r4+an+1=0.The root of equation (9) is r1,2,3,4 = ±t1 ± t1 · i, r5,6,7,8 = ±t2 ± t2 · i,  …,  r4(n−1)−3,4(n−1)−2,4(n−1)−1,4(n−1) = ±tn−1 ± tn−1 · i, and r4n−3,4n−2,4n−1,4n = ±tn ± tn · i. As the expressions t1, t2, …, tn−1, tn are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into equation (9).A special solution for equation (8) is(10)y∗=q1+q2+⋯+qn−1+qnk1.Then, the general solution for equation (8) is(11)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A3e−x⋅t1cosx⋅t1+A4e−x⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+A7e−x⋅t2cosx⋅t2+A8e−x⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−1−1e−x⋅tn−1cosx⋅tn−1+A4n−1e−x⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+A4n−1e−x⋅tncosx⋅tn+A4ne−x⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1.Asx ⟶ −∞ and y1(x) ⟶ (q1 + q2 + ⋯ + qn−1 + qn)/k1. Thus, it is only satisfied when A3, A4, A7, A8, …, A4(n−1)−1, A4(n−1), A4n−1, A4n all have a value of 0, and the deflection curve equation y1(x) is obtained. Substituting this into equation group (6), the deflection curve equations of the different layers of KS on the side of the coal mass are shown in equation (12).As the expressionsB2,1, B2,2, …, B2,n−1, B2,n; …; Bn−1,1, Bn−1,2, …, Bn−1,n−1, Bn−1,n; Bn,1, Bn,2, …, Bn,n−1, Bn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to E1, E2, …, En−1, En; I1, I2, …, In−1, In; and k1, k2, …, kn−1, kn. We can obtain them by substituting these parameters into group (6). The deflection curve equation group of the different layers of KS on the side of the coal mass has 2n unknowns, A1, A2, A5, A6, …, A4(n−1)−3, A4(n−1)−2, A4n−3, A4n−2:(12)y1x=A1ex⋅t1cosx⋅t1+A2ex⋅t1sinx⋅t1+A5ex⋅t2cosx⋅t2+A6ex⋅t2sinx⋅t2+⋯+A4n−1−3ex⋅tn−1cosx⋅tn−1+A4n−1−2ex⋅tn−1sinx⋅tn−1+A4n−3ex⋅tncosx⋅tn+A4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1,y2x=B2,1A1ex⋅t1cosx⋅t1+B2,1A2ex⋅t1sinx⋅t1+B2,2A5ex⋅t2cosx⋅t2+B2,2A6ex⋅t2sinx⋅t2+⋯+B2,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+B2,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+B2,nA4n−3ex⋅tncosx⋅tn+B2,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2,⋮ym−1x=Bm−1,1A1ex⋅t1cosx⋅t1+Bm−1,1A2ex⋅t1sinx⋅t1+Bm−1,2A5ex⋅t2cosx⋅t2+Bm−1,2A6ex⋅t2sinx⋅t2+⋯+Bm−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm−1,nA4n−3ex⋅tncosx⋅tn+Bm−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1,ymx=Bm,1A1ex⋅t1cosx⋅t1+Bm,1A2ex⋅t1sinx⋅t1+Bm,2A5ex⋅t2cosx⋅t2+Bm,2A6ex⋅t2sinx⋅t2+⋯+Bm,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bm,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bm,nA4n−3ex⋅tncosx⋅tn+Bm,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qm−1+qm+⋯+qn−1+qnkm−1+qm+⋯+qn−1+qnkm,⋮yn−1x=Bn−1,1A1ex⋅t1cosx⋅t1+Bn−1,1A2ex⋅t1sinx⋅t1+Bn−1,2A5ex⋅t2cosx⋅t2+Bn−1,2A6ex⋅t2sinx⋅t2+⋯+Bn−1,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn−1,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn−1,nA4n−3ex⋅tncosx⋅tn+Bn−1,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1,ynx=Bn,1A1ex⋅t1cosx⋅t1+Bn,1A2ex⋅t1sinx⋅t1+Bn,2A5ex⋅t2cosx⋅t2+Bn,2A6ex⋅t2sinx⋅t2+⋯+Bn,n−1A4n−1−3ex⋅tn−1cosx⋅tn−1+Bn,n−1A4n−1−2ex⋅tn−1sinx⋅tn−1+Bn,nA4n−3ex⋅tncosx⋅tn+Bn,nA4n−2ex⋅tnsinx⋅tn+q1+q2+⋯+qn−1+qnk1+q2+⋯+qn−1+qnk2+⋯+qn−1+qnkn−1+qnkn.Solving the equations in equation group (7) gives(13)bnd4n−m+1ycmxdx4n−m+1+bn−1d4n−mycmxdx4n−m+⋯+bm+1d8ycmxdx8+bmd4ycmxdx4+bn+1ycmx=bn+2.As the expressionsbm, bm+1, …, bn, bn+1, bn+2 are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; kcm, km+1, …, kn−1, kn; and qm, qm+1, …, qn−1, qn. We can obtain them by substituting these parameters into equation group (7).The characteristic equation corresponding to equation (13) is(14)bnr4n−m+1+bn−1r4n−m+⋯+bm+1r8+bmr4+bn+1=0.The root of equation (14) is r4m−3,4m−2,4m−1,4m = ±sm ± sm · i, …, r4(n−1)−3, 4(n−1)−2, 4(n−1)−1, 4(n−1) = ±sn−1 ± sn−1 · i, r4n−3,4n−2,4n−1,4n = ±sn ± sn · i. As the expressions sm, …, sn−1, sn are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into equation (14).A special solution for equation (13) is(15)y∗=qm+⋯+qn−1+qnkcm.Then, the general solution for equation (13) is(16)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,(17)ycmx=C4m−3ex⋅smcosx⋅sm+C4m−2ex⋅smsinx⋅sm+C4m−1e−x⋅smcosx⋅sm+C4me−x⋅smsinx⋅sm+⋯+C4n−1−3ex⋅sn−1cosx⋅sn−1+C4n−1−2ex⋅sn−1sinx⋅sn−1+C4n−1−1e−x⋅sn−1cosx⋅sn−1+C4n−1e−x⋅sn−1sinx⋅sn−1+C4n−3ex⋅sncosx⋅sn+C4n−2ex⋅snsinx⋅sn+C4n−1e−x⋅sncosx⋅sn+C4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm,⋮ycn−1x=Dn−1,mC4m−3ex⋅smcosx⋅sm+Dn−1,mC4m−2ex⋅smsinx⋅sm+Dn−1,mC4m−1e−x⋅sm·cosx⋅sm+Dn−1,mC4me−x⋅smsinx⋅sm+⋯+Dn−1,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn−1,n−1C4n−1−2ex⋅sn−1sinx⋅sn−1+Dn−1,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn−1,n−1·C4n−1e−x⋅sn−1sinx⋅sn−1+Dn−1,nC4n−3ex⋅sncosx⋅sn+Dn−1,nC4n−2ex⋅snsinx⋅sn+Dn−1,nC4n−1e−x⋅sncosx⋅sn+Dn−1,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1,ycnx=Dn,mC4m−3ex⋅smcosx⋅sm+Dn,mC4m−2ex⋅smsinx⋅sm+Dn,mC4m−1e−x⋅smcosx⋅sm+Dn,mC4me−x⋅smsinx⋅sm+⋯+Dn,n−1C4n−1−3ex⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1−2·ex⋅sn−1sinx⋅sn−1+Dn,n−1C4n−1−1e−x⋅sn−1cosx⋅sn−1+Dn,n−1C4n−1e−x⋅sn−1sinx⋅sn−1+Dn,nC4n−3ex⋅sncosx⋅sn+Dn,nC4n−2ex⋅snsinx⋅sn+Dn,nC4n−1e−x⋅sncosx⋅sn+Dn,nC4ne−x⋅snsinx⋅sn+qm+⋯+qn−1+qnkcm+qm+1+⋯+qn−1+qnkm+1+⋯+qn−1+qnkn−1+qnkn.Substituting equation (16) into equation group (7), the deflection curve equations of the different layers of KS on the side of the goaf are shown in equation (17). As the expressions Dm+1,m, …, Dm+1,n−1, Dm+1,n; …; Dn−1,m, …, Dn−1,n−1, Dn−1,n; Dn,m, …, D n,n−1, Dn,n are relatively complex, the actual equations are not listed here, but they are all real constants related to Em, Em+1, …, En−1, En; Im, Im+1, …, In−1, In; and kcm, km+1, …, kn−1, kn. We can obtain them by substituting these parameters into group (7). The deflection curve equation group of the different layers of KS on the side of the goaf has 4(n−m+1) unknowns, C4m−3, C4m−2, C4m−1, C4m; …; C4(n−1)−3, C4(n−1)−2, C4(n−1)−1, C4(n−1);C4n−3, C4n−2, C4n−1, C4n.In the paper,θ1(x), θ2(x), …, θm−1(x), θm(x), …, θn−1(x), θn(x), and θcm(x), …, θcn−1(x), θcn(x) represent the rotation angles experienced by different layers of KS on the side of the coal mass and the side of the goaf, respectively. From the relationship among the rotation angle, bending moment, and shearing force of the KS with its deflection curve equation in equations (12) and (17) [15, 30, 32], by substituting the deflection curve equations of the different layers of KS into equation (18), we obtain θ1(x), …, θn−1(x), θn(x); M1(x), …, Mn−1(x), Mn(x); Q1(x), …, Qn−1(x), Qn(x); θcm(x), …, θcn−1(x), θcn(x); Mcm(x), …, Mcn−1(x), Mcn(x); and Qcm(x), …, Qcn−1(x), Qcn(x):(18)θx=dyxdx,Mx=−EIdy2xdx2,Qx=−EIdy3xdx3,(19)ym0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qm−1+qm+⋯+qn−1+qnkm−1−qm+⋯+qn−1+qnkm=ycm0,⋮yn−10−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1=ycn−10−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1,yn0−q1+q2+⋯+qn−1+qnk1−q2+⋯+qn−1+qnk2−⋯−qn−1+qnkn−1−qnkn=ycn0−qm+1+⋯+qn−1+qnkm+1−⋯−qn−1+qnkn−1−qnkn,θm0=θcm0,…,θn−10=θcn−10,θn0=θcn0,Mm0=Mcm0,…,Mn−10=Mcn−10,Mn0=Mcn0,Qm0=Qcm0,…,Qn−10=Qcn−10,Qn0=Qcn0,Q10=l1q12,Q20=l2q22,…,Qm−10=lm−1qm−12;M10=0,M20=0,…,Mm−10=0,θcmL2=0,…,θcn−1L2=0,θcnL2=0,QcmL2=0,…,Qcn−1L2=0,QcnL2=0.There are a total of 4(n − m + 1) + 2n unknowns in equations (12) and (17), which are the deflection curve equations of the KS on the side of the coal mass and on the side of the goaf, respectively. The different curve subsidence zone KS on the side of the coal mass and the side of the goaf are continuous at x = 0. At x = 0, the KS on both sides have equal subsidence, rotation angle, bending moment, and shearing force. In the fissure zone, different KS on the side of the coal mass have a bending moment of 0 at x = 0, and the shearing force at x = 0 can be determined using equation (6). The model is symmetric about x = L/2; thus, the shearing force and rotation angle of the KS at x = L/2 is 0. Considering that in equations (12) and (17) (q1 + q2 + ⋯ + qn−1 + qn)/k1, (q2 + ⋯ + qn−1 + qn)/k2, …, (qm−1 + qm + ⋯ + qn−1 + qn)/km−1, (qm + ⋯ + qn−1 + qn)/km, …, (qn−1 + qn)/kn−1, qn/kn are caused by stratal compression before the coal seam is mined, the continuity and boundary conditions of the different KS are shown in equation (19).There are a total of 2n + 4(n − m + 1) equations that can be used to obtain the unknowns in equations (12) and (17) and obtain the deflection curve equations of the different layers of KS on the side of the coal mass and the side of the goaf. In particular, these unknown values will no longer vary with the increase of L after the land surface has fully subsided. The foundation pressure equations of the different layers of KS on the side of the coal mass and the side of the goaf are(20)p1x=k1y1x,p2x=k2y2x−y1x,⋮pm−1x=km−1ym−1x−ym−2x,pmx=kmymx−ym−1x,⋮pn−1x=kn−1yn−1x−yn−2x,pnx=knynx−yn−1x,x≤0,pcmx=kcmycmx,⋮pcn−1x=kn−1ycn−1x−ycn−2x,pcnx=knycnx−ycn−1x,x≥0.


## 3.2. Method for Calculating Abutment Pressure

### 3.2.1. Calculation Method for Coal Mass Abutment Pressure

According to Section3.1, the force exerted by the different layers of KS on the coal mass is p1(x). Considering the weight of the rock mass on the side of the caving zone, q0, the coal mass stress is(21)σ1x=p1x+q0.When the coal lateral pressure in the proximity of the mining boundary is zero, i.e., when the area is under uniaxial compression, as the plastic failure [7, 29, 35] moves inward, the coal will gradually become triaxially stressed, thereby enhancing its compressive strength [7]. If the coal strength on the coal mass side near the mining boundary is smaller than the abutment pressure value in the corresponding location, which can be calculated using (21), a plastic failure will occur in this part of the coal. In this case, the equation for the abutment pressure in the plastic failure region (equation (22)), can be obtained from Ref. [7]:(22)σpx=N0e−2fx/M1+sinφ/1−sinφ,x≤0.The symbolσp(x) represents the abutment pressure of the coal plastic region, M is the mining height of the seam, φ is the internal angle of friction between the seam and the roof/floor slate stratum, f is the friction factor between the stratum surfaces, and N0 is the supporting capacity of the coal wall.As (22) is equal to (21), the width |x1| of the coal plastic region can be estimated, along with the peak abutment pressure σ(x1) on the interface between the coal plastic region and the elastic region. The abutment pressure in the plastic region and in the elastic region can be calculated using (21) and (22), respectively. Hence, the coal mass abutment pressure is(23)σx=σ1x,x≤x1,σpx,x1≤x≤0.Ifσ(x1) = K(q0 + q1 + q2 + ⋯ + qn−1 + qn), then K is the stress concentration coefficient. If σ1(xd) = q0 + q1 + q2 + ⋯ + qn−1 + qn, then |xd| is the abutment pressure width of the side of the coal mass.


### 3.2.2. Calculation Method for Goaf Abutment Pressure

According to Section3.1, the stress distribution of the curve subsidence zone KS in the goaf is pcm(x). By superimposing it on the stress distribution caused by the fissure zone and the caving zone fractured rock masses in the goaf, we can obtain the distribution of abutment pressure in the goaf. It is assumed that the KS fractured blocks in the fissure zone far away from the mining boundary transfer the entire weight to the goaf [9], and the KS fractured blocks at the mining boundary A1, …, Am−1 transfer half of their weight to the goaf [9, 23]. The load distribution that they form in the goaf increases linearly [19]. Thus, the stress distribution formed by the fissure zone KS in the goaf is as follows:(24)σfx=q1l1+q2l2+⋯+qm−1lm−1x,0<x≤l1,q2l2+⋯+qm−1lm−1x+q1,l1≤x≤l2,⋮⋮qm−1lm−1x+q1+q2+⋯+qm−2,lm−2≤x≤lm−1,q1+q2+⋯+qm−2+qm−1,x≥lm−1.The distribution of the abutment pressure in the goaf can be calculated as follows:(25)σcx=pcmx+σfx+q0,x>0.


## 3.2.1. Calculation Method for Coal Mass Abutment Pressure

According to Section3.1, the force exerted by the different layers of KS on the coal mass is p1(x). Considering the weight of the rock mass on the side of the caving zone, q0, the coal mass stress is(21)σ1x=p1x+q0.When the coal lateral pressure in the proximity of the mining boundary is zero, i.e., when the area is under uniaxial compression, as the plastic failure [7, 29, 35] moves inward, the coal will gradually become triaxially stressed, thereby enhancing its compressive strength [7]. If the coal strength on the coal mass side near the mining boundary is smaller than the abutment pressure value in the corresponding location, which can be calculated using (21), a plastic failure will occur in this part of the coal. In this case, the equation for the abutment pressure in the plastic failure region (equation (22)), can be obtained from Ref. [7]:(22)σpx=N0e−2fx/M1+sinφ/1−sinφ,x≤0.The symbolσp(x) represents the abutment pressure of the coal plastic region, M is the mining height of the seam, φ is the internal angle of friction between the seam and the roof/floor slate stratum, f is the friction factor between the stratum surfaces, and N0 is the supporting capacity of the coal wall.As (22) is equal to (21), the width |x1| of the coal plastic region can be estimated, along with the peak abutment pressure σ(x1) on the interface between the coal plastic region and the elastic region. The abutment pressure in the plastic region and in the elastic region can be calculated using (21) and (22), respectively. Hence, the coal mass abutment pressure is(23)σx=σ1x,x≤x1,σpx,x1≤x≤0.Ifσ(x1) = K(q0 + q1 + q2 + ⋯ + qn−1 + qn), then K is the stress concentration coefficient. If σ1(xd) = q0 + q1 + q2 + ⋯ + qn−1 + qn, then |xd| is the abutment pressure width of the side of the coal mass.


## 3.2.2. Calculation Method for Goaf Abutment Pressure

According to Section3.1, the stress distribution of the curve subsidence zone KS in the goaf is pcm(x). By superimposing it on the stress distribution caused by the fissure zone and the caving zone fractured rock masses in the goaf, we can obtain the distribution of abutment pressure in the goaf. It is assumed that the KS fractured blocks in the fissure zone far away from the mining boundary transfer the entire weight to the goaf [9], and the KS fractured blocks at the mining boundary A1, …, Am−1 transfer half of their weight to the goaf [9, 23]. The load distribution that they form in the goaf increases linearly [19]. Thus, the stress distribution formed by the fissure zone KS in the goaf is as follows:(24)σfx=q1l1+q2l2+⋯+qm−1lm−1x,0<x≤l1,q2l2+⋯+qm−1lm−1x+q1,l1≤x≤l2,⋮⋮qm−1lm−1x+q1+q2+⋯+qm−2,lm−2≤x≤lm−1,q1+q2+⋯+qm−2+qm−1,x≥lm−1.The distribution of the abutment pressure in the goaf can be calculated as follows:(25)σcx=pcmx+σfx+q0,x>0.


## 3.3. Determining Model Parameters

E
1, E2, …, Em−1, Em, …, En−1, En and H1, H2, …, Hm−1, Hm, …, Hn−1, Hn as well as L and M can be determined using the physical mechanical parameter tests on drill core samples and the specific mining conditions.


### 3.3.1. Determination ofI1, I2, …, Im−1, Im, …, In−1, In

I
m is used as an example to illustrate the determination of the moment of inertia for different layers of KS. Im is derived using the following equation:(26)Im=112bHm3.


### 3.3.2. Determination ofq1, q2, …, qm−1, qm, …, qn−1, qn

q
m−1 is used as an example to illustrate the respective load for different KS and the load layer weight they control. If the adjacent KS interlayer rock layers are named z1, z2, …, zi−1, zi, with a total of i layers, as shown in Figure 8, then qm−1 is(27)qm−1=γm−1Hm−1+∑i=1iγziHzi.Figure 8
Diagram of adjacent KS.In Figure8, γz1, γz2, …, γzi−1, γzi represent the density of the rock mass z1, z2, …, zi−1, zi, respectively; Ez1, Ez2, …, Ezi−1, Ezi represent the elastic modulus of the rock mass z1, z2, …, zi−1, zi, respectively; Hz1, Hz2, …, Hzi−1, Hzi represent the thickness of the rock mass z1, z2, …, zi−1, zi, respectively. In this study, γ1, γ2, …, γm−1, γm, …, γn−1, γn represent the density of the different layers of KS.


### 3.3.3. Determination ofk1, k2, …, km−1, km, …, kn−1, kn

k
m is used as an example to illustrate the calculation method of k1, k2, …, km−1, km, …, kn−1, kn. The foundation coefficient [15, 16, 20, 30] formed by the total of i layers and KSm−1 shown in Figure 8 is(28)km=1Hm−1/Em−1+∑i=1iHzi/Ezi.In particular, the calculation ofk1 should consider the strata from KS1 to the thick hard rock formation in the coal seam floor [30].


### 3.3.4. Determination ofkcm

Because the model assumes that the fractured rock in the caving zone and fissure zone satisfies Winkler’s elastic foundation hypothesis [15, 16, 20, 30, 32, 33], we can calculate the foundation coefficient kcm of KSm on the goaf side if we know the subsidence value of KSm and the corresponding load bearing of the foundation [30]. When the land surface is under full subsidence [31], the level of subsidence in the curve subsidence zone is consistent [25], and the subsidence of KSm reaches the maximum subsidence [31], δM. At this point, the load of the foundation that corresponds to the maximum subsidence of the side of the goaf, KSm, returns to primary rock stress [7, 11] qm + qm+1 + ⋯ + qn−1 + qn. Hence, the foundation coefficient of the KSm on the side of the goaf, kcm, is(29)kcm=qm+qm+1+…+qn−1+qnδM,where δ represents the subsidence coefficient [30, 31], which can be found in the literatures [30, 31], and M is the mining thickness of the seam.


## 3.3.1. Determination ofI1, I2, …, Im−1, Im, …, In−1, In

I
m is used as an example to illustrate the determination of the moment of inertia for different layers of KS. Im is derived using the following equation:(26)Im=112bHm3.


## 3.3.2. Determination ofq1, q2, …, qm−1, qm, …, qn−1, qn

q
m−1 is used as an example to illustrate the respective load for different KS and the load layer weight they control. If the adjacent KS interlayer rock layers are named z1, z2, …, zi−1, zi, with a total of i layers, as shown in Figure 8, then qm−1 is(27)qm−1=γm−1Hm−1+∑i=1iγziHzi.Figure 8
Diagram of adjacent KS.In Figure8, γz1, γz2, …, γzi−1, γzi represent the density of the rock mass z1, z2, …, zi−1, zi, respectively; Ez1, Ez2, …, Ezi−1, Ezi represent the elastic modulus of the rock mass z1, z2, …, zi−1, zi, respectively; Hz1, Hz2, …, Hzi−1, Hzi represent the thickness of the rock mass z1, z2, …, zi−1, zi, respectively. In this study, γ1, γ2, …, γm−1, γm, …, γn−1, γn represent the density of the different layers of KS.


## 3.3.3. Determination ofk1, k2, …, km−1, km, …, kn−1, kn

k
m is used as an example to illustrate the calculation method of k1, k2, …, km−1, km, …, kn−1, kn. The foundation coefficient [15, 16, 20, 30] formed by the total of i layers and KSm−1 shown in Figure 8 is(28)km=1Hm−1/Em−1+∑i=1iHzi/Ezi.In particular, the calculation ofk1 should consider the strata from KS1 to the thick hard rock formation in the coal seam floor [30].


## 3.3.4. Determination ofkcm

Because the model assumes that the fractured rock in the caving zone and fissure zone satisfies Winkler’s elastic foundation hypothesis [15, 16, 20, 30, 32, 33], we can calculate the foundation coefficient kcm of KSm on the goaf side if we know the subsidence value of KSm and the corresponding load bearing of the foundation [30]. When the land surface is under full subsidence [31], the level of subsidence in the curve subsidence zone is consistent [25], and the subsidence of KSm reaches the maximum subsidence [31], δM. At this point, the load of the foundation that corresponds to the maximum subsidence of the side of the goaf, KSm, returns to primary rock stress [7, 11] qm + qm+1 + ⋯ + qn−1 + qn. Hence, the foundation coefficient of the KSm on the side of the goaf, kcm, is(29)kcm=qm+qm+1+…+qn−1+qnδM,where δ represents the subsidence coefficient [30, 31], which can be found in the literatures [30, 31], and M is the mining thickness of the seam.


## 4. Verification and Application of the Method for Calculating Abutment Pressure

### 4.1. Theoretical Calculation

Our estimation method was applied to calculate the distribution of the abutment pressure of working surface 207 in the Tingnan coal mine. The location of the overburden KS and the height of the fissure zone in the Tingnan coal mine are shown in Figure2, based on previous studies [7, 11–13, 17]. There are two KS in both the fissure zone and the curve subsidence zone. In the method used to calculate abutment pressure, n=4 and m = 3, the parameters are shown in Tables 1 and 2, and the results of the calculation are shown in Figure 9. The peak abutment pressure was 40.25 MPa, the stress concentration coefficient was 2.81, the width of the plastic zone was 13.03 m, the abutment pressure width was 283 m, and the peak abutment pressure of the goaf was 11.24 MPa.Table 1
KS parameters.
KS number
E (Gpa)
k (MN/m3)
q (MPa)
H (m)
l (m)
k
c3 (MN/m3)
KS1
E
1 = 8.47
k
1 = 37.1
q
1 = 0.97
H
1 = 12
l
1 = 13.5
—
KS2
E
2 = 19.3
k
2 = 199
q
2 = 1.48
H
2 = 8
l
2 = 17.3
—
KS3
E
3 = 19.3
k
3 = 114
q
3 = 1.03
H
3 = 25.9
—
k
c3 = 2.51
KS4
E
4 = 11.5
k
4 = 366
q
4 = 10.37
H
4 = 222.14
—
—Table 2
Other parameters (apart from KS) in the method used to calculate abutment pressure.
N
0 (MPa)
f
φ (°)
M (m)
δ
q
0 (MPa)
L (m)
6.14
0.29
16
7
0.65
0.45
890Figure 9
Abutment pressure distribution of working surface 207.


### 4.2. Verification

#### 4.2.1. Mining Stress Measurement

It should be noted that it is nearly impossible to measure real in situ stress using stress sensors, and most stress monitoring systems thus record the vertical stress increment relative to the initial pressure of oil pumped into the sensors for coupling with the monitoring boreholes [1, 9, 36–44]. Influenced by the coupling degree between stress meters and boreholes, there may be some disparities among the stress increments measured from different boreholes [1, 36, 37]. The KJ21A real-time stress monitoring system in the Tingnan coal mine was used to monitor the vertical stress increment, which represents the stress variation in the coal. The range of the stress meter is 0–30 MPa, the error is ±0.6 MPa, and the applicable diameter of the borehole is 48–50 mm. More details and applications of this stress monitoring system can be found in [38–41]. It is widely believed in China that the distribution pattern of front abutment pressure is similar to that of side abutment pressure, based on engineering experience and numerous field measurements [9]. Stress meter Nos. 1, 2, and 3 were installed in boreholes at distances of 208 m, 223 m, and 300 m, respectively, in front of the 207 working surface; they were 9 m, 9 m, and 15 m deep, respectively (Figure 1).The initial pressure (P0) of oil pumped into the sensors for coupling with the monitoring boreholes was the initial value. In general, after a few days, the pressure increment value P1 tended to be stable. When the stress meter was at a different distance d from the working face, the pressure value of the stress meter is expressed by P1(d). The measurement results in Figure 10 are the vertical stress increment (P1(d) – P0). If P1(d) – P0 > P1, it is in the increasing stress zone; if P1(d) – P0 < P1, it is in the decreasing stress zone; and if P1(d) – P0 = P1, it is in the initial stress zone. The widths of the increasing and decreasing stress zones can be obtained based on the measurement results. The sum of the widths of the increasing and decreasing stress zones is the abutment pressure width, and the distance between the peak of stress increment (P1(d) – P0) and mining face is the width of the coal plastic region. Comparing these values with the theoretically computed abutment pressure width and the coal plastic region width can verify the validity of the proposed calculation method.Figure 10
Variation curves of vertical stress increment.Figure10 shows that maximum stress occurs when the distance between the working surface and the stress meters are 13 m, 15.4 m, and 12.1 m, respectively, indicating plastic damage in the coal mass between the working surface and stress meters. The widths of the plastic zone obtained by the three stress meters are consistent with the computed results for the width of the plastic zone (13.03 m). In addition, the sums of the monitoring widths of the increasing and decreasing stress zones from stress meter Nos. 1 and 2 were 208 m and 223 m, respectively. Stress meter Nos. 1 and 2 were too close to the working surface and thus could not measure accurate abutment pressure widths, but their measurements indicate that the abutment pressure widths were greater than 208 m and 223 m, respectively. The sum of the monitoring widths of increasing and decreasing stress zones from stress meter No. 3, placed 300 m in front of the working surface, was 267 m, which indicates that the monitoring abutment pressure width was 267 m. Therefore, the results of in situ abutment pressure width monitoring are consistent with the theoretically computed result of a 283 m abutment pressure width. Therefore, the results of stress monitoring are consistent with the calculated results, which show the validity of the proposed prediction method.


#### 4.2.2. Microseismic Monitoring

A correlation is observed between the occurrence of microseismic events and the presence of abutment pressure, which has allowed for examination of the distribution of abutment pressure [9, 42–45]. Many studies [9, 42–45] have determined abutment pressure width through the distribution range of microseismic events, which has been verified by mining practices. The locations of microseismic events relative to working face 207 during the stoping process from April 1 to April 25, 2018, are shown in Figure 11. During this period, the aggregate stoping distance of working face 207 was 103.8 m. The circles in the figure represent the energy ranges of the recorded microseismic events: blue represents the range of 1,000–10,000 J, and red represents events with energy larger than 10,000 J.Figure 11
Locations of microseismic events relative to working face 207.Figure11 shows that microseismic events mainly occurred within a region 276 m lateral to and 283 m in front of working face 207, heralding rock failure in the coal and rocks under this range of abutment pressure [9, 42–45]. Based on the distribution range of microseismic events in working face 207, the along-strike abutment pressure width of working face 207 was close to 283 m and its lateral abutment pressure width was close to 276 m, which agrees with the computed result (283 m) for the abutment pressure width, as shown in Figure 9. This result validates the accuracy of our calculation method.


#### 4.2.3. Numerical Simulations

Based on geological mining conditions, a numerical model was established using UDEC, as shown in Figure12, and the abutment pressure distribution of working surface 207 was examined. The width and height of the model were 2,300 and 570 m, respectively. The rock layer composition of the model was consistent with the integrated stratigraphic column, and the lower boundary, left boundary, and right boundary of the model were fixed. The model was calculated using the Mohr–Coulomb criterion, and the abutment pressure-monitoring line was placed at the lower boundary of the coal seam. Before extraction, the simulation was run to achieve the initial elastic equilibrium and to generate a premining stress field under gravity alone. The original rock stress of the coal seam was 14 MPa, and the mining height was 7.0 m. The positions and dimensions of working face 204, working face 205, and working face 206 are shown in Figure 12, and the UDEC numerical calculation model was calibrated using surface subsidence after the three faces were mined. Table 3 shows the mechanical parameters used in the numerical model.Figure 12
Model geometries.Table 3
Mechanical parameters of rock mass used in the numerical simulation.
Bulk modulus (GPa)
Shear modulus (GPa)
Cohesion (MPa)
Friction angle (°)
Tensile strength (MPa)
Normal stiffness of joint (GPa)
Shear stiffness of joint (GPa)
Loess
3.3
1.5
0.11
10
0
2.7
2.5
Fine sandstone
10.4
7.5
10.9
43.9
11.9
6.2
4.7
Coarse sandstone
8.8
6.0
9.8
42
8.9
5.4
3.8
Mudstone
6.0
3.6
5.1
38
4.7
3.2
2.9
Coarse sandstone gravel
15
8.6
12.4
43.5
13.8
9.5
6.4
Conglomerate
18
11.8
15.4
48
15.5
10.9
8
Medium sandstone
9.2
6.5
10.2
42.5
9.2
5.6
4.3
Sandy mudstone
7.4
4.5
6.2
39
6.6
4.5
3.2
Siltstone
8.5
5.2
9.4
41
8.0
5.2
3.5
Coal
3.5
2.1
3.6
37
1.8
2.5
1.6The results of the UDEC numerical model and the theoretical estimation of the abutment pressure for working face 207 are shown in Figure13. To compare the numerical model results to the theoretical results, we obtained the root mean square errors [46] (0.85 MPa) and relative standard errors [46] (2%) of the theoretical results relative to the numerical results. The root mean square error (0.85 MPa) is much smaller than peak abutment pressure (40.63 MPa) derived from the numerical model, and the relative standard error was less than 5%, which indicates the forms of the stress development curves obtained from the numerical model are in good agreement with those calculated from our proposed method. The numerical model results show a peak abutment pressure of 40.63 MPa, a plastic zone width of 14.2 m, an abutment pressure width of 292 m, and a goaf peak abutment pressure of 11.55 MPa. The theoretical results show a peak abutment pressure of 40.25 MPa, a plastic zone width of 13.03 m, an abutment pressure width of 283 m, and a goaf peak abutment pressure of 11.24 MPa. The theoretical results are in good agreement with the numerical model results.Figure 13
Comparison of stress distributions predicted by UDEC and the proposed method.


### 4.3. Field Application of Computed Abutment Pressure in Determining the Stopping Line Position on the Working Surface and the Width of the Protection Coal Pillar on the Roadway

Determining the stopping position and the width of the protection coal pillar must consider the characteristics of the overburden KS. On the basis of the integrated stratigraphic column in the Tingnan coal mine, the location of the overburden KS and the height of the fissure zone should be calculated based on previous researches [7, 11–13, 17]. Then, the distribution of the abutment pressure of working surface 207 in the Tingnan coal mine can be calculated according to the calculation method described in Section 3. By determining the stopping position and the width of the protection coal pillar, this ensures that the main roadway will be outside the abutment pressure width of working face 207.The results of our calculations show that the abutment pressure width of working surface 207 after excavation is 283 m. If a protection coal pillar on the one side of the stopping line of working surface 207 was 200 m wide, according to the original design, the roadway would be significantly deformed and damaged in a manner similar to what happened at working surface 206. Consequently, working face 207 should stop mining outside of 283 m from the main roadway and leave no less than a 283 m-wide protective coal pillar to ensure that the main roadway is outside the abutment pressure influence width of the working face, as shown in Figure1. As of May 23, 2019, working face 207 is 300 m from the main roadway and only 17 m from the reasonable stop line. The main roadway has not experienced such phenomena as two-side cracking and floor heave, which shows the validity and effectiveness of our proposed method.Setting the KS4 thickness H4 (as given in Table 1) to 135 m, 170 m, and 222.14 m, reflects how changes in KS4 thickness impact the distribution of abutment pressure. The results of these calculations are shown in Figure 14 and show that as the thickness of KS4 increases, the abutment pressure width increases. The thickness of KS4 has a significant influence on the distribution of abutment pressure. When H4 is 222.14 m, 170 m, and 135 m, the abutment pressure widths are 283 m, 225 m, and 178 m, respectively. Therefore, when the KS4 is 135 m thick, maintaining a protection coal pillar 200 m wide protects the main roadway. However, the thickness of KS4 at the Tingnan coal mine is 222.14 m. Although a 200 m-wide protective coal pillar had been set up on the stopping line side of the excavated working surfaces in the Tingnan coal mine, serious deformation and damage still occurred on the main roadway after the excavation of working surface 206. This was due to neglecting the impact of a 222.14 m thick, hard sandstone KS4 in the overburden on the distribution of abutment pressure, which resulted in a significant increase in the actual abutment pressure width compared to Qian et al.’s [7] or Peng and Chiang’s [21] empirical value.Figure 14
Abutment pressure distributions under different KS4 thicknesses.


## 4.1. Theoretical Calculation

Our estimation method was applied to calculate the distribution of the abutment pressure of working surface 207 in the Tingnan coal mine. The location of the overburden KS and the height of the fissure zone in the Tingnan coal mine are shown in Figure2, based on previous studies [7, 11–13, 17]. There are two KS in both the fissure zone and the curve subsidence zone. In the method used to calculate abutment pressure, n=4 and m = 3, the parameters are shown in Tables 1 and 2, and the results of the calculation are shown in Figure 9. The peak abutment pressure was 40.25 MPa, the stress concentration coefficient was 2.81, the width of the plastic zone was 13.03 m, the abutment pressure width was 283 m, and the peak abutment pressure of the goaf was 11.24 MPa.Table 1
KS parameters.
KS number
E (Gpa)
k (MN/m3)
q (MPa)
H (m)
l (m)
k
c3 (MN/m3)
KS1
E
1 = 8.47
k
1 = 37.1
q
1 = 0.97
H
1 = 12
l
1 = 13.5
—
KS2
E
2 = 19.3
k
2 = 199
q
2 = 1.48
H
2 = 8
l
2 = 17.3
—
KS3
E
3 = 19.3
k
3 = 114
q
3 = 1.03
H
3 = 25.9
—
k
c3 = 2.51
KS4
E
4 = 11.5
k
4 = 366
q
4 = 10.37
H
4 = 222.14
—
—Table 2
Other parameters (apart from KS) in the method used to calculate abutment pressure.
N
0 (MPa)
f
φ (°)
M (m)
δ
q
0 (MPa)
L (m)
6.14
0.29
16
7
0.65
0.45
890Figure 9
Abutment pressure distribution of working surface 207.


## 4.2. Verification

### 4.2.1. Mining Stress Measurement

It should be noted that it is nearly impossible to measure real in situ stress using stress sensors, and most stress monitoring systems thus record the vertical stress increment relative to the initial pressure of oil pumped into the sensors for coupling with the monitoring boreholes [1, 9, 36–44]. Influenced by the coupling degree between stress meters and boreholes, there may be some disparities among the stress increments measured from different boreholes [1, 36, 37]. The KJ21A real-time stress monitoring system in the Tingnan coal mine was used to monitor the vertical stress increment, which represents the stress variation in the coal. The range of the stress meter is 0–30 MPa, the error is ±0.6 MPa, and the applicable diameter of the borehole is 48–50 mm. More details and applications of this stress monitoring system can be found in [38–41]. It is widely believed in China that the distribution pattern of front abutment pressure is similar to that of side abutment pressure, based on engineering experience and numerous field measurements [9]. Stress meter Nos. 1, 2, and 3 were installed in boreholes at distances of 208 m, 223 m, and 300 m, respectively, in front of the 207 working surface; they were 9 m, 9 m, and 15 m deep, respectively (Figure 1).The initial pressure (P0) of oil pumped into the sensors for coupling with the monitoring boreholes was the initial value. In general, after a few days, the pressure increment value P1 tended to be stable. When the stress meter was at a different distance d from the working face, the pressure value of the stress meter is expressed by P1(d). The measurement results in Figure 10 are the vertical stress increment (P1(d) – P0). If P1(d) – P0 > P1, it is in the increasing stress zone; if P1(d) – P0 < P1, it is in the decreasing stress zone; and if P1(d) – P0 = P1, it is in the initial stress zone. The widths of the increasing and decreasing stress zones can be obtained based on the measurement results. The sum of the widths of the increasing and decreasing stress zones is the abutment pressure width, and the distance between the peak of stress increment (P1(d) – P0) and mining face is the width of the coal plastic region. Comparing these values with the theoretically computed abutment pressure width and the coal plastic region width can verify the validity of the proposed calculation method.Figure 10
Variation curves of vertical stress increment.Figure10 shows that maximum stress occurs when the distance between the working surface and the stress meters are 13 m, 15.4 m, and 12.1 m, respectively, indicating plastic damage in the coal mass between the working surface and stress meters. The widths of the plastic zone obtained by the three stress meters are consistent with the computed results for the width of the plastic zone (13.03 m). In addition, the sums of the monitoring widths of the increasing and decreasing stress zones from stress meter Nos. 1 and 2 were 208 m and 223 m, respectively. Stress meter Nos. 1 and 2 were too close to the working surface and thus could not measure accurate abutment pressure widths, but their measurements indicate that the abutment pressure widths were greater than 208 m and 223 m, respectively. The sum of the monitoring widths of increasing and decreasing stress zones from stress meter No. 3, placed 300 m in front of the working surface, was 267 m, which indicates that the monitoring abutment pressure width was 267 m. Therefore, the results of in situ abutment pressure width monitoring are consistent with the theoretically computed result of a 283 m abutment pressure width. Therefore, the results of stress monitoring are consistent with the calculated results, which show the validity of the proposed prediction method.


### 4.2.2. Microseismic Monitoring

A correlation is observed between the occurrence of microseismic events and the presence of abutment pressure, which has allowed for examination of the distribution of abutment pressure [9, 42–45]. Many studies [9, 42–45] have determined abutment pressure width through the distribution range of microseismic events, which has been verified by mining practices. The locations of microseismic events relative to working face 207 during the stoping process from April 1 to April 25, 2018, are shown in Figure 11. During this period, the aggregate stoping distance of working face 207 was 103.8 m. The circles in the figure represent the energy ranges of the recorded microseismic events: blue represents the range of 1,000–10,000 J, and red represents events with energy larger than 10,000 J.Figure 11
Locations of microseismic events relative to working face 207.Figure11 shows that microseismic events mainly occurred within a region 276 m lateral to and 283 m in front of working face 207, heralding rock failure in the coal and rocks under this range of abutment pressure [9, 42–45]. Based on the distribution range of microseismic events in working face 207, the along-strike abutment pressure width of working face 207 was close to 283 m and its lateral abutment pressure width was close to 276 m, which agrees with the computed result (283 m) for the abutment pressure width, as shown in Figure 9. This result validates the accuracy of our calculation method.


### 4.2.3. Numerical Simulations

Based on geological mining conditions, a numerical model was established using UDEC, as shown in Figure12, and the abutment pressure distribution of working surface 207 was examined. The width and height of the model were 2,300 and 570 m, respectively. The rock layer composition of the model was consistent with the integrated stratigraphic column, and the lower boundary, left boundary, and right boundary of the model were fixed. The model was calculated using the Mohr–Coulomb criterion, and the abutment pressure-monitoring line was placed at the lower boundary of the coal seam. Before extraction, the simulation was run to achieve the initial elastic equilibrium and to generate a premining stress field under gravity alone. The original rock stress of the coal seam was 14 MPa, and the mining height was 7.0 m. The positions and dimensions of working face 204, working face 205, and working face 206 are shown in Figure 12, and the UDEC numerical calculation model was calibrated using surface subsidence after the three faces were mined. Table 3 shows the mechanical parameters used in the numerical model.Figure 12
Model geometries.Table 3
Mechanical parameters of rock mass used in the numerical simulation.
Bulk modulus (GPa)
Shear modulus (GPa)
Cohesion (MPa)
Friction angle (°)
Tensile strength (MPa)
Normal stiffness of joint (GPa)
Shear stiffness of joint (GPa)
Loess
3.3
1.5
0.11
10
0
2.7
2.5
Fine sandstone
10.4
7.5
10.9
43.9
11.9
6.2
4.7
Coarse sandstone
8.8
6.0
9.8
42
8.9
5.4
3.8
Mudstone
6.0
3.6
5.1
38
4.7
3.2
2.9
Coarse sandstone gravel
15
8.6
12.4
43.5
13.8
9.5
6.4
Conglomerate
18
11.8
15.4
48
15.5
10.9
8
Medium sandstone
9.2
6.5
10.2
42.5
9.2
5.6
4.3
Sandy mudstone
7.4
4.5
6.2
39
6.6
4.5
3.2
Siltstone
8.5
5.2
9.4
41
8.0
5.2
3.5
Coal
3.5
2.1
3.6
37
1.8
2.5
1.6The results of the UDEC numerical model and the theoretical estimation of the abutment pressure for working face 207 are shown in Figure13. To compare the numerical model results to the theoretical results, we obtained the root mean square errors [46] (0.85 MPa) and relative standard errors [46] (2%) of the theoretical results relative to the numerical results. The root mean square error (0.85 MPa) is much smaller than peak abutment pressure (40.63 MPa) derived from the numerical model, and the relative standard error was less than 5%, which indicates the forms of the stress development curves obtained from the numerical model are in good agreement with those calculated from our proposed method. The numerical model results show a peak abutment pressure of 40.63 MPa, a plastic zone width of 14.2 m, an abutment pressure width of 292 m, and a goaf peak abutment pressure of 11.55 MPa. The theoretical results show a peak abutment pressure of 40.25 MPa, a plastic zone width of 13.03 m, an abutment pressure width of 283 m, and a goaf peak abutment pressure of 11.24 MPa. The theoretical results are in good agreement with the numerical model results.Figure 13
Comparison of stress distributions predicted by UDEC and the proposed method.


## 4.2.1. Mining Stress Measurement

It should be noted that it is nearly impossible to measure real in situ stress using stress sensors, and most stress monitoring systems thus record the vertical stress increment relative to the initial pressure of oil pumped into the sensors for coupling with the monitoring boreholes [1, 9, 36–44]. Influenced by the coupling degree between stress meters and boreholes, there may be some disparities among the stress increments measured from different boreholes [1, 36, 37]. The KJ21A real-time stress monitoring system in the Tingnan coal mine was used to monitor the vertical stress increment, which represents the stress variation in the coal. The range of the stress meter is 0–30 MPa, the error is ±0.6 MPa, and the applicable diameter of the borehole is 48–50 mm. More details and applications of this stress monitoring system can be found in [38–41]. It is widely believed in China that the distribution pattern of front abutment pressure is similar to that of side abutment pressure, based on engineering experience and numerous field measurements [9]. Stress meter Nos. 1, 2, and 3 were installed in boreholes at distances of 208 m, 223 m, and 300 m, respectively, in front of the 207 working surface; they were 9 m, 9 m, and 15 m deep, respectively (Figure 1).The initial pressure (P0) of oil pumped into the sensors for coupling with the monitoring boreholes was the initial value. In general, after a few days, the pressure increment value P1 tended to be stable. When the stress meter was at a different distance d from the working face, the pressure value of the stress meter is expressed by P1(d). The measurement results in Figure 10 are the vertical stress increment (P1(d) – P0). If P1(d) – P0 > P1, it is in the increasing stress zone; if P1(d) – P0 < P1, it is in the decreasing stress zone; and if P1(d) – P0 = P1, it is in the initial stress zone. The widths of the increasing and decreasing stress zones can be obtained based on the measurement results. The sum of the widths of the increasing and decreasing stress zones is the abutment pressure width, and the distance between the peak of stress increment (P1(d) – P0) and mining face is the width of the coal plastic region. Comparing these values with the theoretically computed abutment pressure width and the coal plastic region width can verify the validity of the proposed calculation method.Figure 10
Variation curves of vertical stress increment.Figure10 shows that maximum stress occurs when the distance between the working surface and the stress meters are 13 m, 15.4 m, and 12.1 m, respectively, indicating plastic damage in the coal mass between the working surface and stress meters. The widths of the plastic zone obtained by the three stress meters are consistent with the computed results for the width of the plastic zone (13.03 m). In addition, the sums of the monitoring widths of the increasing and decreasing stress zones from stress meter Nos. 1 and 2 were 208 m and 223 m, respectively. Stress meter Nos. 1 and 2 were too close to the working surface and thus could not measure accurate abutment pressure widths, but their measurements indicate that the abutment pressure widths were greater than 208 m and 223 m, respectively. The sum of the monitoring widths of increasing and decreasing stress zones from stress meter No. 3, placed 300 m in front of the working surface, was 267 m, which indicates that the monitoring abutment pressure width was 267 m. Therefore, the results of in situ abutment pressure width monitoring are consistent with the theoretically computed result of a 283 m abutment pressure width. Therefore, the results of stress monitoring are consistent with the calculated results, which show the validity of the proposed prediction method.


## 4.2.2. Microseismic Monitoring

A correlation is observed between the occurrence of microseismic events and the presence of abutment pressure, which has allowed for examination of the distribution of abutment pressure [9, 42–45]. Many studies [9, 42–45] have determined abutment pressure width through the distribution range of microseismic events, which has been verified by mining practices. The locations of microseismic events relative to working face 207 during the stoping process from April 1 to April 25, 2018, are shown in Figure 11. During this period, the aggregate stoping distance of working face 207 was 103.8 m. The circles in the figure represent the energy ranges of the recorded microseismic events: blue represents the range of 1,000–10,000 J, and red represents events with energy larger than 10,000 J.Figure 11
Locations of microseismic events relative to working face 207.Figure11 shows that microseismic events mainly occurred within a region 276 m lateral to and 283 m in front of working face 207, heralding rock failure in the coal and rocks under this range of abutment pressure [9, 42–45]. Based on the distribution range of microseismic events in working face 207, the along-strike abutment pressure width of working face 207 was close to 283 m and its lateral abutment pressure width was close to 276 m, which agrees with the computed result (283 m) for the abutment pressure width, as shown in Figure 9. This result validates the accuracy of our calculation method.


## 4.2.3. Numerical Simulations

Based on geological mining conditions, a numerical model was established using UDEC, as shown in Figure12, and the abutment pressure distribution of working surface 207 was examined. The width and height of the model were 2,300 and 570 m, respectively. The rock layer composition of the model was consistent with the integrated stratigraphic column, and the lower boundary, left boundary, and right boundary of the model were fixed. The model was calculated using the Mohr–Coulomb criterion, and the abutment pressure-monitoring line was placed at the lower boundary of the coal seam. Before extraction, the simulation was run to achieve the initial elastic equilibrium and to generate a premining stress field under gravity alone. The original rock stress of the coal seam was 14 MPa, and the mining height was 7.0 m. The positions and dimensions of working face 204, working face 205, and working face 206 are shown in Figure 12, and the UDEC numerical calculation model was calibrated using surface subsidence after the three faces were mined. Table 3 shows the mechanical parameters used in the numerical model.Figure 12
Model geometries.Table 3
Mechanical parameters of rock mass used in the numerical simulation.
Bulk modulus (GPa)
Shear modulus (GPa)
Cohesion (MPa)
Friction angle (°)
Tensile strength (MPa)
Normal stiffness of joint (GPa)
Shear stiffness of joint (GPa)
Loess
3.3
1.5
0.11
10
0
2.7
2.5
Fine sandstone
10.4
7.5
10.9
43.9
11.9
6.2
4.7
Coarse sandstone
8.8
6.0
9.8
42
8.9
5.4
3.8
Mudstone
6.0
3.6
5.1
38
4.7
3.2
2.9
Coarse sandstone gravel
15
8.6
12.4
43.5
13.8
9.5
6.4
Conglomerate
18
11.8
15.4
48
15.5
10.9
8
Medium sandstone
9.2
6.5
10.2
42.5
9.2
5.6
4.3
Sandy mudstone
7.4
4.5
6.2
39
6.6
4.5
3.2
Siltstone
8.5
5.2
9.4
41
8.0
5.2
3.5
Coal
3.5
2.1
3.6
37
1.8
2.5
1.6The results of the UDEC numerical model and the theoretical estimation of the abutment pressure for working face 207 are shown in Figure13. To compare the numerical model results to the theoretical results, we obtained the root mean square errors [46] (0.85 MPa) and relative standard errors [46] (2%) of the theoretical results relative to the numerical results. The root mean square error (0.85 MPa) is much smaller than peak abutment pressure (40.63 MPa) derived from the numerical model, and the relative standard error was less than 5%, which indicates the forms of the stress development curves obtained from the numerical model are in good agreement with those calculated from our proposed method. The numerical model results show a peak abutment pressure of 40.63 MPa, a plastic zone width of 14.2 m, an abutment pressure width of 292 m, and a goaf peak abutment pressure of 11.55 MPa. The theoretical results show a peak abutment pressure of 40.25 MPa, a plastic zone width of 13.03 m, an abutment pressure width of 283 m, and a goaf peak abutment pressure of 11.24 MPa. The theoretical results are in good agreement with the numerical model results.Figure 13
Comparison of stress distributions predicted by UDEC and the proposed method.


## 4.3. Field Application of Computed Abutment Pressure in Determining the Stopping Line Position on the Working Surface and the Width of the Protection Coal Pillar on the Roadway

Determining the stopping position and the width of the protection coal pillar must consider the characteristics of the overburden KS. On the basis of the integrated stratigraphic column in the Tingnan coal mine, the location of the overburden KS and the height of the fissure zone should be calculated based on previous researches [7, 11–13, 17]. Then, the distribution of the abutment pressure of working surface 207 in the Tingnan coal mine can be calculated according to the calculation method described in Section 3. By determining the stopping position and the width of the protection coal pillar, this ensures that the main roadway will be outside the abutment pressure width of working face 207.The results of our calculations show that the abutment pressure width of working surface 207 after excavation is 283 m. If a protection coal pillar on the one side of the stopping line of working surface 207 was 200 m wide, according to the original design, the roadway would be significantly deformed and damaged in a manner similar to what happened at working surface 206. Consequently, working face 207 should stop mining outside of 283 m from the main roadway and leave no less than a 283 m-wide protective coal pillar to ensure that the main roadway is outside the abutment pressure influence width of the working face, as shown in Figure1. As of May 23, 2019, working face 207 is 300 m from the main roadway and only 17 m from the reasonable stop line. The main roadway has not experienced such phenomena as two-side cracking and floor heave, which shows the validity and effectiveness of our proposed method.Setting the KS4 thickness H4 (as given in Table 1) to 135 m, 170 m, and 222.14 m, reflects how changes in KS4 thickness impact the distribution of abutment pressure. The results of these calculations are shown in Figure 14 and show that as the thickness of KS4 increases, the abutment pressure width increases. The thickness of KS4 has a significant influence on the distribution of abutment pressure. When H4 is 222.14 m, 170 m, and 135 m, the abutment pressure widths are 283 m, 225 m, and 178 m, respectively. Therefore, when the KS4 is 135 m thick, maintaining a protection coal pillar 200 m wide protects the main roadway. However, the thickness of KS4 at the Tingnan coal mine is 222.14 m. Although a 200 m-wide protective coal pillar had been set up on the stopping line side of the excavated working surfaces in the Tingnan coal mine, serious deformation and damage still occurred on the main roadway after the excavation of working surface 206. This was due to neglecting the impact of a 222.14 m thick, hard sandstone KS4 in the overburden on the distribution of abutment pressure, which resulted in a significant increase in the actual abutment pressure width compared to Qian et al.’s [7] or Peng and Chiang’s [21] empirical value.Figure 14
Abutment pressure distributions under different KS4 thicknesses.


## 5. Conclusions

Abutment pressure distribution is closely related to the overburden KS. Neglecting the influence of the overburden KS on the distribution of abutment pressure often leads to coal mine hazards, such as rockbursts, and so on. However, current calculation methods for working surface abutment pressure do not consider the influence of the overburden KS. Using KS theory and Winkler’s elastic foundation beam theory, we consider the fissure zone KS on the coal mass side and the curve subsidence zone KS as many elastic foundation beams interact with each other. A method to calculate the abutment pressure of the coal mass and the goaf was then established, considering the influence of the overburden KS. The abutment pressure distribution of working surface 207 in the Tingnan coal mine was then calculated using our method and was verified using measurements from borehole stress meters and microseismic monitoring systems, as well as numerical simulations.Results of using the proposed calculation method show that as the thickness of the KS4 in the Tingnan coal mine increased, the abutment pressure width significantly increased. The working faces 205 and 206 neglected the impact of 222.14 m thick, hard sandstone KS4 in the overburden on the distribution of abutment pressure, and only a 200 m-wide protective coal pillar was left, causing serious deformation and damage to the roadway. Based on the results that show the abutment pressure width to be 283 m for working surface 207, the working face must stop 283 m from the main roadway and leave a, no less than 283 m wide, protection coal pillar.The abutment pressure calculation method established in this paper can calculate the abutment pressure distribution of a working face under different overburden KS conditions, providing guidance for forecasting and preventing mine dynamic hazards, controlling the surrounding rock in mining roadways, and designing the layout of mining roadways. However, when the overburden strata have a thick alluvium layer, the thick alluvium layer will also have a strong impact on the abutment pressure. It is necessary to further study the calculation method of abutment pressure under the condition that both KS and thick alluvium layers exist in the overburden.



---

*Source: 7678327-2019-07-15.xml*