A Shaft Pillar Mining Subsidence Calculation Using Both Probability Integral Method and Numerical Simulation

In order to prolong the life cycle of the coal mine, Jinggezhuang (‘JGZ’) coal mine decided to excavate the shaft pillar. The first panel 0091 was designed near the pillar boundary as an experiment in shaft pillar mining. Both probability integral method (PIM) and FLAC were used to evaluate the influence on the shaft safety. PIM parameters were obtained from previous surface subsidence station. The rock property is based on the lab mechanical test. A simulated FLAC model containing shafts and a panel was built based on stratigraphic information. Surface subsidence results of PIM show that the 0091-panel excavation has no influence on the shafts. The simulated results show that the subsidence of the main shaft and air shaft is small and can be ignored, but it could cause the auxiliary shaft 220 mm horizontal displacement. So, the stress and displacement of the underground part of shaft were analyzed, it shows that the stress changes, subsidence and displacement are mainly located at the top part of the shafts. According to the stress and movement of the simulated shafts, 0091 was decided to be excavated and a surface monitor line was built and measured. In comparison of PIM, FLAC, and measured data, the PIM results fit the surface subsidence better. And the FLAC results have smaller maximum subsidence and greater influence area than measured. But FLAC can provide more details such as displacement, subsidence, stress and strain of both surface and underground. So, for a planned mining excavation, both methods should be used especially for the evaluation of deformation of underground constructions. In the future, with the development of the rock numerical computation technology, the numerical simulation method will be recommended first. The research shows compare of two methods of the coal mine subsidence calculation and provides a solution method for shaft pillar mining.


Introduction
After excavation of the underground coal, the goaf space makes a new unbalance of the overlaying rock. The rock will collapse, results in deformation and redistribution until a new balance attained. The deformation will spread to surface if the excavation space big enough. This process and phenomenon are called 'mining subsidence' [Zhao and Meng (2010)]. Mining subsidence can cause the surface constructions such as house, railway, water body, etc. to be destroyed and also create environment problems. People have recognized this problem very early; the history of mining subsidence research can be traced back to the 15 th century. In the first years of the 20 th century, a lot of surface deformation surveys were done by coal mine companies, and mining subsidence became a new area of research. Then, with the massive energy demand of industry, the research of mining subsidence increased, and a lot of mining subsidence prediction theories and methods were proposed based on both continuum and discontinuous mechanics theory [Cui and Deng (2017)]. Those methods played important roles in reducing the damage caused by coal mine subsidence. In China, a stochastic medium theory from Poland was introduced and used in coal mine subsidence calculation in the 1960s [Litwiniszyn (1957[Litwiniszyn ( ,1974; Liu and Liao (1965)]. The method was improved into an easier model named probability integral method ('PIM') [Liu and Liao (1965); He, Yang, Ling et al. (1991); Cui and Deng (2017)]. PIM is an influence function to calculate the surface subsidence based on the geometry integration of excavation space. There are 8 parameters of PIM that can be obtained from back-analysis of the site measured data. It is the most used method in China [Coal Industry Bureau of People's Republic of China (2000)]. But PIM is based on the influence function, it is a semi-empirical function which depends on the parameters what used but not on the deformation mechanism of the rock mass. So, it can calculate only the ground surface displacement but cannot be used to calculate deformation and stress distribution of under lying rock. With development of the rock simulation theories and computer performance, numerical simulation technology become an important way to deal with coal mine engineering problems [Costa, Zingano and Koppe (2000); Nuric, Nuric, Kricak et al. (2000); Pereira, Costa, Salvadoretti et al. (2012)]. Simulation of mining subsidence also have been used widely since the 1980s [Chrzanowski, Monahan, Roulston et al. (1997); Szostak and Chrzanowski (1991)]. Fast Lagrangian Analysis of Continua in three dimensions (FLAC 3D ) is a numerical modeling and simulation software for rock, groundwater, constructions, coal mine and ground support [Bock (2015); Zhao, Jiang and Lan (2015); Shibata, Zarlin, Shimada et al. (2014)]. FLAC 3D can define both linear and nonlinear constitutive models for the simulation elements. The models can calculate deformation and displacement with large deformation mode. The software uses the explicit Lagrange difference to deal with 3D, large area plastic and flowing deformation problems with small computer memories. JGZ coal mine is located in the east of China; the coal resource is in a basin-like syncline. It is about 3.5 Km from North to South and 3.4 Km East to West, the area is 9.23 Km 2 . The coal mine is in the alluvial valley plain area and has great traffic network. The elevation is +38.9 m in the north and +23.9 m in the south, average +35.0 m. The coal mine was first built in 1958 and began production in 1978. It has 140 million tons of coal resource of which 72.68 million tons are recoverable reserves. The designed production is 1.2 million tons per year. The coal mine development method is mainly multi-level based on vertical shaft. There are three shafts located in the industrial site. With about 30 years excavation, JGZ coal mine faced the problem of resource exhausting. The coal mine manager had a plan for mining the industrial site pillar. In order to maintain the shaft safety, an experimental panel was designed first. This paper predicted the mining subsidence of the experimental panel with both PIM and FLAC 3D , to evaluate excavation influence on the shaft, and then results comparison with the data measured was done in the paper.

Introduction of PIM model
PIM is an influence function of the mining surface subsidence prediction model based on stochastic medium theory [Litwiniszyn (1957); Liu and Liao (1965)]. In order to calculate the surface subsidence, the coal mine influenced rock can be hypothesized as stochastic, discontinuous elements as shown in Fig. 1(a).  [He, Yang, Ling et al. (1991)] The figure shows that if one element of the first layer is excavated, one of two elements in the second layer will drop into the excavated space, and it is random which of the two elements will drop in. With this hypothesis, suppose that the subsidence will spread to the surface after one element of the first layer excavation; the surface subsidence curvature and distribution are similar with the Gauss bell function as shown in Fig. 1(b) [Yang, Liu and Wang (2004); Zhang (1998)]. For the large space excavation, the subsidence can be calculated by the geometry integral. As shown in Fig. 2, surface point A (x, y)'s movement can be decomposed into subsidence W (x, y) and horizontal displacement U (x, y, φ). Dj is the area of j th panel which have influence to point A. According to the theory of PIM, the surface subsidence W (x, y) and horizontal displacement U (x, y, φ) can be calculated by Eq. (1)-Eq. (2) [He, Ling, Yang et al. (1991)]. In China, PIM is the most used function for coal mine subsidence prediction and plays an important role in reduction the loss of mining subsidence.
where (x, y) are coordinate of the surface point, W, U are mining subsidence and horizontal displacement, φ is direction angle of surface movement; W0 is the maximum ground subsidence, given by W0=mqcosα; We(x, y) is subsidence of a small unit mining; H is mining depth; m is mining thickness; D is calculation mining area of the panel, n is total number of panels, j is the number of panel; α is dip angle of coal seam; q is subsidence factor; b is displacement factor; tanβ is tangent of major influence angle; r is major influence radius given by r = H/tanβ; dsdt is integration variable of double integral of area Dj.
The PIM function results depend on those parameters. The parameters can be obtained from back analysis of measured data [Li, Peng, Tan et al. (2017)] or by the experience function of overlaying rock mechanical parameters [Coal Industry Bureau of People's Republic of China (2000)].

Introduction of FLAC 3D
FLAC 3D is a simulation software developed by ITASCA company. It is a 3D finite difference solving software. It can be used for rock mechanical simulation and plastic flow in soil, rock and other material. It has a beautiful performance on large deformation [Xie, Zhou, Wang et al. (1999)]. Compared with other simulation systems like ANASYS, ADINA, FLAC 3D can solve the physical instability problems. It can be used in stability analysis of highly nonlinear geometry and physical problems like coal mine rock deformation, tunnel crack, and stress field of surrounding rock [Li, Tan, Gu et al. (2015); Rawamurthy (2014)].
FLAC 3D uses the dynamic explicit method to get the elements' time step solution, then the model's failure and crack can be traced. It is important to research time and space effect of the rock deformation. FLAC 3D has 11 different material constitutive models like null model, elastic model, and mohr-coulomb model, etc. Those materials can be used to simulate the complex geological and mining environment. The software also has powerful pre-and postprocessing functions. Complicated models can be built with its basic grid and a lot of results like stress, strain, and displacement can be print out with values or plots.
There are 6 basic steps of simulation with FLAC 3D as shown in Fig. 3.

Figure 3: Flow chart of simulation with FLAC3D
Step 1. First of all, finite difference grids similar to the initial geological should be built.
The software provides 13 different grids which can simulate arbitrarily topography.
Step 2. Then, the grid should be given a constitutive model and material parameters. And the model also needs a boundary and initial condition.
Step 3. Next step is to get the initial rock stress state.
Step 4. The model excavated or changed according to the actual engineer conditions.
Step 5. Calculation to get a new balance for the engineering changes.
Step 6. The results can be printed out and analyzed after the rock is re-balanced through calculation.

Mining subsidence calculation and its parameters 3.1 Basic geological of 0091 panel
There are four coal seams in the research area; the 9 # coal seam is which investigated.  The strata of this area from bottom to top are Carboniferous upper series (C3), Permian lower series (P1), and Quaternary loose deposit (Q). The thickness of the Quaternary loose layer is about 160 m. The geology of investigated panel is simple. The geometry and geology of section 1-1 are shown in Fig. 5. The direct roof of the 9 # coal is mud stone, shale, sandstone, etc. And direct floor of 9 # coal is clay stone. The 0091 panel is the first panel in industrial site pillar; an observation station was set in order to get the displacement law of the ground and protect the constructions in the industrial site.

Parameters used to predict mining subsidence based on PIM
To use the PIM method, there are 8 parameters that need to be determined first. Those parameters are often obtained based on previous measured mining subsidence data [Zha, Feng and Zhu (2001)]. The JGZ coal mine and other collieries in Tangshan area have built a lot of sites monitoring stations to study the ground movement, and PIM parameters are shown in Tab. 2.

Rock mechanical parameters of FLAC 3D model 3.3.1 Boundary condition and initial stress
A reliable boundary condition is the basic requirement to get correct results of the mechanical simulation. In this shaft pillar mining simulation, the model boundary is displacement boundary. It meets the following conditions: 1. Fix the horizontal displacement of the periphery. As shown in Fig. 6, fix the Xdirection displacement in right and left side plane, and fix the Y-direction displacement in front and back side plane. 2. Fix the vertical displacement of the bottom; 3. Make the top surface of the model move free. The model is simplified for ignoring tectonic stress, it just in the hydrostatic state. Gravity stress field is the only thing to be considered. The initial stress of model depends on the weight and property of the rock. The model's boundary conditions are shown in Fig. 6. In addition, to reduce the error of the boundary effect, distances between the model geometry boundary and the edge of the ground subsidence basin are more than 200 m. The model geometry is big enough to simulation for 0091 panel excavation deformation.

The geometry
The simulation geometry of the model depends on the actual geological conditions and its boundary conditions. The 3D model was established by layers according to exploration results, and the 0091-excavation area is represented as the null elements, as shown in Fig.  7. The model elements were built with wedge-shaped mesh which is a built-in grid of FLAC 3D . The null model was used for 0091panel excavation.   (2000)]. As shown in Fig. 8(a), the maximum subsidence located above the center of 0091 panel and the value reduced from center to outside. The subsidence concentrated in a certain distance from panel center. The horizontal displacement is directing to the panel center. The biggest horizontal displacement is located above the panel boundary. The horizontal displacement is smaller in panel center and basin boundary. The horizontal displacement influence area is almost equal to subsidence influence area. In Fig. 8(a), Shafts are not influenced by the 0091coal mine subsidence. It is safety to excavate the pillar panel 0091. In Fig. 8(b), the maximum subsidence located above the center of the panel. The horizontal displacement is directing to the panel center and the biggest horizontal displacement are located above the panel boundary. The horizontal displacement influence area is bigger than the subsidence influence area. As shown in the Fig. 8(b), main shaft and auxiliary shaft are located in the mining subsidence basin area. It shows that those shafts are influenced by the 0091 excavations. It needs more analysis to decide if it is safe to excavate the panel 0091. More information can be exported from the FLAC 3D . By comparing two figures, more details can be found. First, the basic appearance of subsidence and horizontal displacement are similar. But the FLAC 3D subsidence basin area is bigger than PIM results. The ratio between subsidence and horizontal displacement is different; FLAC 3D 's horizontal displacement influence area is bigger than PIM results. The influence results on shafts are different, PIM results have no influence to shafts and FLAC 3D have two shafts influenced. The differences of basic principles make those results difference, more information and analysis are needed for decision.

Shaft displacement and stress
As the FLAC 3D can simulate the mechanical behavior of rock mass, more information for engineering evaluation such as stress can be exported from the model. In order to maintain the shaft safety, the subsidence (z-axis displacement), x-axis displacement, y-axis displacement, stress-ZZ(σzz), stress-YY(σyy), stress-XX(σxx) are exported and shown in Fig. 9. As shown in Tab. 4, the maximum subsidence of shafts is on the top of the auxiliary shaft, the subsidence is 18 mm. The subsidence of the main shaft and air shaft is less than 10 mm. The horizontal displacement of the auxiliary shaft is 187 mm in x-axis and 137 mm in yaxis which is the maximum horizontal displacement of three shafts. But as shown in Fig.  9, the displacement and stress of shafts are all located at the top part which near the surface. The Fig. 9(d) shows that it has little stress concentration on shaft wall. The stress is mainly from gravity effect, but not from excavation of 0091 panel. So, synthetically with subsidence results and stress distribution of the shafts which are obtained from both PIM and FLAC 3D , the management decided to excavate the panel 0091 and measured the surface subsidence to keep the shaft from being damaged. FLAC 3D also can provide more information of internal rock mass, as shown in Fig. 10. They are contours of z-axis displacement, y-axis displacement, x-axis displacement and stress-zz (σzz) of Section 1-1. The deformation of the inner rock is symmetrical of the panel center. The subsidence decreases progressively from the roof of the panel to the surface. The floor of the panel's z-axis displacement is positive which indicates that the floor of panel is rising. It simulated the floor unloading phenomenon after the panel excavation. Fig. 10(d) shows that there are two stress concentration areas above the panel which is located at the panel boundary. It is similar to the actual excavation.

Results comparison and analysis 4.3.1 Results comparison
With the analysis above, the company decided to excavate the 0091 panel. And subsidence measurement of B line in Fig. 4 was surveyed from September 26, 2008 to April 30, 2010. And subsidence of those points was obtained. In order to analyses the results, the PIM results, FLAC 3D simulated results and measured data are plotted in Fig.  11. There are some similarities and differences.

Figure 11: Results comparison
As shown in Fig. 11, all results of both subsidence and horizontal displacement are symmetry over the center of the 0091 panel. The maximum subsidence is located at the center of the panel. The horizontal displacement direction is towards to the center of the goaf and the maximum displacement is located above the boundary of the panel. Both FLAC 3D and PIM results are fit for the basic law of the mining subsidence. But there are some differences, as shown in Fig. 11(a); the subsidence boundary of the PIM calculated results is smaller than measured but bigger than the FLAC 3D results. The FLAC 3D results have larger range of influence than measured. As shown in Fig. 11(a), the angles of draw of three results are 39 o , 47 o , 57 o . The maximum subsidence result is WPIM>WMeasured>WFLAC 3D . The PIM results fit better with measured data than FLAC 3D . The horizontal displacement was not measured; PIM and FLAC 3D horizontal results are compared in Fig. 11

Results analysis
As shown in Fig. 8 and Fig. 11, there are significant differences such as influence area, value and position of maximum subsidence, horizontal displacement value. The reason of differences is originated from the basic theory of two methods. Those two methods are come from two different sects of rock mechanics. PIM is based on discontinuous particle medium mechanics theory. The theory holds that rock is constituted by sand or small rock particle. The particles do not have any relationship with each other and can have relative movement. It assumes that the movement of a great quantity of granular medium is a stochastic process. For the center of the excavation panel, it has little stress between rock bedding. PIM theory coincides with the stochastic medium hypothesis. But at the position of panel boundary, actually the rock is subjected the force of tensile stress to the center of the panel. With this force of tensile stress, the rock above or outside the panel boundary could have horizontal displacement and subsidence movement toward the panel center. For this reason, the calculation deformation results outside the panel are smaller than measured data and the subsidence basin area is also smaller than measured area. FLAC 3D is a fast Lagrangian analysis of continua software. The software solves rock mechanical behavior with explicit finite difference. To make a good simulation model, it needs complex geological exploration and mechanical tests. FLAC 3D is good at solving the problem of continuous medium. But rock body deformation is an extremely complicated problem. Rock body contains different material, soft surface, fault and cracks. To simplify the coal mine model to be a layered continuous rock mass need to ignore the discontinuous property of the actual rock. It will reduce the maximum subsidence value and enlarge the subsidence basin area. As shown in Fig. 8 and Fig. 11, both results of PIM and FLAC 3D have its own superiority and defects. PIM calculate the subsidence depends on the actual measured parameters, the deformation results of ground surface is better. But FLAC 3D can calculate the dynamic change of the rock, it also can provide more information on time effect and space effect of mining rock deformation process, it is very important to safety evaluation, predict and analysis [Xie, Zhou, Wang et al. (1999)].

Discussion
There are four main categories of mining subsidence prediction method based on which theories and functions are used. They are phenomenological method, rock mechanical method, physical simulation method and numerical simulation method [Cui and Deng (2017)]. PIM is one of the influence functions in phenomenological methods. PIM is mostly used in China and it has produced huge economic benefits to reduce the coal mine subsidence damage. There are a lot of excellent reasons to use PIM. 1. PIM is based on stochastic medium theory which is experimentally verified [He, Yang, Ling et al. (1991)]. 2. PIM just needs 8 parameters which can be obtained from measured data easily. The parameters such as subsidence factor are very easy to be understood by engineers in the coal mine.
should be used to evaluate the deformation influence of coal excavation especially for underground facilities. And we recommend that more simulation method should be used with the development of dynamic rock model and numerical simulation technology in the future.

Conclusion
(1) In order to excavate the shaft pillar, PIM and FLAC 3D simulation were used to calculate the surface subsidence and influence on the shafts. The PIM results show that the panel 0091 excavation has no influence on the shafts. And the simulated results of FLAC 3D shows that the subsidence of main shaft is 9 mm and auxiliary shaft is 18 mm. FLAC 3D results have bigger influence area than PIM calculated. With the analysis of shaft and rock displacement and stress, it shows that the shaft can keep safe after 0091 panel excavated.
(2) Comparison shows that the PIM fit better to the measured surface subsidence than FLAC 3D . Mining subsidence influence area of PIM results are smaller than measured but the FLAC 3D results are bigger. PIM is easier to be used and can predict the surface subsidence with 8 parameters; also, the surface deformation results are better than FLAC 3D . But FLAC 3D can provide more information like strain, stress, etc., of the surface and underground rock. It is useful for underground facilities evaluation. So, both PIM and numerical simulation method should be used in a coal mine subsidence evaluation especially for a project of underground deformation calculation and analysis.
(3) With development of rock mechanical theories, more simulation method should be used in the future