Mine Size Effects on Coal Pillar Stress and Their Application for Partial Extraction

Coal is a nonrenewable resource. Hence, it is important to improve the coal recovery ratio and ensure the stability of coal mines for sustainable development of mining cities. Partial extraction techniques, such as strip pillar mining or room-and-pillar mining, are efficient methods to extract coal. Pillar stress is a critical property for pillar design and for the assessment of mine stability after partial extraction. Current pillar stress calculation methods can sometimes overestimate the pillar stress and unnecessarily large coal pillars may be left underground, which leads to a waste of coal resources. In this paper, the size effects of mining activity on the maximum vertical pillar stress were investigated using numerical simulations. Both strip pillar mining and room-and-pillar mining were considered as possible mining scenarios at different mining depths. The results show that the maximum pillar stress of a mine is primarily controlled by four factors: the mine size to mining depth ratio, the mining width to pillar width ratio, the overburden elastic modulus, and the mining depth. The maximum pillar stress of a mine gradually increases to an ultimate value as the mine size increases. Simplified formulas and methodology have been derived for stress calculations under consideration of mine size effects and, therefore, can reduce the waste of coal resources from the overestimation of pillar stress.


Introduction
Underground mining activities can result in severe ground subsidence and damage man-made structures (buildings, railways, or other infrastructure) on the surface [1][2][3][4][5]. Mining subsidence is an important factor that restricts the sustainable development of mining cities [5]. To reduce mining subsidence and protect surface and ground structures, partial extraction methods such as strip pillar mining or room-and-pillar mining are widely adopted [5][6][7]. Figure 1 shows typical layouts of room-and-pillar mining and strip pillar mining techniques. In a partial extraction operation, massive coal pillars will be left underground to support the overburden, where the effectiveness of subsidence control for the partial extraction will depend on the stability of these pillars. Additionally, the failure of a coal pillar may result in violent ground movements [8,9]. Therefore, the stability of the coal pillars in a partial extraction is very important during and after mining activities.
The pillar stress is a critical factor affecting pillar stability, and the pillar stress calculation is an important aspect of mine design. The stress calculation for room-and-pillar mining is a three-dimensional problem because both the pillar length and the pillar width are relatively small (Figure 1a), whereas the stress calculation for strip pillar mining can be simplified to a two-dimensional problem because the pillar length is usually far longer than the pillar width ( Figure 1b). In coal mine design and pillar stability evaluation, the pillar stress is usually determined by using TAT (tributary area theory), PAT (pressure arch theory), or numerical simulations [5,[10][11][12][13][14][15][16]. TAT assumes that the overburden load will be evenly distributed among coal pillars, and the overburden load above the mined room is equally shared by neighboring pillars [5]. On the other hand, PAT assumes that a pressure arch will form between large barrier pillars and carry most of overburden weight, the abutment angle is utilized to describe the shape of the pressure arch [5,15]. For narrow mine panels, TAT overestimates the pillar stresses since pressure arches may carry the weight of the overburden. In such cases, PAT is usually adopted to calculate pillar stress [5,[14][15][16]. However, the abutment angle used by PAT is mostly summarized from longwall mining [5,[10][11][12][13][14][15][16][17][18][19][20], which is not always accurate for partial extraction. On the other hand, the shape and height of the pressure arch may change during excavation [21], indicating that the portions of the overburden load that transfer to the pillars are affected by the mine size. Therefore, the existing methods fail to provide a reliable stress estimation for partial extraction. Inaccurate pillar stress calculations usually result in unnecessarily large coal pillars that are permanently left underground; the coal resource is wasted.
To address this issue, empirically established techniques and methodologies which estimate pillar stress have been suggested [22]. He et al. [23] observed that the mine subsidence and strata movement are affected by the mine size D. Additionally, they found that the mine size D to mining depth H ratio ( Figure 1c, D/H) is an important parameter to determine whether the subsidence can reach its peak value. This observation suggests that the deformation and stress of strata may be a function of the mine size D. Roberts et al. [11] reported similar observations that the pillar stress is affected by D/H and proposed a stress calculation formula for room-and-pillar mines that included mine size effects.
In this paper, the effects of mine size variation on coal pillar stress are investigated for both strip pillar mining and room-and-pillar mining using numerical software UDEC (Universal Distinct Element Code) and 3DEC (3 Dimension Distinct Element Code). Simulation of mines includes different pillar sizes and different mine depths. The mine size refers to the mined area of the coal seam with coal pillars and mined rooms (Figure 1c). Simple formulas to calculate the maximum pillar stress for different mine sizes are derived and proposed for mine design. Based on the proposed formulas, the recovery ratio may be improved if the mine size is not too large, and the waste of coal resource can be reduced. design and pillar stability evaluation, the pillar stress is usually determined by using TAT (tributary area theory), PAT (pressure arch theory), or numerical simulations [5,[10][11][12][13][14][15][16]. TAT assumes that the overburden load will be evenly distributed among coal pillars, and the overburden load above the mined room is equally shared by neighboring pillars [5]. On the other hand, PAT assumes that a pressure arch will form between large barrier pillars and carry most of overburden weight, the abutment angle is utilized to describe the shape of the pressure arch [5,15]. For narrow mine panels, TAT overestimates the pillar stresses since pressure arches may carry the weight of the overburden. In such cases, PAT is usually adopted to calculate pillar stress [5,[14][15][16]. However, the abutment angle used by PAT is mostly summarized from longwall mining [5,[10][11][12][13][14][15][16][17][18][19][20], which is not always accurate for partial extraction. On the other hand, the shape and height of the pressure arch may change during excavation [21], indicating that the portions of the overburden load that transfer to the pillars are affected by the mine size. Therefore, the existing methods fail to provide a reliable stress estimation for partial extraction. Inaccurate pillar stress calculations usually result in unnecessarily large coal pillars that are permanently left underground; the coal resource is wasted.
To address this issue, empirically established techniques and methodologies which estimate pillar stress have been suggested [22]. He et al. [23] observed that the mine subsidence and strata movement are affected by the mine size . Additionally, they found that the mine size to mining depth ratio (Figure 1c, ⁄ ) is an important parameter to determine whether the subsidence can reach its peak value. This observation suggests that the deformation and stress of strata may be a function of the mine size . Roberts et al. [11] reported similar observations that the pillar stress is affected by ⁄ and proposed a stress calculation formula for room-and-pillar mines that included mine size effects.
In this paper, the effects of mine size variation on coal pillar stress are investigated for both strip pillar mining and room-and-pillar mining using numerical software UDEC (Universal Distinct Element Code) and 3DEC (3 Dimension Distinct Element Code). Simulation of mines includes different pillar sizes and different mine depths. The mine size refers to the mined area of the coal seam with coal pillars and mined rooms (Figure 1c). Simple formulas to calculate the maximum pillar stress for different mine sizes are derived and proposed for mine design. Based on the proposed formulas, the recovery ratio may be improved if the mine size is not too large, and the waste of coal resource can be reduced.

Mine Modeling and Simulation
The discrete element modeling program, UDEC, is used in the current study. The scheme of the model is shown in Figure 1c. The ground surface and the coal seam are modeled as horizontal strata, where the thicknesses of the coal seam and the mine floor are modeled as 7 m and 20 m, respectively. The mining depth H is varied from 100 m to 800 m; thus, the models included both shallow and deep mines. The pillar width w p is assumed to range from 7 m to 134 m. The mining width w c of strip pillar mining is usually 0.1 H-0.25 H to control the ground subsidence, thus, the pillars with very large sizes are proposed mainly for the theoretical analysis of strip pillar mining at deep depths. The ratio of mining width to pillar width r m is varied from 0.75 to 2. The geometric parameters of the models are listed in Table 1  In the model, all non-coal strata are assumed to consist of one type of rock. The properties of non-coal strata are chosen to mimic sandstone except in Section 3.3. In Section 3.3, the elastic modulus and Poisson's ratio are changed while the other properties of non-coal strata refer to sandstone to study the effect of different rock properties on pillar stress. Rock and coal properties are taken from the Yulin mine area in China and are listed in Table 2 [24]. The Mohr-Coulomb failure criterion is assumed for the models. Since the joint stiffness only affects the joint deformation, a large stiffness (5 × 10 9 Pa) is assigned to the interface between non-coal strata and the coal seam, so that the maximum stresses in pillars can be calculated. The strength parameters of the interfaces are matched to the coal strength. The model base is fixed in the vertical direction and the lateral boundary of the model is fixed in the horizontal direction. Static equilibrium analysis is conducted with consideration of gravitational loads on the system. Every pillar contained one vertical observation line located at the pillar center. The model is mined one tunnel at a time and from the left-hand side until the mine reached the design size (Table 1). For each time the tunnel number (or mine size D) increases, the maximum vertical pillar stress of the whole mine is recorded. Table 2. Mechanical properties of strata: σ T is the cohesion, σ C is the tensile strength, ϕ is the internal friction angle, E is the elastic modulus, υ is the Poisson's ratio, and ρ is the density.

Analysis of Simulation Results
Mine failures did not occur in all of the numerical simulations. The numerical simulation results indicate that the maximum vertical stress of a single pillar appeared at the pillar bottom, and the maximum vertical pillar stress of the mine is dominated by the mine size D, the mining width to pillar width ratio r m , and the overburden elastic modulus E o . In this study, the maximum vertical pillar stress of the mine is expressed by the stress concentration factor K: where g is the acceleration due to gravity, ρ 1 and ρ 2 are the density of non-coal strata and the density of the coal seam, respectively; H and h are mining depth and pillar height, respectively, and S is the maximum vertical pillar stress of the mine from the simulation.

Effects of the Mine Size on the Maximum Vertical Pillar Stress of the Mine
By fixing r m to unity, the mine size effect on K was investigated. The results ( Figure 2) show that there is an exponential relationship between K and D/H: where K S is the ultimate stress concentration factor when the pillar stress is stable and reaches its peak value ( Figure 2); e is the natural base, and a is a constant dependent on H, r m , and E o .
Sustainability 2018, 10, x FOR PEER REVIEW 4 of 11 pillar center. The model is mined one tunnel at a time and from the left-hand side until the mine reached the design size (Table 1). For each time the tunnel number (or mine size ) increases, the maximum vertical pillar stress of the whole mine is recorded.

Analysis of Simulation Results
Mine failures did not occur in all of the numerical simulations. The numerical simulation results indicate that the maximum vertical stress of a single pillar appeared at the pillar bottom, and the maximum vertical pillar stress of the mine is dominated by the mine size , the mining width to pillar width ratio , and the overburden elastic modulus . In this study, the maximum vertical pillar stress of the mine is expressed by the stress concentration factor : where is the acceleration due to gravity, and are the density of non-coal strata and the density of the coal seam, respectively; and are mining depth and pillar height, respectively, and S is the maximum vertical pillar stress of the mine from the simulation.

Effects of the Mine Size on the Maximum Vertical Pillar Stress of the Mine
By fixing to unity, the mine size effect on was investigated. The results ( Figure 2) show that there is an exponential relationship between and ⁄ : where is the ultimate stress concentration factor when the pillar stress is stable and reaches its peak value ( Figure 2); e is the natural base, and is a constant dependent on , , and .  increases when ⁄ increases, and gradually, become stable and reaches its peak value when ⁄ ≥3 to 4. Therefore, the parameter ⁄ can be an indication to determine if the pillar stress in the mine reaches its maximum value; if ⁄ <3 to 4, the pillar stress of the mine is not fully developed and will increase with further increase of ⁄ . Otherwise, reaches its maximum value and will no longer increase with the increase of ⁄ .  determine if the pillar stress in the mine reaches its maximum value; if D/H <3 to 4, the pillar stress of the mine is not fully developed and K will increase with further increase of D/H. Otherwise, K reaches its maximum value and will no longer increase with the increase of D/H.
The parameters in Equation (2) are listed in Table 3, which indicates that a is dependent on mining depth but is independent of the mining and pillar widths. For each mining depth, the mean value of a can be used to estimate the maximum pillar stress K (Table 4). It is shown in Figure 3 that the average coefficient a increases as the mining depth increases. The parameters in Equation (2) are listed in Table 3, which indicates that is dependent on mining depth but is independent of the mining and pillar widths. For each mining depth, the mean value of can be used to estimate the maximum pillar stress (Table 4). It is shown in Figure 3 that the average coefficient increases as the mining depth increases.  It is found that the maximum pillar stress of the mine can be calculated by the TAT method when ⁄ ≥3 to 4. Therefore, for a two-dimensional plane (Figure 1c), of coal pillars for a mine with ⁄ ≥3 to 4 can be calculated as: and Equation (2) can be expressed as: It is found that the maximum pillar stress of the mine can be calculated by the TAT method when D/H ≥3 to 4. Therefore, for a two-dimensional plane (Figure 1c), K S of coal pillars for a mine with D/H ≥3 to 4 can be calculated as: and Equation (2) can be expressed as:

Effects of r m on the Maximum Vertical Pillar Stress of the Mine
In order to control ground subsidence, the mining width for partial extraction is usually set between 0.1 H and 0.25 H [6]. To study the effects of r m , the mining width is assumed to be in this range, the mining depth H is varied from 100 m to 400 m, and the r m is varied from 0.75 to 2. The results show that the K curves for different r m values are similar to the K curves shown in Figure 2. The parameters for Equation (2) under different r m are listed in Table 5. Tables 3 and 5 indicate that the performance of the fit parameter a is primarily dominated by r m , followed by the effect of mining depth. Furthermore, from Table 5, it is noticed that K S is approximately equal to 1 + r m , indicating that Equation (4) is valid. As shown in Figure 3, the parameter a increases by 0.06 when the mining depth increases by 100 m. In Table 5, the mining depth varies from 100 m to 400 m, thus, the maximum variation of a induced by the change of mining depth is only 0.18. However, a approximately varies from 0.4 to 2 in Table 5. This means that a is mainly dominated by r m and the mining depth effect is negligible. In turn, the average values of a and K S for each r m can be used to represent the r m effect, as shown in Table 6.  Figure 4 shows the relationship between r m and the average parameters K S and a. An increase in r m results in an increase in K S but a decrease in a. As the variation of a is mainly controlled by r m , and the mining depth effect on a can be neglected, Figure 4a shows that a will increase by −1.2723 if r m increases by a value of 1. Theoretically, K S should be equal to 1 + r m . However, Figure 4b shows that when r m is small, K S is slightly smaller than 1 + r m ; and when r m is large, K S is slightly larger than 1 + r m . This is because the pillar stress is not uniformly distributed [5,10]. The outer parts of a pillar support more overburden load than the inner part of the pillar [10]. Thus, when r m is small, the vertical stress at the pillar center is smaller than the average stress on the pillar. When r m is large, the outer parts of the pillar may yield and lose capacity due to the large overburden load, and the effective pillar size decreases, leading to a situation where the stress at the pillar center is larger than the average pillar stress. As a result, K S is not totally equal to the theoretical value of 1 + r m . is large, the outer parts of the pillar may yield and lose capacity due to the large overburden load, and the effective pillar size decreases, leading to a situation where the stress at the pillar center is larger than the average pillar stress. As a result, is not totally equal to the theoretical value of 1 + .

Effects of Overburden Properties on the Maximum Vertical Pillar Stress of the Mine
To study the effects of overburden properties, the model with a mining depth of 200 m, a mining width of 20 m, and an of one (from Table 1) is selected. First, the overburden elastic modulus, , is set at 15 GPa, and the Poisson's ratio, , is varied from 0.1 to 0.49. Then is fixed at 0.15, and is varied from 2 GPa to 100 GPa. Simulation results show that for different overburden properties are consistent with Equation (2). Table 7 lists the parameters of Equation (2) for different overburden properties.  Figure 5 shows the effect of on and . When is small, and are strongly affected by and the overburden elastic modulus cannot be neglected.

Effects of Overburden Properties on the Maximum Vertical Pillar Stress of the Mine
To study the effects of overburden properties, the model with a mining depth of 200 m, a mining width of 20 m, and an r m of one (from Table 1) is selected. First, the overburden elastic modulus, E o , is set at 15 GPa, and the Poisson's ratio, υ, is varied from 0.1 to 0.49. Then υ is fixed at 0.15, and E o is varied from 2 GPa to 100 GPa. Simulation results show that K for different overburden properties are consistent with Equation (2). Table 7 lists the parameters of Equation (2) for different overburden properties. Table 7. Parameters of Equation (2) for different overburden properties. K S /(1 + r m ) is the ratio of simulated ultimate stress concentration factor K S to the theoretical K S , representing a reduction effect of the overburden elastic modulus E o on simulated K S .  Figure 5 shows the effect of E o on a and K S . When E o is small, a and K S are strongly affected by E o and the overburden elastic modulus cannot be neglected.  Figure 5 shows the effect of on and . When is small, and are strongly affected by and the overburden elastic modulus cannot be neglected.  According to Figure 5a, a can be represented by E o as:

Poisson's Ratio Effect
Combining Equations (4) and (5), K can be calculated as: Equation (6) is based on the model with r m = 1 and mining depth H of 200 m. As shown in Figures 3 and 4b, a linearly increases by 0.0006 if H increases by 1 m, and a linearly increases by −1.273 if r m increases by 1. Thus, K for different r m at different mining depth can be estimated by introducing the mining depth and r m effects into Equation (6): Table 8 shows the relative errors of K as calculated by Equation (7) when compared with the simulation results. It is found that when r m is large (i.e., r m = 2), Equation (7) may slightly underestimate the pillar stress. Otherwise, Equation (7) usually overestimates the pillar stress. As stated, this is due to the fact that the stress on the pillars is not uniformly distributed. The stress at the pillar center was acquired in the simulations, and there exist some differences between the stress calculated by Equation (7) and the simulated stress, as previously discussed. Although there exist some differences between calculated stress and simulated stress, the relative errors are small for most of the cases. Hence, Equation (7) can be used to quickly estimate the pillar stress for different mine sizes.

Mine Size Effect on 3D Partial Extraction
As UDEC is only a two-dimensional simulation software, Equation (7) is only suitable for strip pillar mining. For room-and-pillar mines, the coal pillars bear three-dimensional stresses. Therefore, 3DEC is utilized for three-dimensional stress studies. Figure 6 shows a horizontal coal seam in a three-dimensional model, and the coal seam is mined in the north-south and east-west directions.  According to Equation (7), mining in the north direction will contribute to a maximum pillar stress as: and mining in the east direction will contribute to a maximum pillar stress as: The maximum pillar stress of the three-dimensional model is the superposition of and from the two directions, and can be estimated as: (10) Figure 7 shows the simulated by 3DEC and the calculated by Equation (10). The results show that the curves are similar to curves of the two-dimensional simulations; the difference between simulated and calculated is very small. The curve for Figure 6b ( ⁄ >3 to 4 but ⁄ <3 to 4) is larger than the curve for Figure 6a ( ⁄ >3 to 4 and ⁄ >3 to 4). The value will gradually approach the value calculated by the TAT method with bi-directional ⁄ >3.  The coal seam is mined in the north and east directions simultaneously, beginning from the south-west corner of the model ( Figure 6). Assume the mine size in the north direction is D N , the mine size in the east direction is D E . Two scenarios have been investigated: (1) D N /H >3 to 4 but D E /H <3 to 4 (Figure 6a). In this case, D N and D E increase simultaneously, and the increase in D E stops if five pillars form in the east direction, while the increase in D N stops if 20 pillars form in the north direction. It was found that D N /H =4.1 and D E /H =1.1 when the mining activity stopped. (2) D N /H >3 to 4 and D E /H >3 to 4 ( Figure 6b). In this case, D N and D E increase simultaneously until 20 pillars are formed bi-directionally, and D N /H and D E /H >3 were found when mining activity stopped.
According to Equation (7), mining in the north direction will contribute to a maximum pillar stress K N as: and mining in the east direction will contribute to a maximum pillar stress K E as: The maximum pillar stress of the three-dimensional model K 3D is the superposition of K N and K E from the two directions, and can be estimated as: (10) Figure 7 shows the K 3D simulated by 3DEC and the K 3D calculated by Equation (10). The results show that the K 3D curves are similar to K curves of the two-dimensional simulations; the difference between simulated K 3D and calculated K 3D is very small. The K 3D curve for Figure 6b  K 3D value will gradually approach the K 3D value calculated by the TAT method with bi-directional D/H >3. from the two directions, and can be estimated as: (10) Figure 7 shows the simulated by 3DEC and the calculated by Equation (10). The results show that the curves are similar to curves of the two-dimensional simulations; the difference between simulated and calculated is very small. The curve for Figure 6b ( ⁄ >3 to 4 but ⁄ <3 to 4) is larger than the curve for Figure 6a ( ⁄ >3 to 4 and ⁄ >3 to 4). The value will gradually approach the value calculated by the TAT method with bi-directional ⁄ >3.

Discussion
The TAT method can only be used when the mine size is large enough, but the critical mine size for which the TAT method is suitable has not been well addressed prior to the current study. All the numerical simulations in this paper show that the maximum pillar stress of a mine will gradually reach a peak value (the pillar stress calculated by the TAT method) when D/H ≥3 to 4. Thus, D/H of 3 is suggested as an index to determine whether the TAT method can be used. This is important especially for mining at deep depths, because the D/H may be small, and using the TAT method may overestimate the pillar stress, resulting in unnecessarily large pillars and wasted coal resource.
As no mine failures occurred in these simulations, the K values represent the maximum pillar stress in a mine. For strip mining, Equation (7) can be used to calculate the pillar stress when D/H <3 to 4; otherwise, both Equation (7) and TAT can be used. For room-and-pillar mines, the TAT method can be used if D N /H >3 to 4 and D E /H >3 to 4; otherwise, Equations (8)-(10) should be used. By doing this, the waste of coal resource that results from an overestimation of pillar stress can be reduced.
Equations (8)- (10) are proposed based on the analysis of two-dimensional UDEC models. The small difference between 3DEC simulated stress and calculated stress by formula (Figure 7) indicated that the parameters from two-dimensional simulation can be used to predict the pillar stress in a three-dimensional model. Therefore, two-dimensional UDEC simulation can also be conducted to determine the parameters a and K S in Equations (8)-(10) even for complicated mining situations. This is important because three-dimensional numerical simulations are more time-consuming than two-dimensional numerical simulations. A numerical model with D/H ≥3 to 4 is also suggested when analyzing the mine system in full scale.

Conclusions
The widely used TAT method may overestimate the pillar stress since it neglects mine size effects. Essentially, mine size effects result from pressure arch formation, but the arch shape changes during mining and is hard to determine; current PAT methods are also not always reliable in pillar stress calculations. As a result, more precise stress calculation methods are needed.
Numerical studies of partial extraction mining are presented and several simplified pillar stress formulas are presented for both two-dimensional and three-dimensional pillar stress calculations. It is observed that the mine size has the same effects on both strip mining and room-and-pillar mining. For both shallow and deep mines, the maximum vertical pillar stress of a mine is controlled by D/H, r m , E o , and H. A D/H of 3 to 4 has been found to be useful in determining whether the pillar stress can be calculated by the TAT method. As long as D/H <3 to 4 in any mining direction, the TAT method overestimates the pillar stress, and mine size effects should be considered. Thus, the pillar size can be reduced, and the recovery ratio of partial extraction can be improved. Finally, the parameters a and K S for the three-dimensional stress calculation can be achieved through two-dimensional numerical simulation.