Study on Size Design of Shaft Protection Rock/Coal Pillars in Thick Soil and Thin Rock Strata

Abstract

To prevent serious shaft deflection disasters under asymmetric mining conditions, it is urgent to solve the problem of designing shaft protection rock pillar (SPRP) sizes in thick soil and thin rock strata. In this paper, based on the parallel mining model and the perpendicular mining model, a dynamic prediction model that can describe the horizontal movement of the shaft was established by the probability integration method and the Knothe time function. Next, according to the measured data of the shaft deflection in the Guotun Coal Mine, a back analysis was used to calculate the prediction parameters that were suitable for the deep soil strata. Based on the mining model, the variation law of the horizontal deflection displacement of the shaft and SPRP size was obtained. The results showed that the final displacements of the shaft under the two ideal mining models were equal, while the parallel mining model was superior to the perpendicular mining model at the initial stage of mining. The horizontal displacement of the shaft head had a nonlinear negative correlation with the SPRP, and the SPRP size in thick soil and thin rock strata calculated by the parallel mining model was more reasonable. For the Guotun Coal Mine, when the soil movement angle was 57.8% of the actual value, the horizontal displacement of the main shaft head was reduced by 87%. The results have important theoretical and practical value in preventing shaft deflection in thick soil and thin rock strata.

1. Introduction

Since 2002, a large number of vertical shafts have been built in deep soil strata in China. There are 71, 28, 4, and 3 shafts with soil thicknesses greater than 400, 500, 600, and 700 m, respectively, and the maximum soil thickness is 754.98 m (the main shaft of the Wanfu Coal Mine in Shandong Province). Additionally, vertical shafts with soil thicknesses exceeding half the total depth account for 86% of the 71 shafts (see Figure 1). As the shaft is the “throat” of a coal mine, its safety is vital to the survival of the entire mine. However, underground mining activities may cause the shaft to deflect and subside, which significantly affects hoisting safety and poses a great threat to lining safety. To protect the shaft from the effects of mining activities, a shaft protection coal pillar (SPCP) or shaft protection rock pillar (SPRP) must be retained. The design methods for SPCPs are basically similar, and existing SPCPs are mostly designed by experience or angles. Due to the large size of SPCPs designed by empirical methods, the boundary angle or the movement angle has been used to design SPCPs in European coal mines [1] and is stipulated in existing criteria [2] and manuals [3] in China. Much research has been carried out on the design of SPCPs. Zeenke et al. studied the shape and size of SPCPs in the Ostrava–Karvin Coal Mine and discussed the combination of the SPCP boundary delineated by the pillar angle under different dip angles [4]. Haupt et al. calculated the SPCP size under single or multiple seam mining conditions using the Elhart–Shaoer method with predetermined boundary values [5]. Kratzsch believed that an insufficient SPCP size is the main cause of shaft damage during mining and proposed some measures based on different mining conditions and lining types, such as setting a sliding layer, separating seam, and compressible wood cushion in the shaft lining [6]. Considering the predicted shaft deformation insurance coefficient, Sroka determined the movement angles by using the critical deformation value and obtained a formula for calculating the SPCP based on observation data of SPCP mining in Germany and Poland [7]. Wei et al. studied the influence factors and change law of SPCP size and obtained a more realistic nonlinear formula for calculating SPCP size [8]. In addition, based on the vertical section method, Wei et al. analyzed the shortcomings of the existing graphic method and established an analytical mathematical model of designing SPCPs [9]. With comprehensive analysis of a large number of field observation data and mining slip mechanisms and sliding vectors, He et al. established a mathematical model for estimating the SPCP that is suitable for mountainous areas [10]. However, the above research was conducted under thin soil strata conditions. In addition, the existing movement angles both in China and internationally are the observations in the case of thin soil strata [11]. As seen in

Table 1. The movement angles of some mine areas in China.

Name of Mine	Mining Depth (m)	Movement Angles (°)
		Rock Strata				Soil Strata
		Downhill	Uphill	Strike	Dip
Jiaohe	35–110	75‒0.8α	75	75		45
Xuzhou	90–140	75−0.82α; 70 −0.72α	75; 70	75; 70		45; 36
Shuangyashan	30–220	75 − 0.3α	68	70		45
Pingyuan	100–330	72	72	68		55
Huainan	<180	75 − 0.65α; 53 − 0.1α	75	75	55; 75 − 0.65α; 53 − 0.1α	40–45
Fengfeng	<260	73 − 0.6α	73	73		58
Weizhou	<310	73 − 0.6α	73	73		45–55
Fuxin	<400	73; 83 − 0.9α	75	72		40–50
Fushun	<540	59 − 0.2α	62	65		45
Kailuan	<600	72 − 0.67α (≥30)	35–72	70		35–45
Zaozhuang	<600	86.6 − α	76	76		45
Changzhi	<600	68 − 73	71–74	71–74		55–66
Jining	>600	65 − 75	75	75		40–45
Juye	>600	75 ‒ 0.3α	70–75	70–75		40–45

Note: α is the dip angle of coal seam.

, for different mining conditions, the movement angles of some mine areas in China are basically unchanged. The result is that the SPRP size designed by existing movement angles is obviously insufficient in deep soil strata, resulting in the shaft being within the range of mining influence, and mining activity causes the soil around the shaft at different depths to produce different displacement values. Under asymmetric mining conditions, the shaft will be deflected as same as that in the Guotun Coal Mine. To prevent serious shaft deflection disasters, it is urgent to solve the problem of designing SPRP sizes in thick soil and thin rock strata.

We addressed this issue in this study by, first, establishing two more ideal mining models than actual mining according to the general mining law. Based on this, a dynamic prediction model that can describe the horizontal movement of the shaft was obtained by combining the probability integral method with the Knothe time function. Second, for the main shaft deflection in the Guotun Coal Mine, a back analysis was adopted to calculate prediction parameters that are suitable for deep soil strata. Finally, the relationship between the SPRP size and the deflection displacement of the main shaft was obtained by using the ideal mining model. By analyzing their relationships, a more reasonable SPRP size than that from the specification design was obtained. The results have important guiding significance for protecting lining safety and preventing shaft deflection disasters, and they can also provide a scientific and theoretical basis for mining SPCPs.

2. Project Case

The main shaft, the auxiliary shaft, and the air shaft of the Guotun Coal Mine in Shandong Province are located to the side of the coal-free area. The design parameters are as follows: The width of the enclosure belt is 20 m, and the movement angles of the soil and the rock are 45° and 70°, respectively. Therefore, the minimum distance from the main shaft to the coal pillars is approximately 858 m. The auxiliary shaft is 92 m ES49° away from the main shaft. The shafts in the Guotun Coal Mine have unexpectedly experienced serious deflections within the soil section (see

Table 2. Design parameters and deflective situations of shafts in the Guotun Coal Mine.

Parameters	Main Shaft	Auxiliary Shaft
Shaft depth (m)	853	882
Soil thickness (m)	587.4	582.7
Date	$\frac{\mathrm{Maximum} \mathrm{deflection} \mathrm{value} \left(\mathrm{mm}\right)}{\mathrm{Azimuth} \mathrm{angle} (°)}$
July 2015	$\frac{349}{95}$	$\frac{322}{77}$
August 2017	$\frac{359}{79}$	$\frac{318}{75}$

Note: The maximum deflection value in the table refers to the measured displacement of the main (auxiliary) shaft head relative to the shaft depth at the intersection of the soil and rock layers [12].

). This proves that the movement angles specified in the current code are not suitable for the design of SPRPs in thick soil and thin rock strata.

3. Models and Methods

3.1. Calculation Methods

As early as the 1950s, the stochastic medium theory was first introduced by the Polish scholar Litwiniszyn for the study of rock movement [13]. The method was further developed into a probability integral method by the Chinese scholar Liu et al. [14]. Until now, the probability integral method has been the most widely used method for predicting mining subsidence [15,16,17]. The above prediction methods can only predict the ultimate deformation. However, with the mining of the working face, the stratum structures are inevitably subjected to dynamic deformations associated with the formation of a surface subsidence basin [18]. Currently, several models for predicting dynamic subsidence of the surface have been developed. The commonly used models are as follows: (1) Time function models were established by Knothe based on the Mitscherlich growth law [19], which has been widely used by many scholars around the world [20,21,22]. (2) Empirical models are based on the measurement of surface deformation [23,24]. (3) Models proposed by Peng and Luo and further improved by Luo are based on the normal distribution [18,25]. (4) Models deriving surface subsidence over time are based on the rheological characteristics of overburden strata [26]. Significantly, to determine the time function of a mining area’s subsidence, the corresponding parameters of the time function for the mining area must be accurately obtained. Considering that the above dynamic prediction models describe the movement of one point on the surface, the displacement of arbitrary points in the strata can be analyzed by using the probability integral method.

3.2. Mining Models

Based on the requirements of the problem, two ideal mining models —the parallel mining model (see Figure 2a) and the perpendicular mining model (see Figure 2b) —were established according to the relationship between the mining direction and the central line of the mining area. The mining area is completely mined, and the shaft is located on the central line. An independent coordinate system (ξCη) and a rectangular mining area (expressed by I and II, respectively) are arranged along the north and south sides of the central line. Each mining area is composed of n identical horizontal working faces and is mined in order from small to large according to the face number. To ensure the deflection direction of the shaft remains constant during the mining process, the same face numbers are selected for simultaneous mining. Assume that the strike length and the tendency length of the working face are l and s, respectively; the average advancing rate along the strike is v; the mining depth is H; and the minimum SPRP size is B. The coordinates on both sides of the central line are calculated according to the same-side coordinate system. For the sake of convenient analysis, the superscripts I and II were adopted to indicate the parameters of the corresponding mining area.

3.3. Model Solution

First, taking the parallel mining model as an example, when only mining area I is mined, a surface coordinate system is established and coincides with the horizontal projection of the face coordinate system shown in Figure 3. Assuming that the coordinates of unit E in the i-th working face relative to the face coordinate system are (${\xi }_E^{\mathrm{I}}$, ${\eta }_E^{\mathrm{I}}$, $H^{\mathrm{I}}$), the coordinates of point A in the strata relative to the surface coordinate system are ($x_A^{\mathrm{I}}$, $y_A^{\mathrm{I}}$, $z_A^{\mathrm{I}}$). The points A’ and E’ are the horizontal projections of A and E at the surface, respectively. Therefore, the ultimate subsidence displacement of point A caused by the mining of unit E can be calculated by the probability integral theory. In addition, considering the spatial and temporal characteristics of the mining surface subsidence, Knothe obtained a dynamic time function that is related only to the lithological time coefficient c by assuming that the subsidence velocity is directly proportional to the difference between the maximum subsidence value and the dynamic subsidence value at a certain moment. The dynamic time function T(t) can be expressed as

$$T\left(t\right)=1-{\mathrm{e}}^{-ct}.$$

(1)

Since any point in the subsidence strata has a similar movement law to the surface subsidence in time, under the same conditions, the movement laws of all points can be approximately expressed by the same time function. If t_i denotes the time of the i-th working face from the initial mining to a certain time, it can be assumed to be an integral multiple of l/v to enable convenient calculations. Thus, the subsidence value of point A caused by mining the i-th working face at a certain moment can be derived as

$$w_{Ai}^{\mathrm{I}}=w_0\left(1-{\mathrm{e}}^{-c{t}_i}\right){\displaystyle {\int }_0^l{\displaystyle {\int }_{\left(n-i\right)s}^{\left(n-i+1\right)s}\frac{1}{r_{z_A}^{\mathrm{I}}{}^2}{\mathrm{e}}^{-\mathsf{\pi }\frac{{(x_A^{\mathrm{I}}-{\xi }_E^{\mathrm{I}})}^{{}_2}+{(y_A^{\mathrm{I}}-{\eta }_E^{\mathrm{I}})}^2}{r_{z_A}^{\mathrm{I}}{}^2}}}}\mathrm{d}\eta \mathrm{d}\xi ,$$

(2)

where $r_{z_A}^{\mathrm{I}}$ is the main influence radius at depth $z_A^{\mathrm{I}}$.

The relationship between the main influence radius $r_{z_A}^{\mathrm{I}}$, the mining depth H^I, and the main influence angle β is shown in Figure 4. Based on the probability integral theory, the expression of $r_{z_A}^{\mathrm{I}}$ can be defined as

$$r_{z_A}^{\mathrm{I}}=r_0{\left(\frac{H^{\mathrm{I}}-z_A^{\mathrm{I}}}{H^{\mathrm{I}}}\right)}^{\lambda }=\frac{H^{\mathrm{I}}}{\tan \beta }{\left(\frac{H^{\mathrm{I}}-z_A^{\mathrm{I}}}{H^{\mathrm{I}}}\right)}^{\lambda },$$

(3)

where r₀ is the main influence radius at the surface and λ is a constant.

Based on the superposition principle of displacement, the subsidence value of point A caused by mining n working faces at the same moment can be written as

$$w_A^{\mathrm{I}}={\displaystyle \sum _{i=1}^n{w}_0\left(1-{\mathrm{e}}^{-c{t}_i}\right){\displaystyle {\int }_0^l{\displaystyle {\int }_{\left(n-i\right)s}^{\left(n-i+1\right)s}\frac{1}{r_{z_A}^{\mathrm{I}}{}^2}{\mathrm{e}}^{-\mathsf{\pi }\frac{{(x_A^{\mathrm{I}}-{\xi }_E^{\mathrm{I}})}^{{}_2}+{(y_A^{\mathrm{I}}-{\eta }_E^{\mathrm{I}})}^2}{r_{z_A}^{\mathrm{I}}{}^2}}}}\mathrm{d}\eta \mathrm{d}\xi }.$$

(4)

After integration, Equation (4) can be rewritten as

$$w_A^{\mathrm{I}}={\displaystyle \sum _{i=1}^n\frac{w_0\left(1-{\mathrm{e}}^{-c{t}_i}\right)}{4}\left\{\mathrm{erf}\left(k_z^{\mathrm{I}}{x}_A^{\mathrm{I}}\right)-\mathrm{erf}\left[k_z^{\mathrm{I}}\left(x_A^{\mathrm{I}}-l\right)\right]\right\}}\left\{\mathrm{erf}\left[k_z^{\mathrm{I}}\left[y_A^{\mathrm{I}}-\left(n-i\right)s\right]\right]-\mathrm{erf}\left[k_z^{\mathrm{I}}\left[y_A^{\mathrm{I}}-\left(n-i+1\right)s\right]\right]\right\},$$

(5)

where w₀ is the maximum subsidence value of the strata when the working face is fully mined along the strike, and dip and erf represent the probability integral function, which can be defined by

$$w_0=m\cdot q, \mathrm{erf}\left(k_z^{\mathrm{I}}{x}_A^{\mathrm{I}}\right)=\frac{2}{\sqrt{\mathsf{\pi }}}{\displaystyle {\int }_0^{k_z^{\mathrm{I}}{x}_A^{\mathrm{I}}}{e}^{-u^2}\mathrm{d}u}, k_z^{\mathrm{I}}=\frac{\sqrt{\mathsf{\pi }}}{r_{z_A}^{\mathrm{I}}},$$

(6)

where m and q are the mining thickness and the coefficient of stratal subsidence, respectively.

Suppose that s₀ is the displacement of the inflection point and l_c is the goaf critical size. According to the probability integral theory, l_c = 2r₀ + 2s₀. It is well known that when the opening reaches its critical size along the strike and dip, the maximum value of the surface subsidence is approximately equal to 0.98w₀ [27]. Therefore, when the working face is fully mined along only the strike, the maximum value of surface subsidence can be calculated to be equal to $0.98w_0\mathrm{erf}[\sqrt{\mathsf{\pi }}s/\left(2r_0\right)]$. The advancing time is l_c/v in the critical mining state, and when z = 0, the lithology time coefficient is equal to

$$c=-\frac{v}{2r_0+2s_0}\ln 0.02.$$

(7)

To study the horizontal movement of the strata and its deformation caused by coal mining, φ is defined as the angle rotated counterclockwise from the x axis in the positive direction to the specified direction at point A’. Thus, the inclination ${\chi }_A^{\mathrm{I}}$ of point A in the φ^I direction can be expressed as the derivative of $w_A^{\mathrm{I}}$ in that direction with the unit length, which can be expressed as

$${\chi }_A^{\mathrm{I}}={\displaystyle \sum _{i=1}^n\left(\cos {\phi }^{\mathrm{I}}\frac{\partial w_A^{\mathrm{I}}}{\partial x}+\sin {\phi }^{\mathrm{I}}\frac{\partial w_A^{\mathrm{I}}}{\partial y}\right)}.$$

(8)

In [28], the relationship between the horizontal displacement, $u_A^{\mathrm{I}}$, and ${\chi }_A^{\mathrm{I}}$ at point A is as follows:

$$u_A^{\mathrm{I}}={\displaystyle \sum _{i=1}^{n}b{r}_{z_A}^{\mathrm{I}}{\chi }_A^{\mathrm{I}}},$$

(9)

where b is the horizontal movement factor.

Consequently, Equation (9) is the horizontal movement function of point A in the φ^I direction caused by mining of mining area I. The specific expression is shown in Equation (10):

$$u_A^{\mathrm{I}}={\displaystyle \sum _{i=1}^n\frac{b{w}_0(1-{\mathrm{e}}^{-c{t}_i})}{2}\left(f_1^{\mathrm{I}}{g}_2^{\mathrm{I}}\cos {\phi }^{\mathrm{I}}+f_2^{\mathrm{I}}{g}_1^{\mathrm{I}}\sin {\phi }^{\mathrm{I}}\right)},$$

(10)

where

$$\begin{array}{l}{f}_1^{\mathrm{I}}{=\mathrm{e}}^{-k_z^{\mathrm{I}}{}^2{x}_A^{\mathrm{I}}{}^2}-{\mathrm{e}}^{-k_z^{\mathrm{I}}{}^2{\left(x_A^{\mathrm{I}}-l\right)}^2},\\ f_2^{\mathrm{I}}{=\mathrm{e}}^{-k_z^{\mathrm{I}}{}^2{\left[y_A^{\mathrm{I}}-\left(n-i+1\right)s\right]}^2}-{\mathrm{e}}^{-k_z^{\mathrm{I}}{}^2{\left[y_A^{\mathrm{I}}-\left(n-i\right)s\right]}^2},\\ g_1^{\mathrm{I}}=\mathrm{erf}\left(k_z^{\mathrm{I}}{x}_A^{\mathrm{I}}\right)-\mathrm{erf}\left[k_z^{\mathrm{I}}\left(x_A^{\mathrm{I}}-l\right)\right],\\ g_2^{\mathrm{I}}=\mathrm{erf}\left[k_z^{\mathrm{I}}\left[y_A^{\mathrm{I}}-\left(n-i\right)s\right]\right]-\mathrm{erf}\left[k_z^{\mathrm{I}}\left[y_A^{\mathrm{I}}-\left(n-i+1\right)s\right]\right].\end{array}$$

When only mining area II is mined, assume that the coordinates of point A relative to the surface coordinate system in mining area II are ($x_A^{\mathrm{II}}$,$y_A^{\mathrm{II}}$,$z_A^{\mathrm{II}}$). According to the face occurrence conditions in Figure 2, the following equations exist:

$$x_A^{\mathrm{I}}=x_A^{\mathrm{II}},\hspace{0.33em}{y}_A^{\mathrm{I}}=2ns-y_A^{\mathrm{II}},\hspace{0.33em}{z}_A^{\mathrm{I}}=z_A^{\mathrm{II}},\hspace{0.33em}{H}_A^{\mathrm{I}}=H_A^{\mathrm{II}}=H,\hspace{0.33em}{r}_{z_A}^{\mathrm{I}}=r_{z_A}^{\mathrm{II}},\hspace{0.33em}{k}_{z_A}^{\mathrm{I}}=k_{z_A}^{\mathrm{II}},\hspace{0.33em}{\phi }^{\mathrm{I}}={\phi }^{\mathrm{II}}-\mathsf{\pi }.$$

(11)

Since mining areas I and II are distributed symmetrically, the vertical and horizontal movement functions of point A caused by mining of mining area II can be obtained by substituting Equation (11) into Equations (5) and (10).

Based on the calculations above, when mining areas I and II are mined simultaneously, the displacement of point A caused by mining is equal to

$$w_A=w_A^{\mathrm{I}}+w_A^{\mathrm{II}},\hspace{0.33em}{u}_A=u_A^{\mathrm{I}}+u_A^{\mathrm{II}}.$$

(12)

Combined with the displacement solution process of point A under the parallel mining model, the displacement under the perpendicular mining model can be obtained by the same method. The result shows that the displacement expression of point A under the perpendicular model is the same as that of the parallel model. However, the difference is the relative coordinates of point A in both model coordinate systems and φ values in the same direction.

4. Solution of the Prediction Model for Shaft Deflection

Considering that the shaft has a large slenderness ratio and small lateral bending resistance, the displacement of the point on the shaft axis can be considered to be approximately equal to the soil at the same position. Therefore, assume that point A is the point on the shaft axis, and the displacement of point A can be regarded as that of the shaft at the same depth based on the assumption of a planar cross section. Combined with the location and the deflection direction of the shaft in Figure 2, when z is used to replace z_A, the coordinates of point A can be expressed as follows: under the parallel mining model,

$$x_A^{\mathrm{I}}=x_A^{\mathrm{II}}=l+B,\hspace{0.33em}{y}_A^{\mathrm{I}}=y_A^{\mathrm{II}}=ns,\hspace{0.33em}{z}_A^{\mathrm{I}}=z_A^{\mathrm{II}}=z,\hspace{0.33em}{\phi }^{\mathrm{I}}={\phi }^{\mathrm{II}}=\mathsf{\pi };$$

(13)

under the perpendicular mining model,

$$x_A^{\mathrm{I}}=x_A^{\mathrm{II}}=l,\hspace{0.33em}{y}_A^{\mathrm{I}}=y_A^{\mathrm{II}}=ns+B,\hspace{0.33em}{z}_A^{\mathrm{I}}=z_A^{\mathrm{II}}=z,\hspace{0.33em}{\phi }^{\mathrm{I}}={\phi }^{\mathrm{II}}=\mathsf{\pi }/2.$$

(14)

The subscripts P and V represent the displacement parameters of the shaft under the parallel mining model and the perpendicular mining model, respectively. The displacement expression of the shaft for the two models can be obtained by substituting Equations (13) and (14) into Equation (12), respectively, which can be derived as

$$w_P={\displaystyle \sum _{i=1}^n\frac{w_0\left(1-{\mathrm{e}}^{-c{t}_i}\right)}{2}\left\{\mathrm{erf}\left[k_z\left(l+B\right)\right]-\mathrm{erf}\left(k_{z}B\right)\right\}}\left\{\mathrm{erf}\left(i{k}_{z}s\right)-\mathrm{erf}\left[k_z\left(i-1\right)s\right]\right\},$$

(15)

$$w_V=\mathrm{erf}\left(k_{z}l\right){\displaystyle \sum _{i=1}^n\frac{w_0\left(1-{\mathrm{e}}^{-c{t}_i}\right)}{2}\left\{\mathrm{erf}\left[k_z\left(is+B\right)\right]-\mathrm{erf}\left[k_z\left[\left(i-1\right)s+B\right]\right]\right\}},$$

(16)

$$u_P=-b{w}_0\left({\mathrm{e}}^{-k_z^2{\left(l+B\right)}^2}-{\mathrm{e}}^{-k_z^2{B}^2}\right){\displaystyle \sum _{i=1}^n(1-{\mathrm{e}}^{-c{t}_i})\left\{\mathrm{erf}\left(i{k}_{z}s\right)-\mathrm{erf}\left[\left(i-1\right)k_{z}s\right]\right\}},$$

(17)

$$u_V=b{w}_0\mathrm{erf}\left(k_{z}l\right){\displaystyle \sum _{i=1}^n(1-{\mathrm{e}}^{-c{t}_i})\left\{{\mathrm{e}}^{-k_z^2{\left[\left(i-1\right)s+B\right]}^2}-{\mathrm{e}}^{-k_z^2{\left(is+B\right)}^2}\right\}}.$$

(18)

Equations (15)–(18) are the dynamic prediction models for the shaft movement in the vertical and horizontal directions caused by multiface mining under the conditions of parallel mining and perpendicular mining. The above prediction parameters should be obtained based on the inversion of the measured data. For those parameters that cannot be inverted, the empirical values can be obtained by referring to [2].

5. Size Design of SPRP in Thick Soil and Thin Rock Strata

5.1. Parameters Solution

The Guotun Coal Mine, which has a production capacity of 2.4 Mt/y, is located in the middle of the Juye Coalfield in Shandong Province. The coal is mined by the longwall caving method. The dip angle of the coal seam is generally 5°–8°, which is a nearly horizontal coal seam. The mined face layout of the first mining area as of August 2017 is shown in Figure 5. It shows that the face layout is between the two mining models shown in Figure 2. Thus, an independent coordinate system was established for each working face and coincided with the horizontal projection of the corresponding surface coordinate system; the calculation process of the shaft deflection in the Guotun Coal Mine was consistent with the above. Considering that the vertical displacement of the shaft caused by mining has little influence on lining safety and hoisting safety, it was not necessary to analyze this. When the horizontal displacement of the shaft in the Guotun Coal Mine is represented by u_G, then its expression can be calculated as

$$u_G={\displaystyle \sum _{i=1}^n\frac{b{w}_{0i}(1-{\mathrm{e}}^{-c_i{t}_i})}{2}\left(f_{1i}\cdot g_{2i}\cos {\phi }_i+f_{2i}\cdot g_{1i}\sin {\phi }_i\right)},$$

(19)

where w_0i is the maximum subsidence value of the strata when the i-th working face is fully mined, c_i is the lithological time coefficient of the i-th working face, and φ_i is the value of φ relative to the i-th surface coordinate system:

$$\begin{gathered} f_{1i}\left(x_i,z_i\right){=\mathrm{e}}^{-k_{zi}^2{x}_i^2}-{\mathrm{e}}^{-k_i^2{\left(x_i-l_i\right)}^2}, \\ {f}_{2i}\left(y_i,z_i\right){=\mathrm{e}}^{-k_{zi}^2{y}_i^2}-{\mathrm{e}}^{-k_i^2{\left(y_i-s_i\right)}^2}, \\ {g}_{1i}\left(x_i,z_i\right)=\mathrm{erf}\left(k_{zi}{x}_i\right)-\mathrm{erf}\left[k_{zi}\left(x_i-l_i\right)\right], \\ {g}_{2i}\left(y_i,z_i\right)=\mathrm{erf}\left(k_{zi}{y}_i\right)-\mathrm{erf}\left[k_{zi}\left(y_i-s_i\right)\right], \\ {k}_{zi}=\frac{\sqrt{\mathsf{\pi }}\tan \beta }{H_i}\cdot {\left(\frac{H_i}{H_i-z_i}\right)}^{\lambda }, \end{gathered}$$

where x_i, y_i, and z_i are the coordinates of the shaft relative to the i-th surface coordinate system, respectively, and l_i, s_i, and H_i are the mining length, mining width, and mining depth of the i-th working face, respectively. The specific parameter values are shown in

Table 3. Calculation parameters of the working faces and main shaft coordinates.

i	Face Number	si (m)	li (m)	mi (m)	Hi (m)	Coordinates		North	Start Date	End Date
						xi (m)	yi (m)	φi (°)	(Year. Month)	(Year. Month)
1	1302	200	490	2.8	780	810	950	178	January 2010	July 2010
2	1304	150	870	3.2	770	785	1378	202	August 2010	March 2011
3	1301	227	1768	3.17	840	2525	692	335	November 2010	January 2013
4	1308	157	888	3.08	740	1280	1385	168	July 2011	March 2012
5	1303	230	1675	2.99	845	2535	937	335	June 2012	June 2015
			1695							August 2015
6	1310	190	790	3.2	730	796	2065	198	December 2012	August 2013
7	1305	245	985	3.41	820	2531	1194	335	September 2013	June 2015
			1250							August 2016
8	1312-1	130	640	3	730	655	2313	208	March 2014	June 2014
9	1312-2	130	1340	3	735	934	2137	188	July 2014	May 2015
10	1306	110	460	2.8	745	942	1913	208	July 2015	December 2015
11	1307	240	2100	3.3	770	2836	1445	335	April 2016	August 2017
12	1315	110	360	3.1	740	1090	1605	169	July 2016	September 2016
13	1311	136	266	3.1	785	2060	349	105	April 2017	July 2017

Note: 1312-1 and 1312-2 represent the two parts divided by the 1312 working face at the turning point.

.

In this paper, the relevant prediction parameters in deep soil strata were solved based on the main shaft deflection of the Guotun Coal Mine. The prediction parameters were as follows: According to a research report [29], q = 1.0; in light of the criterion in [2], s_0i = 0; the other parameters were obtained by inverting the measured data. First, according to the westward deflection values of the main shaft in June 2015 (see Figure 6), the inversion curve of the horizontal displacement in that direction was obtained (see curve I in Figure 6), which was basically consistent with the measured values and was suitable for Equation (19). The same parameters were then used to calculate the southward horizontal displacement (see curve II in Figure 6). The average difference between the calculated and measured values in the soil section was only 7 mm, which was less than the measurement error (8 mm). Accordingly, the prediction model could clearly invert the horizontal displacement at different shaft depths, and the reverse result showed high precision. Based on the known inversion prediction parameters, the horizontal displacement of the main shaft at different depths in 2015 (see curve I in Figure 7) and 2017 (see curve II in Figure 7) were calculated to further verify the reliability of the prediction effect. The results showed that the two theoretical values of the main shaft were largely consistent with the measured values at different depths. The prediction results were highly reliable. Based on the inversion results calculated above, the prediction parameters were λ = 0.7, b = 0.67, and tanβ = 0.51.

5.2. Influence of the Mining Model on Shaft Deflection

5.2.1. Analysis Process

Since the ideal mining models are mined more fully than the Guotun Coal Mine, the SPRP size calculated by the former models is more reliable. Therefore, referring to the parameters in Table 2, the model parameters that are suitable for thick soil and thin rock strata were selected as follows: H = 800 m, l = 1600 m, s = 200 m, v = 1200 m/y, m = 3 m, λ = 0.7, q = 1, s₀ = 0, b = 0.67, and tanβ = 0.51. When B took different values, substituting the parameters into Equations (17) and (18), respectively, the relationship between the horizontal displacement of the shaft u and the number of working face n under the two mining models could be obtained. When n took different values, the horizontal displacement curves could be calculated, as shown in Figure 8, Figure 9, Figure 10 and Figure 11. The results showed that the mining had little effect on the bedrock section of the shaft and the deflection displacement at the shaft head was the largest, expressed by u_m. To more clearly show the influence of the different mining models on the shaft deflection, the variation curves of u_m with the mining time were obtained, as shown in Figure 12, and the difference of u_m between the two models with n was calculated, as shown in Figure 13.

5.2.2. Calculation Results

As seen in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, when $n<l/s$, the horizontal displacement of the shaft under the perpendicular mining model was greater than that of the parallel mining model, and the displacement difference of the shaft head increased first and then decreased with the mining time; when n = l/s = 8, the horizontal displacement of the shaft under the two mining models was basically the same; when $n>l/s$, the horizontal displacement of the shaft was basically unchanged, indicating that the mining area had reached a state of full mining. Considering that the premature deflection is against the shaft safety, compared with the parallel mining model, the perpendicular mining model will cause greater horizontal displacement of the shaft at the initial stage of mining. Thus, the parallel mining model is superior to the perpendicular mining model.

5.3. The SPRP Designed by the Prediction Model

For a shaft constructed by the freezing method, deflection acceptance for each stage height needs to be carried out during the process of excavation and construction [30], so the deflection rate of the constructed shaft is approximately equal to 0. Therefore, when the opening reached its critical size (n = 8), the values of u_m under different values of B could be calculated based on the horizontal mining model, as shown in Figure 14.

As shown in Figure 14, with sufficient mining, the horizontal displacement value of the main shaft head in the Guotun Coal Mine was approximately 774 mm; u_m was the nonlinear negative correlation with B, and the fit expression was derived as

$$u_m=1970{\mathrm{e}}^{-{\left(\frac{B-3.707}{883.9}\right)}^2}.$$

(20)

Therefore, for the Guotun Coal Mine or others with similar mining geological conditions, the SPRPs under different shaft head displacements can be obtained by the inverse calculation of Equation (20).

5.4. The SPRP Designed by Movement Angle

5.4.1. Solution of Theoretical Soil Movement Angle

When the shaft is located at the industrial square, its SPRP size is equal to the square protection rock pillar size plus the distance from the shaft to the industrial square boundary. It is known that the main shaft is 140 m away from the industrial square boundary, the width of the enclosure belt is 20 m, and the mining depth under the mining model is 800 m. The soil strata in the first mining area are divided into the Quaternary System and the Neogene System, which have maximum thicknesses of 149.75 and 504.45 m, respectively. According to the actual movement angles in the Guotun Coal Mine, when the SPRP size was equal to 858 m, the average thickness of the soil strata under the mining model was calculated to be 640 m, which corresponded to the actual situation. Since the mining had little influence on the bedrock section of the shaft, the rock movement angle was still equal to 70°. Therefore, assume the angle that is suitable for the design of the SPRP in deep soil strata and, similar to the soil movement angle, is defined as the theoretical soil movement angle, which is expressed by α_s. According to the design requirements in the criterion [2], the main SPRP size was calculated by the parallel mining model as follows:

$$B=\frac{640}{\tan {\alpha }_s}+\frac{160}{\tan 70°}+160.$$

(21)

Substituting Equation (21) into Equation (20), the expression of α_s can be written as

$${\alpha }_s=\arctan \left(2.98/\left[4.12\sqrt{\ln \left(\frac{1970}{u_m}\right)}-1\right]\right).$$

(22)

Based on Equation (22), the relationship between the horizontal displacement of the main shaft head u_m and the theoretical soil movement angle α_s was obtained as shown in Figure 15.

5.4.2. Analysis Results

As shown in Figure 15, when the value of the theoretical movement angle for the Guotun Coal Mine was 77.8% and 57.8% of the actual movement angle (i.e., 45°), the horizontal displacement of the main shaft head was reduced by about 50% and 87%, respectively. This indicates that the shaft deflection can be effectively attenuated by appropriately reducing the soil movement angle used to design the SPRP in thick soil and thin rock strata.

6. Conclusions

Based on the probability integral method and the Knothe time function, a dynamic prediction model that can describe the horizontal movement of a shaft was established. For the main shaft deflection in the Guotun Coal Mine, a back analysis was adopted to calculate the prediction parameters that were suitable for the deep soil strata. By studying the relationship between the deflection displacement of the shaft and the SPRP, the following conclusions were obtained:

(1)

With full mining, the final displacements of the shaft under the two ideal mining models are equal, while the parallel mining model is superior at the initial stage of mining. The horizontal displacement of the shaft head has a nonlinear negative correlation with the SPRP, and the pillar size in thick soil and thin rock strata calculated by the parallel mining model is more reasonable.

(2)

Combined with the specification design method of the SPRP, the shaft deflection can be effectively attenuated by appropriately reducing the soil movement angle. For the Guotun Coal Mine, when the soil movement angle was 57.8% of the actual value, the horizontal displacement of the main shaft head was reduced by 87%.

(3)

According to the current production situation of the Guotun Coal Mine, filling mining is recommended for the coal seam around the shaft; for similar newly built shafts, symmetrical or filling mining should be adopted as much as possible according to the actual situation, and lining structures that can adapt to certain bending deformations should be selected. In addition, the inside of the freezing shaft should be constructed to be large on the top and small on the bottom, which is similar to that of the drilling shaft.