Abstract
The stability of the goaf support system is the key to safe production in gypsum mines. Therefore, this study constructed a pillar-beam support system which contained pillar plastic zones. In this support system, the beam and pillar were taken as energy releaser and energy dissipater, respectively. Through establishing a cusp catastrophe model based on energy theory, the new criterion for instability was obtained which is related with geometric stiffness and system energy dissipation. The results indicate the instability of the support system is caused by the incompatibility of energy release, dissipation, and geometric deformation. When K > 1, the energy released by the support system is compatible with geometric deformation. The support system experiences a quasistatic process from the static state in bottom page to the static state in top page along Path I. When K < 1, the energy released by the support system cannot be in tune with geometric deformation. The support system experiences a catastrophe process along Path II. The evolution from the static state in bottom page to the static state in top page is not progressive, but catastrophic. The redundant energy released in this process leads to mechanical instability of the support system. This study provided theoretical foundation for the mining and treatment of mines. Based on actual engineering examples, the sensitivity of the geometric parameters of the support system was analyzed as well. These parameters are ranked by their sensitivity from high to low, as is shown below: beam thickness, plastic zone width, room span, pillar width, and pillar height. Then, the goaf was classified according to the geometric parameters. Energy catastrophe theory was applied to analyze the stability of the support system in different classes of goaf. The analysis results showed that Class D goaf should be labeled as the unstable zone, which was consistent with the result of field research. To conclude, energy catastrophe theory can be used to demonstrate the nonlinear mechanical mechanism of support system instability in room-pillar mining goaf.
1. Introduction
Gypsum and anhydrite are sulfate minerals which are generated in the widely existing surface water. The compressional movement of Eurasian plate and Pacific Plate led to the expansion of land area and the shrink of the Pacific Ocean. During this process, many marine faces, lacustrine facies, and closed or semiclosed salt basins came into being, which in turn were generated into sulfate ore deposit by evaporation and deposition. Therefore, the reserves of gypsum are quite abundant in eastern coastal China.
The gypsum mine is mainly exploited with room and pillar mining. As the consequence of constant mining, a large amount of goaf appears, which poses potential threats to the safety of mining production and the life of residents nearby. Therefore, it is of significance to conduct theoretical researches on the failure mechanism of the pillar-roof support system of gypsum mines.
In the course of gypsum mining, it is a common practice to reserve a roof of certain thickness which integrates with the rock pillars to form an overall structure. This structure is different from the overlapping form of roof and pillar in coal mines, thus forming the unique “pillar-roof” support system of gypsum mines. However, due to the lack of standardization of mining design and production, a large quantity of goaf is either too high or too wide; the roof is in some cases too thin, the room span is too large, and the pillars are too narrow. All these constitute an unstable support structure. Great efforts have been made in the field of goaf stability by scholars both at home and abroad [1–5]. Zhao et al. [6] employed catastrophe theory and safety-factor strength reduction method to evaluate the stability of goaf roof and studied the safety stock of overlapping roof. Henley [7] used energy theory to analyze the dynamic process of system instability of gypsum mines. Pan and Wang [8] applied catastrophe theory to analyze the factors that lead to the structural instability of the narrow pillar support system. Xu et al. [9] explored the rock-burst mechanism of the coal pillar under a hard roof from the perspective of cusp catastrophe theory and summed up the occurrence criteria of rock burst. Upon this basis, they also discussed the factors that influence rock burst as well as their respective degree of influence. Qin and Wang [10] established a physical model to predict the instability evolution of the pillar-hard roof mechanical system by means of catastrophe theory. They also put forward the mechanical criteria for the necessary and sufficient condition of instability and the formula for deformation jump of instability. They applied catastrophe theory to the research on roadway instability, water-inrush from floor, plate rock mass instability, and other problems about stope instability [11–14]. Zheng [15] explored nonlinear engineering problems such as the instability of high dam and tunnel surrounding rock from the perspective of catastrophe theory. Other methods related to rock damage were also used to study the goaf pillar, such as numerical methods, FEM, X-FEM, ED-FEM, and phase field [16–18]. These abovementioned methods, however, cannot reflect the pillar damage from the nature of gypsum rock damage. Therefore, the constitutive law of gypsum damage is proposed in this paper to analyze the stability of goaf in combination with catastrophic theory.
When the goaf of gypsum comes into being, under the pressure of overlying strata, the roof begins to accumulate resilience and transmits it toward the pillar. Once the resilience accumulated within the pillar-roof support system reaches its critical point, it will be released instantaneously, causing pillar instability and roof collapse. This occurrence, due to its abruptness, can therefore be analyzed by catastrophe theory.
Previous studies also considered the influence of the pillar plastic zone, but the mechanical models established, more often than not, took coal pillars as the research object while ignoring the cooperative deformation of the coal pillar and the roof. Even though some catastrophe models of pillar-roof cooperative deformation were proposed, the influence of the plastic zone and its geometric parameters were excluded. But in fact, the size of plastic zones and the geometric parameters play a critical role in the stability of the support system. Therefore, when creating the catastrophe model of pillar-roof cooperative deformation support system, it is closer to the actual engineering situation to integrate the geometric parameters of mining and the plastic zone of surrounding rock.
2. The Mechanical Model of Support System and Energy Dissipation Characteristics
Since the room of room-pillar mining is relatively symmetric, the mining floor is assumed not to deform for simplicity. The roof of the long working face can be taken as the rock beam. The self-weight of the beam and pressure of overlying strata can be simplified as a uniform load whose intensity is q. The middle part of the roof which is located between two pillars is taken as a structural unit for analysis (as is marked off with the blue dashed line in Figure 1). Among them, h is the height of the ore pillar, a is the width of the pillar, and b is the span of the mine). Under the pressure of the overlying load, the roof and the pillar deform cooperatively (the red part represents the deformation u). The rock beam remains elastic throughout the whole process of pillar deformation and failure.
According to the theory of A. H. Wilson, due to the loading of overlying strata, a yield zone (plastic zone) comes into being along the pillars. The inward rock mass stress in the yield zone does not exceed the yield point, which complies with the yielding rules. Surrounded and confined by the yield zone, this area is in a three-dimensional stress state, thus being called the elastic core zone (as shown in Figure 2), Among them, a is the width of the ore pillar, and Y is the area of the unilateral plastic zone, σ1 is the stress at the junction of the elastic-plastic zone.
The constitutive relation curve in the elastic core zone is different from that in the yield zone. The curve appears to be linear in the elastic core zone, while it appears to be nonlinear in the yield zone; similarly, it shows strain-hardening behavior in the former, while strain-softening behavior in the latter. Once reaching the peak stress, the pillar which appears to be strain softening tends to unload soon. With reference to the studies conducted by Wang et al. [19], the relation of yield zone stress, strain, and damage parameters can be expressed as follows:where denotes the corresponding strain of peak stress and denotes the initial elastic modulus. Accordingly, the force that the pillar unit exerts upon the roof in the course of deformation (the resultant force of the elastic core zone and the yield zone) or the relation of pillar load-deformation can be expressed in the following equation:where denotes the deformation value corresponding to the peak load.
If there were not the support of the pillars, the rock beam would reach a static equilibrium under the pressure of self-weight and overlying load (as shown in Figure 2). However, the supportive force of the pillar exerts restriction upon the displacement of every point on the skew curve of the rock beam. Under the action of the uniform load q, the rock beam accumulates elastic deformation energy. With the accumulation of energy, the rock beam deforms gradually. The points on the skew curve of the rock beam tend to change toward the static equilibrium position (as shown in Figure 3). In this process, the rock beam releases energy constantly while the pillar absorbs the energy released constantly. If the geometric size of the mining area is appropriate, the strain-softening deformation of the pillar is in tune with the elastic resilience of the rock beam. In other words, the energy release is in tune with the structural deformation. In this case, the support system experiences a progressive deformation. If the geometric size of the mining area is inappropriate, the energy release is not going to be in tune with the energy consumed in the process of deformation. Once the energy released exceeds the energy needed for pillar deformation, the excessive energy will turn into system kinetic energy and be released instantaneously, thus causing dynamic instability of the whole support system.
3. Energy Function of Support System and Catastrophe Model Construction
3.1. Energy Function of Support System
According to the stress condition of the support system (as is marked off with the blue dashed line in Figure 1), the total energy of the system V is composed of three parts: the elastic deformation energy released by roof flexural deflection , the energy dissipated by microcrack growth and connection within the pillar , and the external potential energy produced by the displacement of overlying strata , which is shown as follows:where denotes the width of the pillar and indicates the span of the room; denotes the external force from the weight of overlying strata and rock beam; refers to the internal force generated in the process of pillar deformation; stands for the average amount of compression; denotes the roof deflection function; stands for the elastic modulus of gypsum; and denotes the bending inertia moment of the roof.
According to the boundary conditions of rock beam deformation, the roof flexural function of the support system of the selected unit can be shown as follows:where is the positive number and denotes the average amount of compression.
Due to the coordinate deformation of the rock beam and the pillar, according to the deformation condition,
It can obtained as follows:
Let , then
Substitute Equations (4)∼(7) into Equation (3), and the total potential function of the support system can be obtained through the following equation:
3.2. Cusp Catastrophe Model of Support System Instability
Set pillar compression displacement as the state variable; then according to cusp catastrophe theory, the equation of the support system catastrophe model for equilibrium surface can be obtained (as shown in the following equation):
Apparently, Equation (9) is the equilibrium equation to analyze the stability of the pillar-roof support system in gypsum mines. According to the smoothing property of equilibrium surface, at the cusp, , then
From Equation (10), it can be obtained that at the cusp. In order to construct the cusp catastrophe model, the equilibrium surface equation is calculated by the Taylor series expansions at the cusp. Letand expand the cubic term to triple terms using the Taylor series expansion. After simplification, this can be described as follows:
Substitute dimensionless quantity as the state variable, and as the control variables, and then let
From Equations (11∼13), the standard form of the equilibrium surface equation of cusp catastrophe can be obtained. As is shown in Equation (16), denotes the state variable and and refer to the control variables [20]:
When the system is in a critical state,
The bifurcation set equation on the control plane can be shown as follows:
The cusp catastrophe model and bifurcation of the support system are shown in Figure 4. The bifurcation set has a cusp at the coordinate (0, 0). The equilibrium surface can be divided into three pages: the top page, the middle page, and the bottom page. The bottom page stands for the prenatal stage of instability as the elastic potential energy increases; the middle page stands for the unstable state of catastrophe; the top page denotes the new stable state after instability. The polygonal line divides the control variables into two domains (I and II). When the control variables change along Path I, both control variable and state variable change progressively so does the system instability. When the control variables change along Path II, even a slight change of control variables would lead to a sudden sharp increase in the state variable at the edge of the bifurcation set. Catastrophe occurs when the control variables move from the bottom page up to the top page and leaps over the bifurcation set. In this case, the support system loses its stability.
4. Analysis of Support System Instability Mechanism
4.1. Necessary Condition for System Instability
From the control plane in Figure 4, it can be seen that the equilibrium point is likely to span the bifurcation set only when α < 0. Similarly, it is only when α < 0 that structural instability of support system occurs. Therefore, the necessary condition for instability isAdd a parameter K, and is set to denote the geometric stiffness of the support system. If and the values of and are correlated only with the geometric parameters of the support system, with the lithology of the mining area being constant, the necessary condition of system instability is related with beam thickness, room span, pillar height, pillar width, and size of elastic zones along the pillar. As can be concluded according to Equation (19), K < 1 is the necessary condition for system instability.
4.2. Sufficient Condition for System Instability
According to cusp catastrophe theory, when the system state satisfies the equation of the bifurcation set, the system is the most sensitive to external disturbance, thus being in a critical state. Therefore, the equation of the bifurcation set can be taken as the sufficient condition for system instability.
When Δ = 0, Equation (18) has three real roots, one of which is in an unstable state while the rest two in a stable state. Only when the stable state leaps from one bifurcation to another, will the instability of support system occur. The equation of the bifurcation set can equal 0 only when α < 0 or K < 1, which indicates that the sufficient condition and necessary condition of system instability share highly uniformity. When K > 1, the deformation compatibility of geometric parameters is relatively high, and the geometric stiffness is constant. In this case, the accumulating energy drives the support system to change along Path I in Figure 4; both control variable and state variable change progressively. Consequently, the system instability is a gradual process. When K < 1, the deformation compatibility of geometric parameters appears to be poor and the energy released exceeds the energy consumed. In this case, the support system changes along Path II, which ultimately leads to catastrophic instability.
5. Verification of an Engineering Case
5.1. General Geology
The mining section in this case, with a relatively stable occurrence, is located in the gypsum mine in Pizhou city, Jiangsu province. Its actual parameters are as follows: Obliquity: 5∼10° Average mining depth: 480 m Deposit thickness:13.15∼18.14 m Average thickness: 14.5 m Average mining height: 8.5 m Average thickness of the roof and the floor: 3 m Room span: 7∼9 m Pillar width: 6∼8 m Hardness factor (f): 2∼4 Average elastic modulus: 4.0 GPa Poisson’s ratio: 0.3
The lithology of the roof and the floor is mainly dark grey gyp-rock, dotted with small amount of gypsum mudstone appearing to be dark purple or brown.
5.2. Sensitivity Analysis of Factors Affecting Stability
Since the lithology of the mining area is unchangeable, the goaf stability is mainly affected by the geometric parameters and the range of the plastic zones. Xia analyzed the sensitivity of the geometric parameters of goaf in gypsum mines [21], but they failed to take the influence of pillar plastic zones into consideration. At present, there are two methods of sensitivity analysis: local sensitivity analysis and global sensitivity analysis. The former is more applicable to simple relation models, while the latter to complex nonlinear models. Since the input and the output of the K-value model constructed in this paper are close to linear, local sensitivity analysis is used for simplicity. To analyze the sensitivity of plastic zones, the correlations between the geometric stiffness K and pillar width a, room span b, pillar height h, beam thickness B, and plastic zone width Y are established and shown in Figures 5∼9.
The sensitivity coefficient can be obtained through the ratio of the relative change rate of influencing factors and that of stability correlation coefficient K. Thus, the sensitivity coefficient of the ith influencing factor can be shown as follows:
According to the sensitivity analysis, the fitting sensitivity of pillar width, room span, pillar height, beam thickness, and plastic zone width is 0.29, 0.33, 0.11, 2.0, and 0.59, respectively. Therefore, the sensitivity order of influencing factors from highest to lowest is as follows: beam thickness, plastic zone width, room span, pillar width, and pillar height (as shown in Figure 10). The sensitivity of plastic zone width is just next to that of beam thickness, which makes it an important factor affecting the goaf stability. Zhou et al. [22] also conducted researches on the influence of the beam, pillar width, and height on the goaf stability. Their analysis results were consistent with this paper, which could be taken as corroborative evidence to support the reliability and feasibility of catastrophe theory and the analysis on this basis. The analysis results provide the basis for the mining design of the later mines and the treatment of the goaf. In the design, it is necessary to ensure the thickness of the reserved rock beam and reduce the span of the room. In the treatment of the goaf, the support of the plastic pillar should be strengthened to reduce the range of the plastic zone, thus ensuring the safety of the project site.
5.3. Testing of Pillar Plastic Zone
A surrounding rock loosing-circle tester was used to test the range of the pillar plastic zone (the size of Y) (as shown in Figure 11). The propagation velocity of the seismic wave differs when it travels through different rocks. Even within the same rock stratum, the seismic wave travels at different speeds due to the different rock strength, porosity, and density. By measuring the propagation velocity of the seismic wave at the boreholes of different depths in the surrounding rock, the range of the plastic zone can be determined. Boreholes were drilled to 3- to 5-meter deep at the sensing points. Then, sensor probes were placed inside the boreholes. The intelligent instrument recorded the time t when the direct wave arrived. The depth that the sensor probe was placed in the borehole was taken as the propagation distance OS. According to the data measured in the boreholes, the V-t and H-t curve chart (time-depth chart) could be drawn (as shown in Figure 12), thus obtaining the different straight slopes (the propagation velocities in different zones). The range of loosing circle could be marked out in accordance with the changes of wave velocity, as is shown in Figure 13.
It can be observed from Figure 13 that a falloff area appears within the distance of 0–1.5 m, indicating that the surrounding rock in this area is loose. The existence of fracture development and relatively small stress proves that this area could be labeled as the pillar plastic zone. As the distance increases, the wave energy appears to be relatively smooth, with no protuberance and areas of stress concentration. Thus, it can be concluded that the inner part is the pillar elastic core zone. Based on a great amount of field tests, the pillar plastic zone is distributed along the pillar within the distance of 1.5∼2 m.
5.4. Classification of Mining Areas and Stability Calculation
Due to the nonstandardized mining, the geometric parameters of mining areas are different. Based on the careful field study, the goaf was divided into four classes (as shown in Figure 14). The range of pillar plastic zones was measured with instruments. The relevant parameters are shown in Table 1. Then, stability analysis was conducted by means of Equation (19) which denoted the sufficient and necessary conditions of instability based on catastrophe theory.
Class A goaf can be taken as an example. In Class A goaf, the pillar width a = 7 m, pillar height h = 8.5 m, the range of the plastic zone along the two sides Y = 1.9 m, room span b = 8 m, and beam thickness B = 3 m. The elastic modulus in this area is 4.0 GPa and Poisson’s ratio is 0.3. Through calculation, it can be obtained that k = −1.02 and the flexural stiffness of the roof I = 2.25. Substitute the above parameters into equation of K, and then, the result is shown in following equation:
Through calculation, it can be known that the geometric stiffness (K) in Classes A and B is larger than 1, indicating that these two areas can be labeled as safety zones; the geometric stiffness (K) in Class C is equal to 1, indicating that it is in a critical state of instability. In other words, it should be labeled as the danger zone and needs isolation processing. The geometric stiffness (K) in Class D is smaller than 1, indicating that this area is the unstable zone. Moreover, as was recorded, collapses occurred for many times in this area in 2012. This indicates that the analysis results of catastrophe theory are consistent with the actual situation in the field, verifying the reliability of the model.
6. Conclusions
This paper, taking the influence of the surrounding rock plastic zone into consideration, constructs a model based on cusp catastrophe theory and analyzes the stability of the support system in the room-pillar gypsum goaf. The following conclusions are drawn:(1)Analysis on room-pillar mining layout was conducted so as to construct the mechanical structure model of the support system in gypsum mines with plastic zones. Based on energy theory, the energy release of the support system and the energy dissipation characteristics were analyzed. Meanwhile, qualitative analysis was conducted by means of geometric relevance to support system deformation. On this basis, it can be concluded that instability of the support system can be attributed to the fact that the energy released exceeds the energy consumed in the process of deformation. It is the redundant energy that leads to mechanical catastrophic instability of the support system.(2)Based on the cusp catastrophe model constructed, the sufficient and necessary conditions of system instability were analyzed. On this basis, the criterion for system instability can be concluded: geometric stiffness K − 1. When the geometric stiffness K is larger than 1, the energy released by the support system is compatible with geometric deformation. The support system experiences a quasistatic process from the static state in the bottom page to the static state in top page along Path I. When the geometric stiffness K is smaller than 1, the energy released by the support system cannot be in tune with geometric deformation. The support system experiences a catastrophe process along Path II. The evolution from the static state in the bottom page to the static state in the top page is not progressive, but catastrophic. The redundant energy released in this process leads to mechanical instability of the support system.(3)Besides theoretical analysis and calculation, the factors affecting the stability of the support system were analyzed. Accordingly, it is concluded that the principal geometric influencing factors are pillar width, room span, pillar height, beam thickness, and pillar plastic zone width. The sensitivity of these influencing factors was analyzed as well. These factors are ranked by their sensitivity from high to low, as is shown below: beam thickness, plastic zone width, room span, pillar width, and pillar height. Meanwhile, a field test on the pillar plastic zone was conducted by means of surrounding rock measuring instruments. According to the range of the plastic zone as well as different geometric parameters, the goaf was divided into four classes. Then, the abovementioned criterion of K − 1 was used to analyze the stability of each class. The results indicate that Classes A and B can be labeled as safety zones, Class C is in a critical state, and Class D should be labeled as an unstable zone. Through field research, it is found that collapses occurred for many times in Class D, which also verifies the reliability of the stability analysis on room-pillar mining goaf based upon the model constructed. Note: the derived criteria should be adjusted to the conditions in any site by using proper values of the parameters used in the formulae (19).
Data Availability
The chart data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Key Research and Development Program of China (2018YFC060470), project supported by the National Natural Science Foundation of China (51874289 and 21506248), and project supported by the Fundamental Research Funds for the Central Universities (2018ZDPY05).