Research on Catastrophic Pillar Instability in Room and Pillar Gypsum Mining

Gypsum mines in China are mostly exploited through room and pillar mining. Due to backward mining technology and a long history of mining, a great number of pillars were left in gypsum mines. Many serious work safety accidents occurred as the result of goaf instability in history, which posed severe threats to the security of people’s lives and property. Based on the characteristics of surrounding rock damage, this research improved the constitutive equation of gypsum rock mass damage by establishing a damage evolution model and introducing a shape parameter. Meanwhile, the cusp catastrophe equation was deduced based on the catastrophe theory and the constitutive equation of gypsum rock mass damage, thus summarizing the criteria for pillar instability; the pillar safety factor was obtained by means of the interrelation between pillar load and pillar strength. Based on the criteria for pillar instability and the pillar safety factor obtained, the necessary and sufficient conditions for pillar stability were concluded. These conclusions are of significance in that they provide theoretic reference for the treatment of gypsum goaf, as well as for further mining.


Introduction
Seventy percent of gypsum mines in China are exploited through underground mining [1,2], or to be specific, room and pillar mining.However, due to insufficient investment in safety measures, and backward mining technology and equipment, serious mining accidents occurred continually in recent years, among which the casualties in production areas caused by large-area collapse in the goaf ranked at the top.With the advancement of recovery, the exposure area becomes larger and larger, and the exposure duration longer and longer.Consequently, the ground pressure increases accordingly and leads to various ground pressure phenomena, such as local collapse in the goaf, pillar deformation, and difficulty in stope maintenance in the adjacent work zones.The locations of high stress (especially tensile stress) or in soft overburden rocks and unfavorable geological structure development are highly likely to become the breakthrough point for large-area sudden collapse and strata movement, thus causing devastating consequences to production and safety.When the number of mined areas exceeds the critical points, the untimely treatment often leads to rock collapse and strata dislocation with shock airflows and noises in the stope and roadway, which cause serious casualties and ground pressure disasters like surface cracking and damage, ground-building subsidence and collapse, and rock burst in the mining area.
The previous studies on room and pillar mining, be it focusing on the stability of stope or the optimization of stope structural parameters [3][4][5], took the pillar stability as the principal research object.Generally speaking, the pillar stability is studied from the following aspects: pillar strength, stress characteristics of pillars, failure mechanism of pillar instability, and the interaction between surrounding rock and pillars.Other methods related to rock damage were also used to study goaf pillars, such as numerical methods, finite element method (FEM), extended FEM (X-FEM), and phase field [6][7][8][9][10][11][12].These above-mentioned methods, however, cannot reflect the pillar damage from the nature of gypsum rock damage.Therefore, the constitutive law of gypsum damage was constructed in this paper to analyze the stability of goaf in combination with catastrophic theory (since catastrophic theory is the analysis of brittle rock damage, this paper does not consider the softening stage of gypsum).
Based on the mechanism of gypsum damage evolution and improved cusp catastrophe theory, this paper improves the equation of damage evolution and establishes the criterion for pillar instability in room and pillar gypsum mining.The aim of this paper was to master the failure laws of pillars and to lay a solid theoretic foundation for feasible solutions to maintain the pillar and even the goaf stability of gypsum mines, thus promoting the sustainable development of mining.

Evolution Equation and Constitutive Equation of Gypsum Rock Mass Damage
In the room and pillar mining process of gypsum, blasting mining and disturbance often cause certain damage to the rock mass of the pillar, thus reducing the effective bearing area of the pillar.As the mining activities advance, the effective bearing area of the pillar decreases correspondingly.When the decrease reaches a certain degree, sudden catastrophic instability occurs within the pillar.This kind of instability is characteristically nonlinear.To study the nonlinear instability process of a gypsum pillar, the damage degree of gypsum rock mass should be studied by firstly introducing a damage variable.To be specific, the variable in this paper refers to the damage factor D, which denotes structure and mechanical properties of the rocks (elastic modulus or strength) caused by crack growth.It is related to crack density and loading state; meanwhile, it is also under the impact of deformation and strain rate.Let the apparent area of the material be A and the actual effective loading area be A ef , and then the damage factor D can be expressed as follows: If there is no damage in the material before loading, D 0 = 0; if there are initial cracks in the material, D 0 = 0; if there is rupture failure in the material, D 0 = 1.Thus, D < 1 indicates that there is some damage of various degrees in the material.
The decrease in effective bearing area leads to the increase in true stress.
Based on the principle of strain equivalence proposed by Lemaitre [13], it is supposed that the strain caused by the stress σ on the damaged material is equivalent to the strain caused by the effective stress σ on the undamaged material, which can be expressed as follows(as shown in Figure 1): and then, where D denotes the damage variable, E refers to the elastic modulus of the undamaged material, and E is the elastic modulus of the damaged material.
Sustainability 2018, 10, x FOR PEER REVIEW 3 of 12 where D denotes the damage variable, E refers to the elastic modulus of the undamaged material, and E′ is the elastic modulus of the damaged material.In the mining process, the gypsum rock mass is damaged to various degrees due to the blasting and pressure release of the surrounding rock.This damage can be shown by the stress-strain curve of the rock mass.Therefore, this weakening behavior can serve as an indicator of the damage inside the gypsum rock mass.In the loading process, the thin fissures and micro-cracks inside the rock mass grow into larger ones, and finally cause the damage to the gypsum rock mass.A damage parameter D is set to represent the volume of cracks in the rock mass.Most random distributions in the natural world are consistent with the rules of Weibull random distribution, as are the micro-cracks inside the gypsum rock mass.In other words, the damage parameter D is consistent with Weibull random distribution, which can be expressed as follows: where  denotes the strain, and m refers to a shape parameter.In the above equation, m and a are nonnegative.
Based on the fundamental equation of continuum damage mechanics, it can be drawn that Before reaching the peak value of the stress-strain curve, there are four geometric conditions (as shown in Figure 2): denotes the peak-load strain value, and max  denotes the peak-load stress value.In the mining process, the gypsum rock mass is damaged to various degrees due to the blasting and pressure release of the surrounding rock.This damage can be shown by the stress-strain curve of the rock mass.Therefore, this weakening behavior can serve as an indicator of the damage inside the gypsum rock mass.In the loading process, the thin fissures and micro-cracks inside the rock mass grow into larger ones, and finally cause the damage to the gypsum rock mass.A damage parameter D is set to represent the volume of cracks in the rock mass.Most random distributions in the natural world are consistent with the rules of Weibull random distribution, as are the micro-cracks inside the gypsum rock mass.In other words, the damage parameter D is consistent with Weibull random distribution, which can be expressed as follows: where ε denotes the strain, and m refers to a shape parameter.In the above equation, m and a are nonnegative.
Based on the fundamental equation of continuum damage mechanics, it can be drawn that Then, Before reaching the peak value of the stress-strain curve, there are four geometric conditions (as shown in Figure 2): ε max denotes the peak-load strain value, and σ max denotes the peak-load stress value.Equation ( 8) can automatically satisfy the above geometric conditions ① and ②.The strain derivative of Equation ( 8) can be calculated using the following equation: Combining geometric condition ③ with Equation ( 8), the following equation can be drawn: Taking natural logarithms of both sides, this can be written as Combining geometric condition ④ with Equation ( 9), the following equation can be drawn: Taking natural logarithms of both sides, this can be written as Comparing Equations ( 11) and ( 12), we get Through Equations ( 12) and ( 14), the following equation can be obtained: Substituting Equation ( 15) into Equation ( 6), the following equation can be obtained: Equation ( 8) can automatically satisfy the above geometric conditions 1 and 2 .The strain derivative of Equation ( 8) can be calculated using the following equation: Combining geometric condition 3 with Equation ( 8), the following equation can be drawn: Taking natural logarithms of both sides, this can be written as Combining geometric condition 4 with Equation ( 9), the following equation can be drawn: Taking natural logarithms of both sides, this can be written as Comparing Equations ( 11) and ( 12), we get Through Equations ( 12) and ( 14), the following equation can be obtained: Substituting Equation ( 15) into Equation ( 6), the following equation can be obtained: Equation ( 16) is the equation of damage evolution.From Equations ( 14) and ( 16), it can be seen that the damage parameter D at any point in the gypsum rock mass is related to the elastic modulus, strength, and peak-load strain of the rock mass, as well as the strain on that point.Substituting the experimental data into Equation ( 14), it can be obtained as m = 1.89.To simplify the calculation, the value of m was set as 2.0.Substituting into Equation ( 16), the following equation can be obtained: Substituting Equation (17) into Equation ( 7), the constitutive equation of gypsum rock mass damage can be obtained as

The Criterion for Pillar Instability Based on Cusp Catastrophe Theory
Gypsum is mostly mined through room and pillar mining, where the pillar is usually strip-shaped.The key to this mining method lies in the stability of the pillar.To some degree, the pillar stability decides the safety and sustainability of this mining method.As the mining activities advance, the stress inside the pillar begins changing.With the internal stress of the pillar increasing, the edge of the pillar begins collapsing.When the collapse reaches a certain degree, catastrophic instability of the pillar occurs [14][15][16][17][18][19].
According to the theory put forward by A. H. Wilson, the zone (ranging from the peak stress to the edge of the pillar) is called the yield zone (or plastic zone), in that the pillar stress exceeds the yield point; the zone (ranging inward from the peak stress) is called the elastic core zone, in that the stress on the rock mass does not exceed the yield point and still complies with the elasticity rule.Setting the width of plastic zone as Y, the width of elastic core zone as W p -2Y, the strip mining depth as H, the gravity density of the overlying strata as γ, the pillar width as W p , and the mining width as W m , then the mechanical model of the strip-shaped pillar can be established (as shown in Figure 3).The actual load of the pillar can be expressed as Sustainability 2018, 10, x FOR PEER REVIEW 5 of 12 Equation ( 16) is the equation of damage evolution.From Equations ( 14) and ( 16), it can be seen that the damage parameter D at any point in the gypsum rock mass is related to the elastic modulus, strength, and peak-load strain of the rock mass, as well as the strain on that point.Substituting the experimental data into Equation ( 14), it can be obtained as m = 1.89.To simplify the calculation, the value of m was set as 2.0.Substituting into Equation ( 16), the following equation can be obtained: Substituting Equation (17) into Equation ( 7), the constitutive equation of gypsum rock mass damage can be obtained as

The Criterion for Pillar Instability Based on Cusp Catastrophe Theory
Gypsum is mostly mined through room and pillar mining, where the pillar is usually stripshaped.The key to this mining method lies in the stability of the pillar.To some degree, the pillar stability decides the safety and sustainability of this mining method.As the mining activities advance, the stress inside the pillar begins changing.With the internal stress of the pillar increasing, the edge of the pillar begins collapsing.When the collapse reaches a certain degree, catastrophic instability of the pillar occurs [14][15][16][17][18][19].
According to the theory put forward by A. H. Wilson, the zone (ranging from the peak stress to the edge of the pillar) is called the yield zone (or plastic zone), in that the pillar stress exceeds the yield point; the zone (ranging inward from the peak stress) is called the elastic core zone, in that the stress on the rock mass does not exceed the yield point and still complies with the elasticity rule.Setting the width of plastic zone as Y, the width of elastic core zone as Wp-2Y, the strip mining depth as H, the gravity density of the overlying strata as γ, the pillar width as Wp, and the mining width as Wm, then the mechanical model of the strip-shaped pillar can be established (as shown in Figure 3).The actual load of the pillar can be expressed as (1 ) The constitutive relationship curve of the pillar yield zone is different from that of the elastic core zone.The former can be expressed as follows: The constitutive relationship curve of the pillar yield zone is different from that of the elastic core zone.The former can be expressed as follows: where , ε max denotes a constant, and E denotes the initial elastic modulus.
Inside the pillar yield zone, the pillar height is h.Then, the relationship between the load P s and the deformation u can be shown using the following equation: where u 0 denotes the deformation value corresponding to the peak load.
On the contrary, in the elastic core zone, the relationship between the load P e and the deformation u can be shown using Equation (22).
According to the above analysis, the strain energy in the yield zone and the elastic potential energy in the elastic core zone can be expressed respectively as follows: The gravity potential energy of the overlying strata can be obtained with the following equation: Based on the above analysis of the mechanical model, the total potential energy function of the model system shown in Figure 3 can be expressed as follows: This equation can be used to conduct the cusp catastrophe analysis with u as a state variable.Evaluating the first derivative of V and setting it as zero, the equation of equilibrium surface V can be obtained as The cusp can be obtained by means of the smooth property of the equilibrium surface, which is shown as follows: From Equation (28), it can be drawn that u = u 1 = √ 3u 0 is the inflection point of the equilibrium surface equation.In order to establish the cusp catastrophe model, Equation ( 27) is expanded at the point u = u 1 = √ 3u 0 according to the Taylor formula, and the cubic term is taken.
Accordingly, the cubic expansion of the Taylor formula of the equation of equilibrium surface V at the point u = u 1 = √ 3u 0 can be obtained.
To simplify the calculation, a dimensionless quantity x is introduced as a state variable, and p and q are introduced as control variables.Then, the relationship can be expressed as follows: From Equations ( 31)-( 34), the equilibrium equation of the standard form of cusp catastrophe can be obtained, with x as a state variable and p and q as control variables.
Solving Equation ( 35) and inflection equation 3x 2 + p = 0, and eliminating x, the control equation can be obtained.The bifurcation is the projection of the crease of the equilibrium surface on the p-q plane.Its control equation can be expressed as follows: There are two sectors in the plane of control variables, as shown in Figure 4.In Sector E, ∆ > 0; thus, the support system remains stable.In Sector J, ∆ ≤ 0; thus, system instability occurs.Substituting Equation (38) into Equation (36), the bifurcation equation of the catastrophe for this system can be obtained, as expressed in Equation (39).Equation (38) into Equation (36), the bifurcation equation of the catastrophe for this system can be obtained, as expressed in Equation (39).
The above equation is the expression of the necessary and sufficient mechanical criterion for pillar instability in strip mining.When Δ ≤ 0, pillar instability occurs; when Δ > 0, the pillar remains stable.

Analysis of Pillar Stability
According to the ultimate strength theory [20][21][22][23], the pillar safety factor f can be obtained through pillar ultimate stress σp and the real stress in it σs [24][25][26][27].To calculate the strength of pillars, The above equation is the expression of the necessary and sufficient mechanical criterion for pillar instability in strip mining.When ∆ ≤ 0, pillar instability occurs; when ∆ > 0, the pillar remains stable.

Analysis of Pillar Stability
According to the ultimate strength theory [20][21][22][23], the pillar safety factor f can be obtained through pillar ultimate stress σ p and the real stress in it σ s [24][25][26][27].To calculate the strength of pillars, many researches were conducted and a dozen calculation methods were summarized, among which the Bieniawski formula is widely accepted and applied (as shown in Equation ( 40)).
where σ denotes the uniaxial compressive strength of gypsum rock mass in the experiment.
Then, the stress that the pillar receives can be expressed as Furthermore, the pillar safety factor f can be obtained through the following equation: From the mechanical analysis of pillar stability based on the catastrophe theory, it can be concluded that, to guarantee pillar stability, catastrophic instability must be avoided.Therefore, the bifurcation equation has to be larger than zero (∆ > 0).
Since the second item in the above equation is larger than zero, the necessary and sufficient condition for ∆ > 0 is as follows: Substituting k into the above Equation (43), the condition can be written as

Engineering Application
The Zhou Lou mining section was used in this case, with a relatively stable occurrence.It is located in the gypsum mine in Pizhou city, Jiangsu, China.With a relatively good stability, its actual parameters are as follows: obliquity = 6-10 • ; deposit thickness = 14.79-22.81m; hardness factor (f ) = 2-4.From the site investigation, it is known that the mean buried depth of the mining area is 400 m, and the unit weight is 2.5 N/m 3 .The mean compressive strength of gypsum is 30 MPa, and the mean height of the pillar is 8 m.Seismic waves travel through rocks of different properties at different speeds.Even within the same stratum, the speed of the seismic wave varies in accordance with rock strength, porosity, and density.When a roadway is excavated, the stress conditions of the surrounding rock change accordingly, which leads to the deformation of the surrounding rock and the change in mechanical properties of the rock mass.Correspondingly, the speed of the seismic wave changes as well.Using a seismic wave tester, the speed of the seismic wave at boreholes of different depths in the surrounding rock was measured, thus determining the width of the plastic zone.The development of the pillar yield zone was about 1-2 m; thus, Y was set as 1.5.Substituting these values into the above equation, the pillar width was found to be larger than 4.3 (w p > 4.3).Therefore, in order to guarantee the pillar stability, its width should be larger than 4.3 m.
The pillar stability is the key factor in mines using the room and pillar mining system; thus, its safety factor has to be larger than 1.5.In summary, there are two criteria for pillar stability.
The mining area was divided into three classifications according to the size of pillars, namely A, B, and C, as shown in Figure 5. Calculations were conducted using the above-obtained pillar stability criteria, and the results of the calculations are shown in Table 1.As is shown in Table 1, the design sizes of pillars in A and B are larger than 4.3 m and the safety factor is larger than 1.5.The theoretical calculation indicates that these two classifications are relatively stable; by contrast, the design size of pillars in C is smaller than 4.3 m, and the safety factor smaller than 1.5, which means that C satisfies the instability conditions.As discovered through field investigation, A and B show a relatively high stability, whereas C witnessed collapse.This observation verifies the practicability of the instability criteria.

Conclusions
Based on the characteristics of surrounding rock damage, a damage evolution model of gypsum rock mass was established by analyzing the damage mechanism of rock mass; meanwhile, the constitutive equation of gypsum rock mass damage was improved by introducing a shape parameter.
Based on the catastrophe theory, the pillar instability mechanism was analyzed.The criteria for pillar instability were deduced by means of the constitutive equation of gypsum rock mass damage.
According to the interrelationship between pillar load and pillar strength, the safety factor of the pillar was calculated.By analyzing the criteria for pillar instability, the necessary and sufficient conditions for pillar instability were deduced.On this basis, two rules for pillar stability were concluded: the pillar safety factor should be larger than 1.5 (f > 1.5) and the pillar width should be larger than 4.3 m (wp > 4.3).The two rules obtained herein were verified in an engineering case-the Zhou Lou mining section of the Pizhou gypsum mine.The calculation results indicate that classifications A and B are safe zones, while classification C is an unstable zone.Meanwhile, the field investigation revealed that collapses occurred in C, which in turn validates the applicability of both rules for analyzing the stability of goaf in the Pizhou gypsum mine.The derived criteria should be As is shown in Table 1, the design sizes of pillars in A and B are larger than 4.3 m and the safety factor is larger than 1.5.The theoretical calculation indicates that these two classifications are relatively stable; by contrast, the design size of pillars in C is smaller than 4.3 m, and the safety factor smaller than 1.5, which means that C satisfies the instability conditions.As discovered through field investigation, A and B show a relatively high stability, whereas C witnessed collapse.This observation verifies the practicability of the instability criteria.

Conclusions
Based on the characteristics of surrounding rock damage, a damage evolution model of gypsum rock mass was established by analyzing the damage mechanism of rock mass; meanwhile, the constitutive equation of gypsum rock mass damage was improved by introducing a shape parameter.
Based on the catastrophe theory, the pillar instability mechanism was analyzed.The criteria for pillar instability were deduced by means of the constitutive equation of gypsum rock mass damage.
According to the interrelationship between pillar load and pillar strength, the safety factor of the pillar was calculated.By analyzing the criteria for pillar instability, the necessary and sufficient conditions for pillar instability were deduced.On this basis, two rules for pillar stability were concluded: the pillar safety factor should be larger than 1.5 (f > 1.5) and the pillar width should be larger than 4.3 m (w p > 4.3).The two rules obtained herein were verified in an engineering case-the Zhou Lou mining section of the Pizhou gypsum mine.The calculation results indicate that classifications A and B are safe zones, while classification C is an unstable zone.Meanwhile, the field investigation revealed that collapses occurred in C, which in turn validates the applicability of both rules for analyzing the stability of goaf in the Pizhou gypsum mine.The derived criteria should be adjusted to the conditions in any site, using proper values of the parameters introduced in Equations ( 45) and (46).

Figure 4 .
Figure 4.The division of control space by bifurcation.

Figure 4 .
Figure 4.The division of control space by bifurcation.

Table 1 .
Geometric parameters of goaf of different classifications.