The stress change law of a collapse column and the failure depth of a coal seam floor before and after mining when the fully mechanized coal mining face passes through the collapse column are investigated. Here, we present the constructed program in FISH language, render the damage variable in FLAC3D to establish the numerical model, and complete the numerical calculation. The 10–115 working faces that pass the collapse column at a coal mine in Tuanbai are identified as the research object. The floor failure is numerically simulated to assess the damage. The following results were obtained: the failure depth of the full floor is stabilized at 14.6 m; the maximum failure depth of the floor near the collapse column is 18.2 m; and the stress concentration coefficient is 1.27 times greater than that of normal mining. The calculated depth failure of the floor of the working face without structural defects is 14.6–14.7 m based on the Hoek–Brown criterion. With the collapse column, the failure depth of the floor is 16.8–17.8 m. According to the water injection test, the maximum failure depth of the floor is 18 m. The three derived values agree well with one another.

1. Introduction

In northern China, collapse columns are widely distributed in the Permian–Carboniferous coalfields. When the compaction and cementation of the fillings in the collapse columns are poor, the columns may be activated and transformed into Ordovician limestone water inrush channels under the action of mining and other external factors, which can lead to water disasters. Collapse columns can directly affect the safe and efficient operations of coal mines [14]. Prior to the full mechanization of a mining face that passes a collapse column, the influence of this collapse column on mining stress and the influence of mining stress on the collapse column and floor failure of the working face both need to be established. Many scholars have studied the stress distributions of collapse columns and their surrounding rocks and the related failure depth of floors. Gao et al. [5] summarized the law of floor water inrush and proposed the preferred plane theory of water inrush. Xu et al. [6] analyzed the mechanism and criterion of the collapse column’s activation and water conduction under weak runoff conditions to improve the accuracy of water inrush forecasts. Wang et al. [710] explained the burst of water inrush from the floor based on cusp catastrophe theory by using mathematical mechanics. Finite difference, finite element, and other numerical simulation methods have been applied in recent years to investigate water inrush in coal seams. Wang and Yin [11] numerically simulated the stress-strain characteristics of the surrounding rocks of a mining field under the influence of collapse columns by using FLAC3D with the finite difference method. Wang and Song [12] studied the randomness and connectivity of coal floor element failures with the renormalization group method. Liu and Xiong [13] and Hao et al. [14] conducted numerical simulation and found that the passing collapse column in a fully mechanized working face mainly undergoes three stages, and the greatest impact can be observed at the side of the goaf area, which are normally formed by lagging water inrush rather than the abovementioned working face. Wang and Park [15] explained the failure and water inrush mechanisms of a coal floor in a goaf area by conducting simulation analysis in FLAC3D. Many domestic and foreign scholars have studied the fracture evolution and the weakening of materials and rocks by applying damage and fracture mechanical methods. Subsequently, the corresponding constitutive relations and strength criteria that can be applied to the collapse column model are established. Zhu and Sun [16] applied self-consistent theory to deduce the equivalent flexibility tensor of rocks and the instantaneous modulus described by Krajcinovic and Fonseka [17], who previously proposed the idea of change in flexibility of damaged bodies due to macromechanical damage effects. Xu et al. [18] and Liu et al. [19] derived the constitutive relations of fractured rock masses from the reciprocal theorem and established the damage evolution equation and the probability model with a strength criterion by employing fracture mechanics. Both authors also verified the applicability of the water inrush law of collapse columns through theoretical analysis and numerical simulation. However, no interrelations were established between the two studies.

Here, we simulate a real terrain to construct a program in FISH language, render the damage variables in FLAC3D, and complete the calculation, in which the damage and the stress and failure depth of the floor are considered, by numerically simulating the floor failure.

2. Engineering Background

The minefield in Tuanbai is located at the central part of the Huoxi coalfield at the southern part of Huozhou mine. The minefield is widely covered by Quaternary loess. The exposed area of the bedrock is approximately 20% of the total area of the minefield, and it is mainly distributed in the ridge and gully areas in the midwestern area. The stratum is the lower Shihezi formation type of Lower Permian (P1x) and upper Shihezi formation type of Upper Permian (P2s). The floor elevation of #10 coal mine in the past 5 years is from +380 m to +460 m, and the average thickness is 6 m. The water level elevations of the K2 and O2 aquifers are 480 m and 520 m, respectively. Therefore, #10 coal is mined under pressure relative to the two aquifers. The fracture water of the Ordovician limestone karst is the main factor affecting #10 coal mine. In addition, the lower coal group is similar to the Ordovician limestone. The average interval between #10 coal mine and the Ordovician limestone is 36 m. The aquiclude is mainly the Taiyuan and Benxi formation types.

Collapse columns are well developed at the Tuanbai minefield with an average of 29.6 columns per square kilometer. They are rarely seen on the surface because of the Cenozoic strata coverage and the slip-and-collapse effect of rocks; in other words, they are mostly located underground and unexposed. The cross sections of the collapse columns are mainly elliptical, but some are circular or irregular. The cross-sectional areas differ in terms of dimensions. The long axis is 20–98 m (mostly 40–60 m), and the short axis is mostly 20–40 m. The largest cross-sectional area is approximately 15,000 m2, and the smallest is 100 m2. Most of the cross-sectional areas are approximately 1000 m2. The presence of collapse columns suggests the demarcation of noncoal zones that can affect the layout of the mining area and its working face. Consequently, tunneling costs and the degree of difficulty of realizing a fully mechanized coal mine are greatly increased. Collapse columns may also indicate a water drench phenomenon, which implies water inrush, and thus can complicate the hydrogeological conditions of the mine.

3. Realization of the Damage Theory in FLAC3D

A rock is a complex material composed of several mineral particles, pores, and cements. Under long-term geological conditions, rock materials inevitably acquire micropores and microcracks that are randomly distributed in different sizes and shapes. In damage mechanics, the macroscopic mechanical effects caused by the formation and development of internal defects of materials and the process and law to explain the material damages are mainly studied. The internal state variable called “damage variable” is also introduced to describe the mechanical effects of materials that contain microdefects and the mechanical behavior of the damaged materials. Furthermore, the morphology, distribution, and evolution law of all types of damages are explored, i.e., from the microscopic scale by using the averaging method to reflecting the findings in terms of their corresponding macroscopical mechanical behavior.

In general, the study of damage mechanics consists of four stages:(1)Selecting the appropriate damage variable.(2)Establishing the damage evolution equation.(3)Constructing the constitutive model by considering the material damage.(4)Solving the stress-strain and damage values of the different points of materials based on their initial conditions and boundary conditions and then evaluating the damage states of these points based on the damage values. If the damage reaches a critical value, then the point is considered damaged. Then, similar calculations are rendered based on the distribution of the new damage and boundary conditions, repeated, and finally terminated when the damage criterion of the component is reached.

Figure 1 simplifies the total stress-strain curve of rocks of the double-line type. The AB segment is the elastic stage without damage, while the BC segment is the strain-softening stage with isotropic damage. Thus,where and are the stress and strain of the peak point, respectively; and are the stress and strain of any point in the softening segment, respectively; is the final stable strain; and is the residual stress.

Based on the stress equivalence principle of damage mechanics, we obtain the following equation where is the damage variable and is the elastic modulus:

The constant , also known as the damage modulus, is introduced and can be expressed as

The expression of the damage variable is obtained by the Gauss formula and the motion equation of the Lagrange finite difference method:

A real material that is not uniform and continuous can still be depicted with an initial damage. This damage variable should have a value larger than that in formula (4). After modifying formula (4), the damage evolution equation is obtained as follows:

To embed formula (2) into the numerical computation model in FLAC3D, its constitutive FISH language needs to be programmed to be able to record the strain values in different units after the array is changed every 100 steps. Then, to obtain the strain value and the stable strain value of the peak point after their corresponding units are changed, the IDs of these units also need to be changed. In other words, the strain values should complement the altered arrays, and the strain change values of the units have the same IDs. The results of the modeling are substituted into formula (5) to obtain the new damage variable. Subsequently, the new damage variable is substituted into formula (2) to complete the numerical calculation. Finally, the numerical simulation of the floor damage under the condition of damage is realized.

4. Numerical Simulation and Prediction

4.1. FLAC3D Model Establishment

The geometric size of the numerical simulation model that contains the collapse columns is 350 m × 240 m × 170 m. The model is divided into seven layers based on the actual coal seam situation at the Tuanbai minefield. The black fifth layer represents the coal seam with a thickness of 6 m. The aquifer is 60 m from the model floor, and its thickness is 38.6 m. The positions of the model and the collapse columns are shown in Figure 2. The coordinates of the center position are (175, 185). Mesh encryption is conducted on the collapse columns during model establishment. The pore radius of the upper end of the collapse column is 10 m and that of the lower end of the collapse column is 20 m, including the surrounding fracture zones. The development height of the collapse column is approximately 115 m and directly passes the low coal group and connects the key layer and the aquifer. The collapse column is also divided into seven layers, in which the top and bottom layers are attributed the same properties as that of the rock formation. The collapse columns at the coal seam also have similar properties as those of the lower-level collapse columns. Thus, the column contains four layers. The collapse is in the shape of a column from top to bottom and surrounded by fracture zones. The rock mechanical parameters and the physical parameters of each rock stratum are shown in Table 1. The physical and mechanical parameters and the fluid mechanical parameters of the fracture zones of the collapse column are shown in Table 2.

The weight of the modeled overlying strata is measured based on the geostress of the Tuanbai coal mine. The maximum principal stress at the level of 423 m (depth of 306.28 m) is 11.13 MPa. Therefore, the self-weight stress of the overlying strata is taken as q = 10 MPa. The stress is loaded into the model in the form of uniform distribution. The pore water pressure in the bearing aquifer at the bottom is set with a fixed value. On the basis of the field test, the stress is 1.2 MPa. The mining engineering face starts from the position 45 m from the right boundary of the model and advances toward the left. The working face is 40 m wide, and the boundary of the model is 40 m from left to right. The simulation is advanced 10 times in total at 20 m each period. New rock materials are added to the goaf area from the previous section at each advancing step to simulate the collapse of the slope. When the working face has advanced to 140 m, its right end is expected to directly pass the collapse column.

4.2. Fluctuating Line Fitting of the Working Face

To simulate the real terrain, the fitting function of the terrain contour is obtained. The baseline is established on the basis of the trend of the terrain. The vertical elevation from the topographic fluctuating curve to the baseline is measured. Furthermore, 100 groups of data points are taken. The horizontal and vertical axes are designated as the X and Y axes, respectively. The data points are recorded in Microsoft Excel. Polynomial fitting is conducted with MATLAB, and the fluctuations in the vertical direction are reduced in proportion to obtain the six-order polynomial approximation fitting function of the topographic fluctuating curve.

FISH language is used for FLAC3D programming. The obtained fitting function is used to calculate the height of a node in the z direction on the surface of the model. Then, the zooming factor k is calculated based on the original height of the node. Finally, the heights of all nodes in the z direction at a certain position are zoomed by the factor k, thereby resulting in the fluctuation of the model strata along the trend that conforms with the fitted terrain fluctuation curve. As a result, a real terrain is simulated. The modified model can effectively simulate real terrain trends, and the numerically simulated results can approximate the real situation (Figure 3).

4.3. Analysis of Simulated Results

To study the scope and depth of the floor of the coal seam as an effect of the mining process, we compare two conditions (normal mining and mining under collapse columns) by using FLAC3D. The failure scope of the floor caused by mining and the eventual disturbance of collapse columns and the original rock stress near the working face are analyzed.

4.3.1. Analysis of the Vertical Stress: Floor Distribution

Figure 4 shows the change in the vertical stress of the surrounding rocks as the working face advances. The following discussions can be derived from the figure:(1)When the working face is more than 25 m away from the collapse column, the stress field of the floor undergoes stress concentration and pressure relief with the advancement of the working face, which is similar to the situation without the collapse column. The supporting pressure at the front end of the working face is concentrated, the pressure relief area of the floor is gradually increased, and the rear end of the working face is filled and compacted and then gradually restored into the original rock stress states.(2)When the working face begins to undergo excavation, the rock area located away from the collapse column is transformed into a pressure relief area. However, with the advancement of the working face (i.e., the working face is 5 m away from the collapse column), the stress field near the working face begins to markedly affect the collapse column. The stress concentration coefficient is 1.15 times higher than that of normal mining and subsequently increased to 1.27 times when advanced to the collapse column. Then, the collapse column is affected by the support pressure at the front end of the working face, the pressure relief area of the floor, and the overall supporting pressure. Finally, as the collapse column passes, the working face reverts to normal recovery.(3)Owing to the relatively high stress concentration coefficient, the stress relief area of the surrounding rocks and the influencing depth of the pressure relief area are both increased when the working face advances from the closing to the passing of the collapse column.

4.3.2. Failure Characteristics of the Plastic Zone of the Floor

Under the influence of mining, the failure of the floor changes with the advancement of the working face. The failure range of the plastic zone of the floor also continues to increase. Figure 5 shows the change in the plastic zone during normal recovery and with the passing of the collapse column with different advancement distances. In this paper, only three stages of passing are presented. Figure 6 shows the following curve change of the floor failure depth when the working face is close to the collapse column:(1)Working face 25 m ahead of the collapse column(2)Working face 5 m ahead of the collapse column(3)Working face above the collapse column(4)Working face 10 m behind the collapse column

The following discussions can be derived from Figures 5 and 6:(1)The area near the front end of the working face undergoes the deepest failure at the floor of the working face. The failure of the middle part of the working face is relatively low. Thus, a failure zone in the inverted saddle shape is formed.(2)The failure depth of the floor is 9.8 m. When the working face advances to approximately 60 m, the failure depth of the floor becomes 14 m, and the failure rate is gradually reduced. The failure depth is finally stabilized at 14.6 m. The failure modes of the area at each advancing step for the lower working face are the same.(3)When the working face is 25 m away from the collapse column, the depth of the floor is stabilized at 14.6 m. Then, the working face advances to the collapse column. The closer the working face is to the collapse column, the greater is the failure depth of the floor. The failure depth is 17.8 m and located 5 m away from the working face. The deepest failure is observed near the edge of the collapse column.(4)When the working face advances to the location above the collapse column, the failure depth of the floor in front of the working face is the largest at 18.2 m. The failure depth of the middle part of the collapse column is relatively small at 14.7 m and is equivalent to the case when the working face is more than 25 m away from the collapse column.(5)The failure depth of the floor is restored to that of the full floor when the working face passes the collapse column.

5. Theoretical Calculation

Hoek–Brown [19] deduced the relationship between the ultimate principle stresses of rock and rock mass failure by trial and error after deriving theoretical research results and then obtained practical experience on rock property characterization based on the analysis and correction of Griffith theory. The Hoek–Brown strength criterion can be expressed aswhere and are the maximum and the minimum principle stresses under rock failure (i.e., compressive stress is positive) and is the uniaxial compressive strength of a rock, which can be measured with uniaxial pressure and point load tests. To determine the softness and hardness of the rock mass, we let to range from 0.001 to 25 and (breaking degree of rock mass) to range from 0 to 1. According to theoretical research, is related not only to the structure of a rock mass but also to its internal friction angle, compressive strength, and tensile strength. By contrast, is only related to the structure of the rock mass.

Hoek–Brown stressed that the material constants and of a rock can be estimated based on Bieniawski’s rock classification index (i.e., RMR). For the disturbed rock mass,

For the undisturbed rock mass,where is the value of intact rock, which can be determined with a triaxial test. In engineering, is 10 for mudstone, shale, and slate, while is 15 for sandstone and quartzite. According to the calculation, the failure depth of the full floor is approximately 14.6–14.7 m. Under the condition of the floor with a collapse column, the failure depth of the floor is increased by almost 3 m and can reach 16.8–17.8 m.

6. Water Injection Test

The water injection test is conducted in accordance with the seepage expression of the protective layer (aquifer) [20]:where and are the permeability coefficients of the aquifer and the aquifuge, respectively; and are the water pressure and the hydraulic gradient, respectively; and and are the seepage length and the section area, respectively. When the rock fracture is opened because of the mining pressure, becomes relatively larger and becomes relatively smaller. Thus, is increased under the same pressure. Conversely, when the rock fracture is closed because of mining pressure, is decreased. When the floor undergoes water conduction because of the ore pressure, initial and should be 0. Therefore, the depth when initial is 0, a value near 0 or no obvious sudden change in pressure relief time during water pressurization represents the direct failure depth of ore pressure.

To safely exploit the #10 coal seam at the Tuanbai coal mine, a field test has been conducted on the site with collapse columns in the first mining 10–115 working faces. The testing results may be used as basis for the safe mining of working faces in the future.

6.1. Test Methods

According to the “Drilling Design for Fracture Depth Detection of Collapse Columns,” five holes can be drilled at the bottom plates of working faces #10–#115 (Figure 7). The end holes are placed at the seam floors at 17.5 m, 18.5 m, 20.5 m, 21.5 m, and 22 m, which are denoted as #1, #2, #3, #4, and #5, respectively. Outcrop pregrouting is accomplished after grouting to ensure the decompression time of longer than 5 min. Water is added into each hole. Uprising time, maximum hydraulic pressure, decompression time, decompression range, and working face distance are recorded. The results of data analysis are shown in Figure 8.

Figure 8 can be interpreted as follows:(1)The sudden changes of decompression time were observed at #1 and #2 holes on April 17th. The time decreased in the range of 1–5 minutes and was stable at approximately 1 minute. The drill holes ending at #1 and #2 were damaged. In this circumstance, the working face was 67.8 m away from the #11 working face and 17.8 m away from the horizontal position of the drill holes ending at the #1 drilling hole. The angle between the fissure point and the working face was nearly 45° and conformed with the fracture development law at the bottom plate.(2)The decompression times at #3, #4, and #5 did not change dramatically and were maintained at >5 minutes. The fissure zone was speculated to be shorter than 20 m.(3)The drilling data of the #3 hole fluctuated on April 19th, but they were subsequently normalized. Following the principle of statistics, these observation data were ignored in the analysis.

On the basis of the comprehensive analysis combined with the water injection test results, the disturbance depth at the bottom plate of the working faces reached 18 m.

7. Conclusions

(1)In the process of model calculation, to reflect the influence of the mining on the properties of the surrounding rocks more comprehensively, we embedded the model with failure conditions in FLAC3D in FISH language. As depicted by the overall stress-strain curve of the tested rock samples, the failure model is the double-line type.(2)The stress distribution of the floor is simulated. When the working face is more than 25 m away from the collapse column, the stress field of the floor with the collapse column is the same as that of the full floor. On the basis of the characterized distribution of the plastic zone of the floor, the failure depth of the full floor is stabilized at 14.6 m. Moreover, when the working face advanced to a location above the collapse column, the failure depth of the floor in front of the working area is the largest at 18.2 m.(3)The failure depth of the floor is calculated with the Hoek–Brown criterion. According to the results, the failure depth of the floor of the working face is 14.6–14.7 m with no structural defects. In terms of the floor with collapse column, the failure depth can increase by 3 m and reach 16.8–17.8 m.(4)The failure depth of the floor of the 10–115 working faces at the Tuanbai coal mine is investigated with the water injection test method. The measured failure depth is approximately 18 m. A comparison of theoretically calculated and numerically simulated values shows that the difference is minimal.

Data Availability

The others can access the data supporting the conclusions of the study from this research article. The nature of the data is the laboratory experimental data, the field observation data, and the theoretical calculation data. The laboratory experimental data used to support the findings of this study are included within the article, mainly the mechanical parameters in Tables 1 and 2. The field observation data used to support the findings of this study are included within the article, and the change of water inflow on different mining face distance is included in Figure 8. The theoretical calculation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

All the authors contributed to publishing this paper. Jinlong Cai conceived and designed the experiments, performed the experiments, and wrote the paper; Wensong Xu analyzed the data; and Ming Tu revised and reviewed the manuscript.


This work is part of the projects financially supported by the Chinese National Natural Science Foundation (NSFC) (Grant no. 51574007). The authors gratefully acknowledge the financial support received from NSFC.