# Numerical Analysis of the Width Design of a Protective Pillar in High-Stress Roadway: A Case Study

Authors: FuZhou Qi; ZhanGuo Ma; Ning Li; Bin Li; Zhiliu Wang; WeiXia Ma
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The width design of protective pillars is an important factor affecting the stability of high-stress roadways. In this study, a novel numerical modeling approach was developed to investigate the relationship between protective pillar width and roadway stability. With the 20 m protective pillar width adopted in the field test, large deformation of roadways and serious damage to surrounding rocks occurred. According to the case study at the Wangzhuang coal mine in China, the stress changes and energy density distribution characteristics in protective pillars with various widths were analysed by numerical simulation. The modeling results indicate that, with a 20 m wide protective pillar, the peak vertical stress and energy density in the pillar are 18.5 MPa and 563.7 kJ/m3, respectively. The phenomena of stress concentration and energy accumulation were clearly observed in the simulation results, and the roadway is in a state of high stress. Under the condition of a 10 m wide protective pillar, the peak vertical stress and energy density are shifted from the pillar to roadway virgin coal region, with peak values of 9.5 MPa and 208.3 kJ/m3, respectively. The decrease in vertical stress and energy density improves the stability of the protective pillar, resulting in the roadway being in a state of low stress. Field monitoring suggested that the proposed 10 m protective pillar width can effectively control the large deformation of the surrounding rock and reduce coal bump risk. The novel numerical modeling approach and design principle of protective pillars presented in this paper can provide useful references for application in similar coal mines.


---

## Body

## 1. Introduction

With increased mining depth, the failure of protective pillars in deep high-stress roadways has received much attention [1]. The influence of the protective pillar design on the stability of roadways has a bearing on worker safety and the effective mining of coal resources. Protective pillar around the gob provides a safe and stable production environment for the roadway and decreases disturbance to the roadway from high external stress. An improper protective pillar design can cause frequent roof sagging, pillar rib bulging, and severe coal bump to occur in a roadway, which will eventually cause the collapse of the surrounding rock structure [2, 3]. Therefore, protective pillar width design is a necessity for maintaining the stability of the surrounding rock. Many methods of designing an appropriate protective pillar width have been proposed, including theoretical analysis calculations and numerical simulations.A variety of calculation models have been established to analyse the proper width of a protective pillar and decrease the risk of pillar failure [4, 5]. Ghasemi et al. [6] used limit equilibrium analysis to explore the relationship between the range of elastic area and the stability of the protective pillar. Cao and Zhou [7] investigated the influence of movement in key blocks of broken strata above the roadway on the stability of the protective pillar and defined a reasonable width for the pillar. Yang et al. [8] showed the scope of stress distribution in the upper and lower strata of a protective pillar by analysing the influence of roof failure on the bearing capacity of the pillar.In reality, many factors, such as the protective pillar’s form, the rock mechanical properties, the stress state, and the panel mining methods, should be considered in designing the width of the protective pillar. Since some of these factors are not considered in theoretical analysis calculations, they are limited in application. Various numerical simulations have been adopted to analyse the effect of the protective pillar width on the roadway stability because of their low cost and high efficiency [9]. Moreover, as numerical simulation includes numerous factors during modeling, the results are often more reliable than those of theoretical calculations. Shabanimashcool and Li [10] established local and global numerical models to study the influence of falling strata on the stress and breaking range of a protective pillar during panel retreat. Bai et al. [11] investigated roof failure characteristics and stress change rules of gob-side entry driving by stages, noting that a rational protective pillar width could enhance the roadway stability. Ma et al. [12] analysed the shear failure behavior for protective pillars and found that the failure mode varied with the pillar width and inclination of the coal seam. Wang et al. [13] adopted FLAC3D software to study the dynamic response and failure mechanism of protective pillars; the results showed that enlarging the elastic core in the pillar increased the risk of coal bump. However, based on previous studies, we found that the numerical simulations seldom considered the following: (1) When the overlying strata in a gob are cut off along a preexisting fissure, the caving materials will fill the gob and provide supporting resistance to the roof strata, which will relieve stress on the protective pillar and decrease its width. (2) The energy stored in the rock mass is closely related to its stress state. The stored energy will vary as the stress changes. Thus, the evolution laws of energy stored in a rock mass can be inferred by the change in stress, decreasing or avoiding the deformation failure of the surrounding rock and the occurrence of coal bump. In view of these limitations in current simulations and considering the supporting feature of gob caving materials and the relationship between the stress and energy of the rock mass, a numerical modeling is developed to analyse failure mechanism for protective pillars.In this paper, the relationship between roadway stability and protective pillar width is investigated based on numerical simulations and field tests. First, the failure characteristics of high-stress roadways are studied when the width of protective pillar is 20 m. Next, a novel modeling approach is proposed. In this modeling, a strain-softening model is adopted to describe the mechanical behaviour of protective pillars, and a double-yield model is used to simulate the gob materials mechanical behaviour. The stress changes and energy density distribution characteristics of roadways with five different widths are analysed. Finally, the effect of the optimal protective pillar on controlling the deformation of the surrounding rock is evaluated by a field test. The modeling approach and design principle presented in this paper can be used to analyse the protective pillar design of high-stress roadways at other similar sites.


## 2. Engineering Background

### 2.1. Geology Conditions

As shown in Figure1, the selected study site is in the Wangzhuang coal mine, Shanxi Province, China. The average depth and thickness of coal seam are 300 m and 4 m. Roadway roof and floor mainly include mudstone, siltstone, and sandy mudstone. A detailed stratigraphic column is shown in Figure 2.Figure 1
Location of the Wangzhuang coal mine.Figure 2
Detailed stratigraphic column for the test site.Panels 8101, 8102, and 8103 are located in mining area #8 of the test coal mine, with a strike length of 1088 m and dip length of 210 m, as shown in Figure3. Panel 8101 is located to the north of panel 8102, where the coal seam has been exploited. Panel 8103 is located to the south of panel 8102, which is unexploited. At present, the width of the protective pillar between the 8101 tailgate and 8102 headgate is 20 m (Figure 3).Figure 3
Layout of roadway and panel at the test site.


### 2.2. Field Monitoring

In order to evaluate the current protective pillar stability, the roadway deformation for 8102 headgate was measured in the process of roadway excavation and panel retreat. Four measuring stations (solid magenta circles in Figure3) were arranged in the 8102 headgate. The distance between measuring stations was 50 m. The field measuring process included the following work:(1)
The pins (red solid circle in Figure4(a)) were fixed to the corresponding regions of the roadway at each measuring station. Floor and roof pins were fixed in the midspan; the distance of pins in two ribs from the roof was 1.7 m.(2)
In the observation process, the convergence of the roof, floor, and two ribs was measured by using a digital deformation instrument and measuring line (Figure4(b)). The digital deformation instrument was accurate to 0.01 mm, which satisfied the needs of deformation observation.(3)
During the measuring period, the field data were measured once every 3 days.Figure 4
Roadway convergence measurement in field. (a) Layout of the measurement station. (b) Measurement devices.
(a)(b)The results of the 8102 headgate displacements at different stages are shown in Figure5. Note that the displacement of the 8102 headgate in Figure 5 is the average displacement of the four measuring stations in the field. The roadway deformation tended to be stable after excavation of 57 days, with a rate of deformation close to zero (Figure 5(a)). The total convergence of the roof, protective pillar rib, and virgin coal rib reached 286 mm, 206 mm, and 158 mm, respectively. As shown in Figure 5(b), compared with the convergence during excavation, the roadway displacement grew rapidly in the panel retreat period. The convergence of roadway roof increased to 991 mm, and the protective pillar rib and virgin coal rib convergence reached 661 mm and 396 mm, respectively. More than 80% of the roadway deformation occurred in a location where measuring stations advanced approximately 40 m from the longwall face. The results show that the deformation of the roof was the largest and that the deformation in the protective pillar rib was the next largest; the deformation on the floor was considerably smaller than that of the roof and ribs. The immediate roof and two ribs of headgate 8102 were located in the coal seam with low strength, and the floor was located in mudstone with relatively high strength. This special engineering geological condition led to a great amount of deformation in the roof and the two ribs once the roadway was excavated.Figure 5
Measured convergences in 8102 headgate. (a) During 8102 headgate excavation. (b) During panel 8102 retreat.
(a)(b)Severe roof sagging and pillar rib bulging (Figures6(a) and 6(b)) occurred frequently in the process of panel 8102 retreat. Approximately 30 m in front of panel 8102, roof fall and small-scale coal bump accidents (solid blue circle in Figure 3) occurred, as shown in Figures 6(c) and 6(d). Fortunately, no injuries occured in the accidents; however, the production efficiency was affected. It took about 20 days to repair the roadway, wasting large amounts of time and money. Professionals did not find the cause of the roof caving and coal bump. It was presumed that the coupling action of high mining-induced stress and low coal strength led to instability of the protective pillar, resulting in the accidents. Therefore, the roof stratum and protective pillar rib were the areas that received attention in panel retreat process.Figure 6
Failure characteristics of 8102 headgate. (a) Roof sagging. (b) Protective pillar bulging. (c) Roof fall. (d) Coal bump.
(a)(b)(c)(d)


## 2.1. Geology Conditions

As shown in Figure1, the selected study site is in the Wangzhuang coal mine, Shanxi Province, China. The average depth and thickness of coal seam are 300 m and 4 m. Roadway roof and floor mainly include mudstone, siltstone, and sandy mudstone. A detailed stratigraphic column is shown in Figure 2.Figure 1
Location of the Wangzhuang coal mine.Figure 2
Detailed stratigraphic column for the test site.Panels 8101, 8102, and 8103 are located in mining area #8 of the test coal mine, with a strike length of 1088 m and dip length of 210 m, as shown in Figure3. Panel 8101 is located to the north of panel 8102, where the coal seam has been exploited. Panel 8103 is located to the south of panel 8102, which is unexploited. At present, the width of the protective pillar between the 8101 tailgate and 8102 headgate is 20 m (Figure 3).Figure 3
Layout of roadway and panel at the test site.


## 2.2. Field Monitoring

In order to evaluate the current protective pillar stability, the roadway deformation for 8102 headgate was measured in the process of roadway excavation and panel retreat. Four measuring stations (solid magenta circles in Figure3) were arranged in the 8102 headgate. The distance between measuring stations was 50 m. The field measuring process included the following work:(1)
The pins (red solid circle in Figure4(a)) were fixed to the corresponding regions of the roadway at each measuring station. Floor and roof pins were fixed in the midspan; the distance of pins in two ribs from the roof was 1.7 m.(2)
In the observation process, the convergence of the roof, floor, and two ribs was measured by using a digital deformation instrument and measuring line (Figure4(b)). The digital deformation instrument was accurate to 0.01 mm, which satisfied the needs of deformation observation.(3)
During the measuring period, the field data were measured once every 3 days.Figure 4
Roadway convergence measurement in field. (a) Layout of the measurement station. (b) Measurement devices.
(a)(b)The results of the 8102 headgate displacements at different stages are shown in Figure5. Note that the displacement of the 8102 headgate in Figure 5 is the average displacement of the four measuring stations in the field. The roadway deformation tended to be stable after excavation of 57 days, with a rate of deformation close to zero (Figure 5(a)). The total convergence of the roof, protective pillar rib, and virgin coal rib reached 286 mm, 206 mm, and 158 mm, respectively. As shown in Figure 5(b), compared with the convergence during excavation, the roadway displacement grew rapidly in the panel retreat period. The convergence of roadway roof increased to 991 mm, and the protective pillar rib and virgin coal rib convergence reached 661 mm and 396 mm, respectively. More than 80% of the roadway deformation occurred in a location where measuring stations advanced approximately 40 m from the longwall face. The results show that the deformation of the roof was the largest and that the deformation in the protective pillar rib was the next largest; the deformation on the floor was considerably smaller than that of the roof and ribs. The immediate roof and two ribs of headgate 8102 were located in the coal seam with low strength, and the floor was located in mudstone with relatively high strength. This special engineering geological condition led to a great amount of deformation in the roof and the two ribs once the roadway was excavated.Figure 5
Measured convergences in 8102 headgate. (a) During 8102 headgate excavation. (b) During panel 8102 retreat.
(a)(b)Severe roof sagging and pillar rib bulging (Figures6(a) and 6(b)) occurred frequently in the process of panel 8102 retreat. Approximately 30 m in front of panel 8102, roof fall and small-scale coal bump accidents (solid blue circle in Figure 3) occurred, as shown in Figures 6(c) and 6(d). Fortunately, no injuries occured in the accidents; however, the production efficiency was affected. It took about 20 days to repair the roadway, wasting large amounts of time and money. Professionals did not find the cause of the roof caving and coal bump. It was presumed that the coupling action of high mining-induced stress and low coal strength led to instability of the protective pillar, resulting in the accidents. Therefore, the roof stratum and protective pillar rib were the areas that received attention in panel retreat process.Figure 6
Failure characteristics of 8102 headgate. (a) Roof sagging. (b) Protective pillar bulging. (c) Roof fall. (d) Coal bump.
(a)(b)(c)(d)


## 3. Protective Pillar Design Principle

### 3.1. Stress Distribution Laws of Protective Pillars

With roadway excavation, roof strata loads, initially borne by solid coal, are transferred to the two roadway ribs. The distribution of vertical stress in the protective pillar varies greatly with the variation in the pillar width. With a protective pillar width of 5–20 m, the vertical stress distribution characteristics are described by a “single peak,” while the vertical stress distribution characteristics are described by a “double peak” with a 20–40 m wide protective pillar, as shown in Figure7.Figure 7
Stress distribution of the protective pillar [35]. (a) Single-peak distribution. (b) Double-peak distribution.
(a)(b)


### 3.2. Relationship between Protective Pillar Width and Roadway Stability

The relationship between the width of protective pillar and roadway convergence was used to evaluate the stability of the roadway, as shown in Figure8.Figure 8
The protective pillar width and roadway convergence.For determining a rational pillar width for a high stress roadway, the following considerations could be made:(1)
The width of the protective pillar could not be too thin, such as the left side of point A in Figure8. Because of its weak bearing capacity, the protective pillar could not maintain roadway stability in this case.(2)
The protective pillar width could not be too wide, indicated by the curve between points B and D. On account of the strain energy accumulation, the protective pillar may lead to coal bumps in this situation.(3)
The width of the protective pillar could not reach the critical width, such as point B. Under this condition, the peak stress would be shifted from the solid coal region to the protective pillar. The bearing strength of the protective pillar could not resist the peak stress effectively, and the roadway was in a state of high stress.A suitable width for protective pillars should provide enough capacity to bear roof strata loads while maintaining the pillar integrity, reducing the pillar width, and increasing coal recovery rate. Therefore, a thin pillar could be chosen, corresponding to the left side of point B (Figure8). The roadway could undertake a certain degree of large deformation in this condition, and its integrity could be maintained.


### 3.3. Determination of the Cutting Height for the Overlying Strata in the Gob

During mining, as the working face constantly moved forward, the strata behind the working face were cut off along the preexisting fissure. The gob caving material had a certain bearing capacity when compressed by the upper strata. The gob caving material could support the overlying strata and reduce the pressure on the pillars. Therefore, when designing protective pillars, it is necessary to account for the supporting characteristics of the caving materials.To effectively utilize the supporting characteristics of the caving materials, it is important to determine the cutting height of the overlying strata in the gob. After roof stratum caving, the volume of caving materials expands. In other words, the pile height of the caving rock mass will be larger than the original height. The volume expansion of the caving materials can be fully described with the bulking factorKp. The intact rock original height is Σh, and the rock mass pile height after roof strata collapse isKpΣh. The stratum void height Δ between the piled rock and the overlying stable rock will be(1)Δ=∑h+M−Kp∑h=M−∑hKp−1,where M denotes the coal seam height andKp denotes the bulking factor.According to equation (1), when Δ = 0, the cutting height of the roof strata is calculated by the following equation [15, 16]:(2)∑h=MKp−1.Based on a number of field observations, the caving height is 4 to 8 times of the coal seam height [17, 18]. For panel 8102, the coal seam height is 4 m, and the bulking factor is determined to be 1.25 [19, 20]. Therefore, the cutting height of the strata above the gob is 16 m.


## 3.1. Stress Distribution Laws of Protective Pillars

With roadway excavation, roof strata loads, initially borne by solid coal, are transferred to the two roadway ribs. The distribution of vertical stress in the protective pillar varies greatly with the variation in the pillar width. With a protective pillar width of 5–20 m, the vertical stress distribution characteristics are described by a “single peak,” while the vertical stress distribution characteristics are described by a “double peak” with a 20–40 m wide protective pillar, as shown in Figure7.Figure 7
Stress distribution of the protective pillar [35]. (a) Single-peak distribution. (b) Double-peak distribution.
(a)(b)


## 3.2. Relationship between Protective Pillar Width and Roadway Stability

The relationship between the width of protective pillar and roadway convergence was used to evaluate the stability of the roadway, as shown in Figure8.Figure 8
The protective pillar width and roadway convergence.For determining a rational pillar width for a high stress roadway, the following considerations could be made:(1)
The width of the protective pillar could not be too thin, such as the left side of point A in Figure8. Because of its weak bearing capacity, the protective pillar could not maintain roadway stability in this case.(2)
The protective pillar width could not be too wide, indicated by the curve between points B and D. On account of the strain energy accumulation, the protective pillar may lead to coal bumps in this situation.(3)
The width of the protective pillar could not reach the critical width, such as point B. Under this condition, the peak stress would be shifted from the solid coal region to the protective pillar. The bearing strength of the protective pillar could not resist the peak stress effectively, and the roadway was in a state of high stress.A suitable width for protective pillars should provide enough capacity to bear roof strata loads while maintaining the pillar integrity, reducing the pillar width, and increasing coal recovery rate. Therefore, a thin pillar could be chosen, corresponding to the left side of point B (Figure8). The roadway could undertake a certain degree of large deformation in this condition, and its integrity could be maintained.


## 3.3. Determination of the Cutting Height for the Overlying Strata in the Gob

During mining, as the working face constantly moved forward, the strata behind the working face were cut off along the preexisting fissure. The gob caving material had a certain bearing capacity when compressed by the upper strata. The gob caving material could support the overlying strata and reduce the pressure on the pillars. Therefore, when designing protective pillars, it is necessary to account for the supporting characteristics of the caving materials.To effectively utilize the supporting characteristics of the caving materials, it is important to determine the cutting height of the overlying strata in the gob. After roof stratum caving, the volume of caving materials expands. In other words, the pile height of the caving rock mass will be larger than the original height. The volume expansion of the caving materials can be fully described with the bulking factorKp. The intact rock original height is Σh, and the rock mass pile height after roof strata collapse isKpΣh. The stratum void height Δ between the piled rock and the overlying stable rock will be(1)Δ=∑h+M−Kp∑h=M−∑hKp−1,where M denotes the coal seam height andKp denotes the bulking factor.According to equation (1), when Δ = 0, the cutting height of the roof strata is calculated by the following equation [15, 16]:(2)∑h=MKp−1.Based on a number of field observations, the caving height is 4 to 8 times of the coal seam height [17, 18]. For panel 8102, the coal seam height is 4 m, and the bulking factor is determined to be 1.25 [19, 20]. Therefore, the cutting height of the strata above the gob is 16 m.


## 4. Establishment of Numerical Model

### 4.1. Modeling Scheme

To analyse the relationship between the roadway stability and protective pillar width, a numerical model was established based on FLAC3D software. The dimensions of the model were 160 m × 59.1 m × 2 m. The model consisted of panels 8101 and 8102 and their roadway system (Figure9). A vertical stress of 7.5 MPa, representing the loads of overlying rock, is applied to the model top boundary. Based on the in situ stress monitoring data, the ratio of horizontal to vertical stress was determined to be 1.1. For the four vertical model planes, the horizontal displacement was constrained. The displacements in the horizontal and vertical directions were restrained in bottom boundary. The Mohr–Coulomb model was adopted to simulate the rock mass, except for the protective pillar and gob caving materials. Table 1 shows the mechanical properties of the rock mass.Figure 9
Numerical simulation model.Table 1
Mechanical properties of the rock mass.
LithologyDensity (kg/m3)Compressive strength (MPa)Elasticity modulus (GPa)Poisson's ration (ν)Cohesion (MPa)Friction (°)Fine sandstone275057.212.50.222.629Mudstone190013.15.00.291.425Sandy mudstone245027.87.10.261.827Siltstone268047.69.80.242.331Mudstone190013.15.00.291.425Coal seam16009.81.50.320.822Mudstone190013.15.00.291.425Fine sandstone275057.212.50.222.629Siltstone268047.69.80.242.331The numerical simulation comprised five steps: (1) model establishment and application of in situ stress, (2) initial stress balance in the model, (3) 8101 tailgate excavation, (4) retreat of panel 8101 and equivalent gob materials simulation, and (5) 8102 headgate excavation with five different pillar widths. On the basis of the 8102 headgate conditions, the protective pillar height was maintained at 3.5 m, and the protective pillar width was simulated at 5, 10, 15, 20, and 25 m, as shown in Figure10.Figure 10
The FLAC3D numerical simulation steps.


### 4.2. The Constitutive Model of Protective Pillars

There are elastic phase, plastic softening phase, and residual phase in the process of protective pillar failure. With pillar yielding, plastic softening phases occur until a residual strength level is obtained. At present, the widely used constitutive model of simulating protective pillar failure is the Mohr–Coulomb strain-softening model, in which pillars are considered to be materials with softening properties [21, 22]. After yielding, the friction angle and cohesion soften with the change in plastic strain. To simulate the strain-softening behaviour of protective pillar, a pillar submodel is established and includes the roof, pillar, and floor, as shown in Figure 11. The mechanical behaviour of the pillar is simulated by the strain-softening model, and the mechanical behaviour of the roof and floor is simulated by the Mohr–Coulomb model.Figure 11
A pillar submodel.Many strength formulas of pillar have been put forward in the recent decades [23–25]. It is noted that Salamon-Munro proposed the following empirical formula, with pillar W/H (width-to-height) ratios from 2 to 20 [23].(3)Pillar Strength=7.716w0.46h0.66,where w and h indicate the width and height of pillar, respectively.The empirical strength formula of the pillar has good consistency with developed average strength formulas. Therefore, the Salamon-Munro strength formula was adopted in this study to verify the pillar numerical model. A commonly used iterative method is adopted [21, 26]. The parameters of the strain-softening model are determined by matching the strength of the pillar acquired by numerical simulation with that determined by the Salamon-Munro formula. Table 2 shows the input parameters for calibrating the pillar strain-softening model. The simulated and calculated strengths of pillar are displayed in Figure 12. The results show that the simulated strength of the pillar matches with the calculated strength of pillar from the Salamon-Munro empirical formula very well, suggesting that the parameters can be adopted to accurately simulate the mechanical behaviour of a protective pillar.Table 2
Strain-softening properties of pillar with plastic strain.
Strain00.00250.0050.00750.01Cohesion (MPa)0.80.680.540.40.28Friction angle (°)2423222121Figure 12
The simulated and calculated strengths of pillar.


### 4.3. Modeling of the Caving Materials in the Gob

After panel retreat, the roof strata behind the panel collapse, and the caving rock in the gob is compacted and consolidated. Regarding the influence of the gob caving rock on the protective pillars, the mechanical behaviour of the gob caving materials was simulated by the double-yield model [21, 27]. The input parameters of the double-yield model included material properties and cap pressure [28]. The cap pressure parameters were determined by the following equation proposed by Salamon:(4)σ=E0ε1−ε/εm,where σ denotes the vertical stress applied to the gob caving materials; ε denotes the volumetric strain and maximum strainεm, respectively; E0 denotes the initial modulus of the gob materials. Parameters E0 and εm could be obtained by the two following equations:(5)εm=Kp−1Kp,(6)E0=10.39σc1.042Kp7.7,Where Kp denotes the bulking factor; and σc denotes the compressive strength of the caving materials in gob. Based on equation (2), the bulking factor Kp and compressive strength σc were determined to be 1.25 and 12 MPa, respectively. According to equations (5) and (6), the maximum strain εm and initial tangent modulus E0 are 0.2 and 24.8 MPa, respectively. The cap pressure parameters of the numerical model are shown in Table 3.Table 3
The double-yield model cap pressures.
Strain00.010.020.030.040.050.060.070.080.09Stress (MPa)00.260.550.881.241.652.132.673.314.06Strain0.100.110.120.130.140.150.160.170.180.19Stress (MPa)4.966.067.449.2111.5714.8819.8428.1044.6494.24To obtain the material parameters of the double-yield model (bulk modulus, shear modulus, friction angle, and dilation angle), a gob element (Figure13) was selected. A 10−5 m/s vertical velocity was used to represent the load and apply to the element top. The horizontal deformation of four vertical faces of the element was restrained, while the vertical deformation at the element bottom was set to zero. A fitting method in previous related research was adopted to match the gob element stress-strain curve with Salamon’s model [21, 29]. Figure 14 compares the stress-strain curves of two models, and they match very well. Based on the comparison result, the final parameters of the gob material in the double-yield model were obtained and used in FLAC3D (Table 4).Figure 13
An element of the gob.Figure 14
The stress-strain curve comparison of simulated and calculated model.Table 4
The caving material parameters in gob.
ParametersDensity (kg/m3)Bulk modulus (GPa)Shear modulus (GPa)Friction (°)Dilation (°)Value12004.82.7255


### 4.4. The Relationship between the Stress and Energy Density of the Surrounding Rock Element

The stress state of the surrounding rock structure is closely related to its stored energy. During excavation, the stress of the rock mass is adjusted with changes of the energy storage structure. As the essential attribute of the deformation and destruction of surrounding rock, the energy change reflects the continuous development of defects, weakening, and loss of strength in the surrounding rock mass [30]. Therefore, failure of the surrounding rock structure is reflected by internal energy change.An element of the surrounding rock structure was chosen, and the size and shape of the element changed with the external forces. Figure15 shows the 3D stress state of the element, which can be considered the superposition of two stress states. On the one hand, it can be described by the average stress σm of the element in three directions:(7)σm=σ1+σ2+σ33,where σ1, σ2, and σ3 are the maximum, intermediate, and minimum principal stress, respectively.Figure 15
Three-dimensional stress state of surrounding rock element.On the other hand, it can be considered that stressesσ1−=σ1−σm, σ2−=σ2−σm, andσ3−=σ3−σm are applied in the three directions of the element, where σ1−, σ2−, and σ3− are considered to be the deviations in stress from the principal stresses σ1, σ2, and σ3 [31].(8)σ1−=2σ1−σ2−σ33,σ2−=2σ2−σ3−σ13,σ3−=2σ3−σ1−σ23.Under the action of the average stressσm, the volume of the element changes, while its shape is invariant. The strain energy density was calculated according to the following formula:(9)υv=1−2μ6Eσ1+σ2+σ32,where υv denotes the strain energy density under the action of average stressσm; and σ1, σ2, and σ3 are the maximum, intermediate, and minimum principal stress, respectively.Under the action of the deviatoric stressesσ1−, σ2−, and σ3−, the shape of the element changes, while its volume is invariant. Thus, the strain energy density was calculated according to the following formula:(10)υd=1+μ6Eσ1−σ22+σ2−σ32+σ3−σ12,where υd denotes strain energy density under the action of the deviatoric stresses σ1−, σ2−, and σ3−.The element strain energy density includes the volume-changed and shape-changed strain energy density. Therefore, the strain energy density of the element can be expressed [32, 33]:(11)υ=υv+υd=1−2μ6Eσ1+σ2+σ32+1+μ6Eσ1−σ22+σ2−σ32+σ3−σ12,where υ denotes the sum of the strain energy density caused by the changes in volume and shape of the element.According to equation (11), it can be clearly seen that energy stored in the surrounding rock element is closely related to its own stress state. A coupled program of the strain energy density and stress of a surrounding rock element was developed in this study, which was used to investigate the relationship between the internal energy and stress of the rock mass.


### 4.5. Verification of the Gob Caving Materials Supporting Capacity

To validate the effectiveness of the supporting capacity of caving materials in gob 8101, Figure16 displays the vertical stress evolution law in the caving area. There is 0.42 MPa of vertical stress at the gob edge. Meanwhile, the vertical stress increases gradually with distance increase from the gob edge. The maximal vertical stress reaches 7.45 MPa at an approximately 73.6 m distance from the gob edge and remains relatively steady. The caving materials in the gob bear 99% of the initial vertical stress (7.45/7.5 MPa). With increasing distance from the gob edge, the caving materials are compressed and the vertical stress in the caving area increases gradually. The vertical stress of the gob reaches the initial stress (7.5 MPa) at approximately 73.6 m from the gob edge, that is, 25% (73.6/300) of the overlying depth. Through extensive field tests, Wilson and Carr [34] and Campoli et al. [35] proposed that the gob vertical stress could reach the initial stress of the rock mass at a distance of 0.2–0.38 times the buried depth of coal seam. In this study, the obtained vertical stress distribution agrees well with the related conclusions [14, 25, 36], indicating that the simulation result of the supporting capacity of the gob caving materials is reliable.Figure 16
The vertical stress distribution in gob 8101.


## 4.1. Modeling Scheme

To analyse the relationship between the roadway stability and protective pillar width, a numerical model was established based on FLAC3D software. The dimensions of the model were 160 m × 59.1 m × 2 m. The model consisted of panels 8101 and 8102 and their roadway system (Figure9). A vertical stress of 7.5 MPa, representing the loads of overlying rock, is applied to the model top boundary. Based on the in situ stress monitoring data, the ratio of horizontal to vertical stress was determined to be 1.1. For the four vertical model planes, the horizontal displacement was constrained. The displacements in the horizontal and vertical directions were restrained in bottom boundary. The Mohr–Coulomb model was adopted to simulate the rock mass, except for the protective pillar and gob caving materials. Table 1 shows the mechanical properties of the rock mass.Figure 9
Numerical simulation model.Table 1
Mechanical properties of the rock mass.
LithologyDensity (kg/m3)Compressive strength (MPa)Elasticity modulus (GPa)Poisson's ration (ν)Cohesion (MPa)Friction (°)Fine sandstone275057.212.50.222.629Mudstone190013.15.00.291.425Sandy mudstone245027.87.10.261.827Siltstone268047.69.80.242.331Mudstone190013.15.00.291.425Coal seam16009.81.50.320.822Mudstone190013.15.00.291.425Fine sandstone275057.212.50.222.629Siltstone268047.69.80.242.331The numerical simulation comprised five steps: (1) model establishment and application of in situ stress, (2) initial stress balance in the model, (3) 8101 tailgate excavation, (4) retreat of panel 8101 and equivalent gob materials simulation, and (5) 8102 headgate excavation with five different pillar widths. On the basis of the 8102 headgate conditions, the protective pillar height was maintained at 3.5 m, and the protective pillar width was simulated at 5, 10, 15, 20, and 25 m, as shown in Figure10.Figure 10
The FLAC3D numerical simulation steps.


## 4.2. The Constitutive Model of Protective Pillars

There are elastic phase, plastic softening phase, and residual phase in the process of protective pillar failure. With pillar yielding, plastic softening phases occur until a residual strength level is obtained. At present, the widely used constitutive model of simulating protective pillar failure is the Mohr–Coulomb strain-softening model, in which pillars are considered to be materials with softening properties [21, 22]. After yielding, the friction angle and cohesion soften with the change in plastic strain. To simulate the strain-softening behaviour of protective pillar, a pillar submodel is established and includes the roof, pillar, and floor, as shown in Figure 11. The mechanical behaviour of the pillar is simulated by the strain-softening model, and the mechanical behaviour of the roof and floor is simulated by the Mohr–Coulomb model.Figure 11
A pillar submodel.Many strength formulas of pillar have been put forward in the recent decades [23–25]. It is noted that Salamon-Munro proposed the following empirical formula, with pillar W/H (width-to-height) ratios from 2 to 20 [23].(3)Pillar Strength=7.716w0.46h0.66,where w and h indicate the width and height of pillar, respectively.The empirical strength formula of the pillar has good consistency with developed average strength formulas. Therefore, the Salamon-Munro strength formula was adopted in this study to verify the pillar numerical model. A commonly used iterative method is adopted [21, 26]. The parameters of the strain-softening model are determined by matching the strength of the pillar acquired by numerical simulation with that determined by the Salamon-Munro formula. Table 2 shows the input parameters for calibrating the pillar strain-softening model. The simulated and calculated strengths of pillar are displayed in Figure 12. The results show that the simulated strength of the pillar matches with the calculated strength of pillar from the Salamon-Munro empirical formula very well, suggesting that the parameters can be adopted to accurately simulate the mechanical behaviour of a protective pillar.Table 2
Strain-softening properties of pillar with plastic strain.
Strain00.00250.0050.00750.01Cohesion (MPa)0.80.680.540.40.28Friction angle (°)2423222121Figure 12
The simulated and calculated strengths of pillar.


## 4.3. Modeling of the Caving Materials in the Gob

After panel retreat, the roof strata behind the panel collapse, and the caving rock in the gob is compacted and consolidated. Regarding the influence of the gob caving rock on the protective pillars, the mechanical behaviour of the gob caving materials was simulated by the double-yield model [21, 27]. The input parameters of the double-yield model included material properties and cap pressure [28]. The cap pressure parameters were determined by the following equation proposed by Salamon:(4)σ=E0ε1−ε/εm,where σ denotes the vertical stress applied to the gob caving materials; ε denotes the volumetric strain and maximum strainεm, respectively; E0 denotes the initial modulus of the gob materials. Parameters E0 and εm could be obtained by the two following equations:(5)εm=Kp−1Kp,(6)E0=10.39σc1.042Kp7.7,Where Kp denotes the bulking factor; and σc denotes the compressive strength of the caving materials in gob. Based on equation (2), the bulking factor Kp and compressive strength σc were determined to be 1.25 and 12 MPa, respectively. According to equations (5) and (6), the maximum strain εm and initial tangent modulus E0 are 0.2 and 24.8 MPa, respectively. The cap pressure parameters of the numerical model are shown in Table 3.Table 3
The double-yield model cap pressures.
Strain00.010.020.030.040.050.060.070.080.09Stress (MPa)00.260.550.881.241.652.132.673.314.06Strain0.100.110.120.130.140.150.160.170.180.19Stress (MPa)4.966.067.449.2111.5714.8819.8428.1044.6494.24To obtain the material parameters of the double-yield model (bulk modulus, shear modulus, friction angle, and dilation angle), a gob element (Figure13) was selected. A 10−5 m/s vertical velocity was used to represent the load and apply to the element top. The horizontal deformation of four vertical faces of the element was restrained, while the vertical deformation at the element bottom was set to zero. A fitting method in previous related research was adopted to match the gob element stress-strain curve with Salamon’s model [21, 29]. Figure 14 compares the stress-strain curves of two models, and they match very well. Based on the comparison result, the final parameters of the gob material in the double-yield model were obtained and used in FLAC3D (Table 4).Figure 13
An element of the gob.Figure 14
The stress-strain curve comparison of simulated and calculated model.Table 4
The caving material parameters in gob.
ParametersDensity (kg/m3)Bulk modulus (GPa)Shear modulus (GPa)Friction (°)Dilation (°)Value12004.82.7255


## 4.4. The Relationship between the Stress and Energy Density of the Surrounding Rock Element

The stress state of the surrounding rock structure is closely related to its stored energy. During excavation, the stress of the rock mass is adjusted with changes of the energy storage structure. As the essential attribute of the deformation and destruction of surrounding rock, the energy change reflects the continuous development of defects, weakening, and loss of strength in the surrounding rock mass [30]. Therefore, failure of the surrounding rock structure is reflected by internal energy change.An element of the surrounding rock structure was chosen, and the size and shape of the element changed with the external forces. Figure15 shows the 3D stress state of the element, which can be considered the superposition of two stress states. On the one hand, it can be described by the average stress σm of the element in three directions:(7)σm=σ1+σ2+σ33,where σ1, σ2, and σ3 are the maximum, intermediate, and minimum principal stress, respectively.Figure 15
Three-dimensional stress state of surrounding rock element.On the other hand, it can be considered that stressesσ1−=σ1−σm, σ2−=σ2−σm, andσ3−=σ3−σm are applied in the three directions of the element, where σ1−, σ2−, and σ3− are considered to be the deviations in stress from the principal stresses σ1, σ2, and σ3 [31].(8)σ1−=2σ1−σ2−σ33,σ2−=2σ2−σ3−σ13,σ3−=2σ3−σ1−σ23.Under the action of the average stressσm, the volume of the element changes, while its shape is invariant. The strain energy density was calculated according to the following formula:(9)υv=1−2μ6Eσ1+σ2+σ32,where υv denotes the strain energy density under the action of average stressσm; and σ1, σ2, and σ3 are the maximum, intermediate, and minimum principal stress, respectively.Under the action of the deviatoric stressesσ1−, σ2−, and σ3−, the shape of the element changes, while its volume is invariant. Thus, the strain energy density was calculated according to the following formula:(10)υd=1+μ6Eσ1−σ22+σ2−σ32+σ3−σ12,where υd denotes strain energy density under the action of the deviatoric stresses σ1−, σ2−, and σ3−.The element strain energy density includes the volume-changed and shape-changed strain energy density. Therefore, the strain energy density of the element can be expressed [32, 33]:(11)υ=υv+υd=1−2μ6Eσ1+σ2+σ32+1+μ6Eσ1−σ22+σ2−σ32+σ3−σ12,where υ denotes the sum of the strain energy density caused by the changes in volume and shape of the element.According to equation (11), it can be clearly seen that energy stored in the surrounding rock element is closely related to its own stress state. A coupled program of the strain energy density and stress of a surrounding rock element was developed in this study, which was used to investigate the relationship between the internal energy and stress of the rock mass.


## 4.5. Verification of the Gob Caving Materials Supporting Capacity

To validate the effectiveness of the supporting capacity of caving materials in gob 8101, Figure16 displays the vertical stress evolution law in the caving area. There is 0.42 MPa of vertical stress at the gob edge. Meanwhile, the vertical stress increases gradually with distance increase from the gob edge. The maximal vertical stress reaches 7.45 MPa at an approximately 73.6 m distance from the gob edge and remains relatively steady. The caving materials in the gob bear 99% of the initial vertical stress (7.45/7.5 MPa). With increasing distance from the gob edge, the caving materials are compressed and the vertical stress in the caving area increases gradually. The vertical stress of the gob reaches the initial stress (7.5 MPa) at approximately 73.6 m from the gob edge, that is, 25% (73.6/300) of the overlying depth. Through extensive field tests, Wilson and Carr [34] and Campoli et al. [35] proposed that the gob vertical stress could reach the initial stress of the rock mass at a distance of 0.2–0.38 times the buried depth of coal seam. In this study, the obtained vertical stress distribution agrees well with the related conclusions [14, 25, 36], indicating that the simulation result of the supporting capacity of the gob caving materials is reliable.Figure 16
The vertical stress distribution in gob 8101.


## 5. Modeling Results and Discussion

### 5.1. Modeling Results

The simulation results of the vertical stress and energy density of the two ribs in the 8102 headgate are shown in Figure17. To analyse the simulation results in detail, Figure 17(a) is taken as an example. The upper diagram indicates vertical stress distribution and the lower diagram shows the energy density distribution of two roadway ribs. The blue and yellow areas in the middle represent elements of surrounding rock in the elastic and yielded states, respectively. The curves on the two sides of the roadway represent the changes of vertical stress and energy density. Note that all the data come from the mid-height on the two roadway ribs. The peak values of the vertical stress and energy density in the protective pillar are designated as σcs and vcd, respectively, while the peak values of the vertical stress and energy density in the virgin coal rib are designated as σvs and vvd, respectively.Figure 17
Vertical stress and energy density distribution on two ribs of 8102 headgate with different pillar widths. (a) The width of 5 m. (b) The width of 10 m. (c) The width of 15 m. (d) The width of 20 m. (e) The width of 25 m.
(a)(b)(c)(d)(e)With a 5 m width, the protective pillar is damaged completely and could not keep the stability of the 8102 headgate, as shown in Figure17(a). Stress concentration and energy accumulation occur in a position 4.6 m from the edge of the protective pillar rib, with σcs of 8.6 MPa and vcd of 176.1 kJ/m3, while at a position 14.3 m from the virgin rib edge, σvs and vvd are 16.7 MPa and 461.5 kJ/m3, respectively. The peak values of the vertical stress and energy density for the virgin coal rib are obviously greater than those of the protective pillar. When the pillar width reaches 10 m, protective pillar is in yield state but is not crushed. The virgin coal rib yield range decreases. Moreover, σcs increases gradually from 8.6 MPa to 9.5 MPa with vcd increasing from 176.1 kJ/m3 to 208.3 kJ/m3, remaining lower than the vertical stress (16.2 MPa) and energy density (434.8 kJ/m3) of the virgin coal rib (Figure 17(b)), showing that the roof strata loads are mainly on the solid coal region. When the protective pillar is 15 m and 20 m wide, the vertical stress and energy density increase sharply and exceed the virgin rib (Figures 17(c) and 17(d)), showing that mining pressure is transferred rapidly to the protective coal pillar, sharply changing the energy density. The excess stress and elastic energy result in severe roadway deformation and the occurrence of coal bumps, which are identical to the failure characteristics of 8102 headgate in the field test. With condition of 25 m wide protective pillar, there is a 6 m elastic zone in the pillar (Figure 17(e)). The yielded range of the two roadway ribs shrinks significantly. A stable elastic energy storage structure is formed in the protective pillar, which is capable of absorbing substantial stress and elastic energy. The distribution of the vertical stress and energy density exhibits a double peak, which is consistent with the conclusion based on Figure 7(b). The modeling results illustrate that adjusting the pillar width improves the stability of the surrounding rock.For the different protective pillar widths, the roadway deformation laws are shown in Figure18. The displacement of two ribs in roadway is obviously larger than that of the roof and floor. The effect of the protective pillar width on the stability of the two ribs is greater than that of the roadway roof and floor. Meanwhile, the roadway deformation of 10 m pillar width is similar to that of 25 m wide protective pillar.Figure 18
Roadway deformation with different pillar width.The results (Figures17 and 18) show that the protective pillar width influences greatly the pillar load-bearing capacity and roadway deformation. Most surprisingly, when the protective pillar widths are 10 m and 25 m, the difference in roadway deformation is not obvious. This demonstrates that the existence of an elastic intact zone enhances the stability of the protective pillar. Consequently, increasing the protective pillar width could enhance the reliability of the pillar load-bearing capacity and improve the roadway stability.


### 5.2. Discussion of the Rational Width for a Protective Pillar

The original balanced state of the main roof above the roadway is strongly disturbed during the panel retreat period. The main roof is broken [37–39], as shown in Figure 19(a). Rock masses A, B, and C interact, which influences the stability of the 8102 headgate. During the roof fracture period, rock mass B above the roadway begins to rotate and subside. The rotation and subsidence have a great impact on the stability of the protective pillar. The rock B fracture position and rotating speed are closely related to the pillar width. According to the ultimate balance theory, the geometry size l of rock B can be derived from the panel dip length S and periodic breakage length L of the main roof, as in the following formula [40, 41]:(12)l=2L1710LS2+102−10LS,where S is the dip length of the panel and L is the periodic breakage length of the main roof. From field measurements, the periodic breakage length of the main roof is 23 m. The dip length of panel 8102 was 210 m. From equation (12), the geometry size lof rock B was 24.5 m.Figure 19
Failure modes of the protective pillar. (a) 5 m wide. (b) 20 m wide.
(a)(b)From the results of the numerical simulation, in the case of a 5 m wide protective pillar, the peak vertical stress of the virgin coal rib occurs 14.3 m from the rib edge. The length between the peak vertical stress location in the virgin rib and the rib edge plus the width of the protective pillar and the roadway equals the geometry size of rock B, indicating that rock mass B will fracture at the location where the peak vertical stress of the virgin coal rib occurs (Figure19(a)). A 5 m wide protective pillar cannot withstand the loads of the overlying strata. Rock mass B rotates violently and the protective pillar is crushed. With a 20 m wide protective pillar, the vertical stress peak in the virgin rib is transferred to the pillar. Based on the geometry size of rock B, it could be concluded that its fracture position is located at the junction of the roof and the virgin coal rib (Figure 19(b)), which is consistent with the field result showing that roof sagging and falling often occur at the place close to the junction between the roof and the virgin coal rib. In the panel retreat period, the superimposed effect of mining stress and high in situ stress surpasses the bearing strength of the protective pillar, and the continued pillar deformation leads to the fracture and rotation movement of rock mass B in the main roof. Finally, the protective pillar is crushed and fails. With a 10 m wide protective pillar, the vertical stress peak is still located in the virgin coal rib; that is, the area of virgin rib is subjected to the main load of the overlying strata. The pillar bears relatively less load, which puts the roadway in a state of low stress. Meanwhile, the protective pillar could provide support to the fractured rock mass B and prevent it from rotary movement, which maintains the stability of the roadway. When the width of protective pillar is 25 m, the bearing capacity is improved due to the expansion of the elastic area in the pillar. The protective pillar can not only maintain its own stability but also prevent the fractured rock mass B of the main roof from rotating, transitioning the roadway to a stable state. However, a greater pillar width wastes more coal resources.A rational pillar width should have the capacity to withstand abutment pressure and maintain a stable roadway, while maximizing the recovery of coal resources. Based on the above discussion, the recommended 10 m wide protective pillar may be more conducive to improving the roadway stability and coal recovery rate, while reducing economic loss due to large deformation and coal bump failure.


## 5.1. Modeling Results

The simulation results of the vertical stress and energy density of the two ribs in the 8102 headgate are shown in Figure17. To analyse the simulation results in detail, Figure 17(a) is taken as an example. The upper diagram indicates vertical stress distribution and the lower diagram shows the energy density distribution of two roadway ribs. The blue and yellow areas in the middle represent elements of surrounding rock in the elastic and yielded states, respectively. The curves on the two sides of the roadway represent the changes of vertical stress and energy density. Note that all the data come from the mid-height on the two roadway ribs. The peak values of the vertical stress and energy density in the protective pillar are designated as σcs and vcd, respectively, while the peak values of the vertical stress and energy density in the virgin coal rib are designated as σvs and vvd, respectively.Figure 17
Vertical stress and energy density distribution on two ribs of 8102 headgate with different pillar widths. (a) The width of 5 m. (b) The width of 10 m. (c) The width of 15 m. (d) The width of 20 m. (e) The width of 25 m.
(a)(b)(c)(d)(e)With a 5 m width, the protective pillar is damaged completely and could not keep the stability of the 8102 headgate, as shown in Figure17(a). Stress concentration and energy accumulation occur in a position 4.6 m from the edge of the protective pillar rib, with σcs of 8.6 MPa and vcd of 176.1 kJ/m3, while at a position 14.3 m from the virgin rib edge, σvs and vvd are 16.7 MPa and 461.5 kJ/m3, respectively. The peak values of the vertical stress and energy density for the virgin coal rib are obviously greater than those of the protective pillar. When the pillar width reaches 10 m, protective pillar is in yield state but is not crushed. The virgin coal rib yield range decreases. Moreover, σcs increases gradually from 8.6 MPa to 9.5 MPa with vcd increasing from 176.1 kJ/m3 to 208.3 kJ/m3, remaining lower than the vertical stress (16.2 MPa) and energy density (434.8 kJ/m3) of the virgin coal rib (Figure 17(b)), showing that the roof strata loads are mainly on the solid coal region. When the protective pillar is 15 m and 20 m wide, the vertical stress and energy density increase sharply and exceed the virgin rib (Figures 17(c) and 17(d)), showing that mining pressure is transferred rapidly to the protective coal pillar, sharply changing the energy density. The excess stress and elastic energy result in severe roadway deformation and the occurrence of coal bumps, which are identical to the failure characteristics of 8102 headgate in the field test. With condition of 25 m wide protective pillar, there is a 6 m elastic zone in the pillar (Figure 17(e)). The yielded range of the two roadway ribs shrinks significantly. A stable elastic energy storage structure is formed in the protective pillar, which is capable of absorbing substantial stress and elastic energy. The distribution of the vertical stress and energy density exhibits a double peak, which is consistent with the conclusion based on Figure 7(b). The modeling results illustrate that adjusting the pillar width improves the stability of the surrounding rock.For the different protective pillar widths, the roadway deformation laws are shown in Figure18. The displacement of two ribs in roadway is obviously larger than that of the roof and floor. The effect of the protective pillar width on the stability of the two ribs is greater than that of the roadway roof and floor. Meanwhile, the roadway deformation of 10 m pillar width is similar to that of 25 m wide protective pillar.Figure 18
Roadway deformation with different pillar width.The results (Figures17 and 18) show that the protective pillar width influences greatly the pillar load-bearing capacity and roadway deformation. Most surprisingly, when the protective pillar widths are 10 m and 25 m, the difference in roadway deformation is not obvious. This demonstrates that the existence of an elastic intact zone enhances the stability of the protective pillar. Consequently, increasing the protective pillar width could enhance the reliability of the pillar load-bearing capacity and improve the roadway stability.


## 5.2. Discussion of the Rational Width for a Protective Pillar

The original balanced state of the main roof above the roadway is strongly disturbed during the panel retreat period. The main roof is broken [37–39], as shown in Figure 19(a). Rock masses A, B, and C interact, which influences the stability of the 8102 headgate. During the roof fracture period, rock mass B above the roadway begins to rotate and subside. The rotation and subsidence have a great impact on the stability of the protective pillar. The rock B fracture position and rotating speed are closely related to the pillar width. According to the ultimate balance theory, the geometry size l of rock B can be derived from the panel dip length S and periodic breakage length L of the main roof, as in the following formula [40, 41]:(12)l=2L1710LS2+102−10LS,where S is the dip length of the panel and L is the periodic breakage length of the main roof. From field measurements, the periodic breakage length of the main roof is 23 m. The dip length of panel 8102 was 210 m. From equation (12), the geometry size lof rock B was 24.5 m.Figure 19
Failure modes of the protective pillar. (a) 5 m wide. (b) 20 m wide.
(a)(b)From the results of the numerical simulation, in the case of a 5 m wide protective pillar, the peak vertical stress of the virgin coal rib occurs 14.3 m from the rib edge. The length between the peak vertical stress location in the virgin rib and the rib edge plus the width of the protective pillar and the roadway equals the geometry size of rock B, indicating that rock mass B will fracture at the location where the peak vertical stress of the virgin coal rib occurs (Figure19(a)). A 5 m wide protective pillar cannot withstand the loads of the overlying strata. Rock mass B rotates violently and the protective pillar is crushed. With a 20 m wide protective pillar, the vertical stress peak in the virgin rib is transferred to the pillar. Based on the geometry size of rock B, it could be concluded that its fracture position is located at the junction of the roof and the virgin coal rib (Figure 19(b)), which is consistent with the field result showing that roof sagging and falling often occur at the place close to the junction between the roof and the virgin coal rib. In the panel retreat period, the superimposed effect of mining stress and high in situ stress surpasses the bearing strength of the protective pillar, and the continued pillar deformation leads to the fracture and rotation movement of rock mass B in the main roof. Finally, the protective pillar is crushed and fails. With a 10 m wide protective pillar, the vertical stress peak is still located in the virgin coal rib; that is, the area of virgin rib is subjected to the main load of the overlying strata. The pillar bears relatively less load, which puts the roadway in a state of low stress. Meanwhile, the protective pillar could provide support to the fractured rock mass B and prevent it from rotary movement, which maintains the stability of the roadway. When the width of protective pillar is 25 m, the bearing capacity is improved due to the expansion of the elastic area in the pillar. The protective pillar can not only maintain its own stability but also prevent the fractured rock mass B of the main roof from rotating, transitioning the roadway to a stable state. However, a greater pillar width wastes more coal resources.A rational pillar width should have the capacity to withstand abutment pressure and maintain a stable roadway, while maximizing the recovery of coal resources. Based on the above discussion, the recommended 10 m wide protective pillar may be more conducive to improving the roadway stability and coal recovery rate, while reducing economic loss due to large deformation and coal bump failure.


## 6. Field Test

A protective pillar width of 10 m was applied in field to maintain the stability of the surrounding rock 8103 headgate. The protective pillar stress distribution and deformation of the 8103 headgate were measured to validate the reliability of the simulation results.


### 6.1. Stress Measurement of the Protective Pillar in the Field

The vertical stress at different distances from the protective pillar rib of the 8103 headgate is plotted in Figure20. The field data are also displayed for comparison with the modeling results. The measured results and simulated results have good agreement, verifying the numerical modeling accuracy.Figure 20
Stress measurement of protective pillar. (a) Layout of stress meters. (b) Distribution law of vertical stress after roadway excavation.
(a)(b)


### 6.2. Roadway Deformation Laws

The measured roadway deformation results at different stages are shown in Figure21. The deformation of the roadway became stable after 50 days of excavation. The convergences of the roof, virgin coal rib, protective pillar rib, and floor were 115 mm, 74 mm, 87 mm, and 16 mm, respectively. During the panel 8103 retreat period, the roadway deformation mainly occurred 60 m in front of the working face; the total convergences at the roof, virgin coal rib, protective pillar rib, and floor were 301 mm, 162 mm, 214 mm, and 24 mm, respectively, which represented reductions of 70%, 59%, 68%, and 77% (Figure 22), compared with the deformation that occurred under the previous protective pillar width. The stability of the surrounding 8103 headgate is improved, as shown in Figure 23. The field observation results indicate that the proposed protective pillar width is feasible for maintaining the stability of the roadway.Figure 21
Measured convergences in 8103 headgate. (a) During 8103 headgate excavation. (b) During panel 8103 retreat.
(a)(b)Figure 22
Comparison of roadway convergences during panel 8103 retreat.Figure 23
The stable environment of surrounding rock in 8103 headgate.


## 6.1. Stress Measurement of the Protective Pillar in the Field

The vertical stress at different distances from the protective pillar rib of the 8103 headgate is plotted in Figure20. The field data are also displayed for comparison with the modeling results. The measured results and simulated results have good agreement, verifying the numerical modeling accuracy.Figure 20
Stress measurement of protective pillar. (a) Layout of stress meters. (b) Distribution law of vertical stress after roadway excavation.
(a)(b)


## 6.2. Roadway Deformation Laws

The measured roadway deformation results at different stages are shown in Figure21. The deformation of the roadway became stable after 50 days of excavation. The convergences of the roof, virgin coal rib, protective pillar rib, and floor were 115 mm, 74 mm, 87 mm, and 16 mm, respectively. During the panel 8103 retreat period, the roadway deformation mainly occurred 60 m in front of the working face; the total convergences at the roof, virgin coal rib, protective pillar rib, and floor were 301 mm, 162 mm, 214 mm, and 24 mm, respectively, which represented reductions of 70%, 59%, 68%, and 77% (Figure 22), compared with the deformation that occurred under the previous protective pillar width. The stability of the surrounding 8103 headgate is improved, as shown in Figure 23. The field observation results indicate that the proposed protective pillar width is feasible for maintaining the stability of the roadway.Figure 21
Measured convergences in 8103 headgate. (a) During 8103 headgate excavation. (b) During panel 8103 retreat.
(a)(b)Figure 22
Comparison of roadway convergences during panel 8103 retreat.Figure 23
The stable environment of surrounding rock in 8103 headgate.


## 7. Conclusions

This study was mainly focused on analysing the effect of protective pillar width on the roadway stability to identify a pillar design principle based on field tests and numerical simulations. The main conclusions in this paper are summarized as follows:(1)
A novel numerical model was established to analyse the failure mechanism of the protective pillar. As an innovative design method of protective pillars, the supporting features of the gob caving materials on overlying strata and the relationship between the internal stress and stored energy of the rock mass were considered in the modeling process. The stress change and distribution characteristics of the energy density are regarded as important conditions in designing the width of a protective pillar and evaluating the stability of the roadway.(2)
The modeling results showed that, with a 20 m wide pillar, the peak vertical stress and energy density in the protective pillar were 18.5 MPa and 563.7 kJ/m3, respectively. Excessive stress and elastic energy resulted in considerable deformation and coal bump failure in headgate 8102. With a pillar width of 10 m, the location of the peak vertical stress and energy density moved from the protective pillar to the virgin rib. The main loads of the overlaying strata were borne by the virgin rib, and the pillar was subjected to a relatively low load. The roadway was in a state of low stress and could maintain its stability.(3)
The results of the field measurements showed that a 10 m wide protective pillar was able to effectively control the roadway deformation and release most of the storing energy in protective pillar. Meanwhile, the proposed modeling approach and protective pillar design principles in this study can provide a useful basis for application in similar coal mines.



---

*Source: 5533364-2021-02-28.xml*