Abstract
In situ fluidization mining is a promising method which shows considerable potential for reducing the ecological influence impacts of coal mining. During fluidization mining, coal gangue will be separated from the ores, filled back, and returned to the excavation region, which, depending on the filling rate, greatly reduces the stratum deformation of strata, depending on the filling rate. Numerical simulations can assist engineers in evaluating stratum deformation and optimizing affordable designs with affordable costs. Accordingly, in this work, a continuous–discontinuous element method is adopted, which is a coupled particle–block finite-discrete element method with an explicit integration formulation, which is adopted for simulating the coal in situ fluidization mining–filling processes. The theory concept of an “equivalent mining height” is utilized to define the backfill mining height. The results partly reveal the influences of the filling rate on the stratum deformation, surface subsidence, and rock burst phenomena, all of which further indicate the advantages of the fluidization mining method.
1. Introduction
With the increasing economic demand for coal resources, the long-term, high-intensity excavation of coal at shallow depths has exhausted the available coal resources near the surface; as a result, coal mining has gradually moved into deeper strata [1], while the mining scope of the current solid mineral mining method is relatively shallow [2, 3]. Traditional coal mining yields large volumes of waste known as gangue, the discharge of which can occupy land, pollute the environment, and endanger surface buildings. To reduce ecological damage, a novel method involving the in situ fluidization mining of deep coal [4] in which the coal gangue is separated from ores and placed back into the excavation region was proposed. Mining, selection and filling as one [5–9], backfill bodies, and coal pillars can help to control the settlement of strata and mitigate surface subsidence. Hence, these approaches can protect the environment while mining and converting deep coal in situ.
Field monitoring, theoretical analysis, and numerical simulation methods have been widely used to control surface subsidence, reduce ore pressures, and effectively reduce the impacts of waste on the surface. Zhou et al. [10] proposed the concept of the filling ratio and studied both the factors affecting the filling ratio. Sun [11] established a solid compaction and backfill mining dynamic model based on the key stratum theory for controlling rock layers, thereby providing a theoretical basis for determining the design of a solid backfill mining working face. Xu et al. [12] conducted numerical simulations to analyze the movement and deformation of the overburden in the stope and concluded that backfill mining can effectively control the failure degree of the overburden and reduce the stress concentration and influence range. Liu et al. [13] proposed a mechanical model of an elastic plate to analyze the coordinated instability of a coal pillar and roof and then specifically described an approach for evaluating the stability of the goaf and the design process for the filling parameters. With parameter studies, Zhao et al. [14] reported how the mechanical properties of the filling body affect the movement and stress of the roof. Feng et al. [15] proposed the concept of a backfill mining collaborative control system based on the specific geological conditions of coal mines and analyzed the supporting characteristics of the overburden strata. Yu et al. [16] backfilled a mining roadway, established a mechanical model of roof settlement, and analyzed its control effects on ground settlement compared to traditional mining methods. Lu et al. [17] analyzed the stress distribution and deformation of backfill bodies through on-site ground monitoring and studied the boundary effect of backfill bodies. Zhao et al. [18] monitored the in situ stress and deformation of the working face and showed that the solid backfill technique could effectively reduce the failure degree of the overburden rock and control roof subsidence, thereby realizing the effective control of deformation in the overburden.
Considering complex geological environments, compared to experimental studies and on-site monitoring, numerical simulations can more economically identify the failure patterns and deformations of strata than can theoretical analyses. However, it is challenging for the conventional finite element method (FEM) to consider the cracking processes of strain-softening materials in a discrete manner. The increasing number of discontinuities affects the numerical stability and efficiency; as a result, in these types of cases, FEM models may not converge. Moreover, although some novel numerical methods, such as the cracking element method [19–24], cracking particle method [25–29], and peridynamics [30–34], can help to capture the initiation and propagation of fractures, in numerical simulations of fluidization mining, masses and gangue should be modeled by blocks and particles, respectively, which can fragment, requiring a numerical tool capable of handling the increased complexity of fragmented blocks and particles.
In this paper, the continuous-discontinuous element method (CDEM) coupled with the particle-block method (hereafter CDPB) is adopted for dynamically simulating the in situ replacement of coal in the goaf with gangue and studying the characteristics of the stope ore pressure and rock deformation. CDPB is a hybrid finite-discrete element method built in an explicit iterative numerical framework. The finite element procedure is used for small deformations of the system, and the discrete element procedure is used for large displacements and large deformations of the model. The concept of the “equivalent mining height” is used to define the in situ converted mining height of deep coal seams. By changing the filling ratio of gangue, surface subsidence can be controlled, which provides a theoretical basis for reducing the mine pressure and managing the resulting effects.
The remainder of this paper is organized as follows: the numerical method and constitutive model employed in this paper are presented in Section 2. Section 3 describes the simulation of in situ fluidization mining. In Section 4, the simulation results are presented and discussed. Finally, the concluding remarks are given in Section 5.
2. Numerical Method and Constitutive Model
2.1. Continuous-Discontinuous Element Method
The CDEM [35] is an explicit iterative numerical method that couples finite elements to discrete elements based on the Lagrange function, as shown in Equation (1). The numerical model is composed of block elements and interface elements: the block elements represent the continuous material characteristics, such as their elasticity/plasticity and damage, whereas the interface elements represent the discontinuous material characteristics, such as fracturing. Virtual springs are introduced between interface elements to transfer mechanical information, as shown in Figure 1. By setting the fracture criterion and mechanical parameters of these springs, the initiation, development, and propagation of multiple cracks can be simulated. Throughout the analysis, the force balance of the system is calculated by setting the unbalance rate: where is the nonconservative force of the system, is the Lagrange function, is the displacement of an element node, and is the node velocity.
The control equation of CDEM is where , , and are the acceleration, velocity, and displacement matrices of the nodes of block elements. and are the displacement and velocity matrices of the nodes of virtual interfaces. , , , , , and are the unit mass, unit damping, unit stiffness, contact surface, contact surface damping, and node external load matrices, respectively. Based on the incremental method, the CDEM uses the explicit Euler predifference method to solve the problem. The unbalance degree of the system force is characterized by the unbalance rate. Equation (2) is solved in three steps: (i) calculating the deformation and damping forces of all elements, (ii) calculating the contact and damping forces of all interfaces, and (iii) calculating the external forces, accelerations, velocities, and displacements of all nodes.
2.2. Elastic Damage-Fracture Constitutive Model
In the numerical calculation, a linear elastic constitutive model is applied to the finite elements [36]. The resulting linear elastic constitutive model of finite elements is expressed by the incremental method: where is the stress tensor, is the incremental stress tensor, is the incremental strain tensor, is the bulk strain increment, is the bulk modulus, is the shear modulus, is the next step, and is the current step.
The Mohr–Coulomb brittle fracture model is applied to the virtual interfaces to calculate the deformation and fracturing. The incremental method is used to calculate the normal and tangential test contact forces in the next step on each virtual interface: where and are the normal and tangential contact forces, respectively; and are the normal and tangential contact stiffnesses, respectively, on the unit area (unit: Pa/m); is the area of the virtual interface; and and are the normal and tangential relative displacement increments, respectively. Formula (5) and Formula (6) are used to judge tensile failure and correct the normal contact force and tensile strength.
If then where is the tensile strength of the virtual interface at the next moment (unit: Pa). Formula (7) and Formula (8) are used to judge shear failure and correct the tangential contact force and cohesion.
If then where is the internal friction angle of the virtual interface and and are the cohesion of the virtual interface at the initial moment and the next moment, respectively (unit: Pa).
2.3. Validation
The Brazilian disk test is used to validate the numerical tool. We used the same setup provided in [37], where the diameter and tensile strength of the disk are taken as 6 cm and 6 MPa, respectively. The disk is built with particles and blocks. The peak load per unit thickness can be determined by
The obtained force-displacement curves and failure patterns are shown in Figure 2. The peak load is generally close to the analytical solution, and the failure pattern is also agreeable.
(a)
(b)
2.4. Coupled Particle-Block System
This work considers the coupling of block and particle discrete elements [38]. The coupling processes include (i) detecting the particle-block contact and (ii) calculating the contact force. The contact force is obtained by the incremental method: where and represent the normal and tangential displacement differences (cracking openings), respectively, between the contact particles.
The formula for calculating the contact torque is as follows: where and are the torque and bending moment, respectively; and are the moment of inertia and polar moment of inertia of the contact surface, respectively; and and are the incremental differences in the torsion and bending angle between particles, respectively.
Particle-block contacts are detected in two steps: preliminary detection and accurate detection. The preliminary detection step is achieved by the subspace method. The scope of retrieval can be reduced by the mapping relationship among particles, cells, and grid areas. The accurate detection step locates the contact between a particle and the boundary of the finite element. The point-edge contact model is adopted, as shown in Figure 3. The creation of a point-edge contact should satisfy two conditions simultaneously: first, the distance from the particle centroid to the boundary edge (as shown in Equation (12)) must be less than or equal to the particle radius ; second, the projected location of the particle centroid on the boundary edge is located inside the edge , where is the relative position vector from to particle and is the outer normal vector of the edge. Once the particle is in contact with an edge, normal and tangential springs are automatically created. The interpolation coefficient of the contact point is solved by Equation (13):
2.5. In Situ Fluidization Mining
In conventional mining, coal mining and production processes are independent. The coal and coal gangue (solid waste) will be transported together to the surface through a pipeline. Then, a “gangue mountain” will be formed on the surface before backfill, potentially causing pollution. In addition, the costs of transportation and backfill will increase dramatically for deep mining. In situ fluidization mining is a method combining the mining and production processes that separates coal and gangue in situ as an economical and eco-friendly mining method. Moreover, compared to conventional caving mining, since gangue is backfilled in site, fluidization mining can reduce the stratum deformations and mitigate the fracturing of the roof.
In situ fluidization mining includes mining, separation, and filling processes (see Figure 4). In the process of unattended mining, an underground tunnel drilling machine is used for mining treatment. The minerals are separated and sorted in the separation tank. Under the action of gravity, the slag mixture is selected to form the backfill material, which is backfilled to the excavated area. In the next sections, considering the equivalent mining height theory, the surface subsidence and ore pressures of fluidization mining will be studied and compared to conventional mining methods.
2.6. Equivalent Mining Height Theory
In the process of filling the mining face with gangue, the coal is converted in situ. While the coal to be mined remains in front of the working face, solid waste materials such as gangue are filled into the goaf, occupying parts of the stope. The actual coal seam thickness is reduced. This process can be regarded as extremely thin coal seam mining. The equivalent mining height [40] depends on the roof displacement and filling body compression. As shown in Figure 5. The maximum excavation height is the equivalent mining height after backfill mining. Furthermore, the deformation of the overlying strata, surface subsidence, and ore pressure of in situ fluidization mining depend on the filling ratio of gangue.
The equivalent mining height of gangue filling is where is the amount of roof subsidence known in advance, is the amount of nonconnection, and is the amount of compressed filling body, which can be obtained by where is the total mining height after backfill mining (m), and is the compression ratio of the filling body. depends on the porosity: where and are the porosities of the filling solid after initial and final compaction, respectively. Hence, is
3. Numerical Simulation of the In Situ Conversion and Mining of Deep Coal Resources
3.1. Engineering Background
The Kailuan Group Tangshan mining area is located in Tangshan, Hebei Province. Within this area, coal seams 5, 8, and 9 are thick, stable, and recoverable. These coal seams are located at depths of 700800 m and thus are considered deep coal seams. The in situ converted goaf is located in the No. 9 coal seam, and the thickness of the coal seam is 3.5 m. The ore pressures are imposed by the buildings, railways, and water bodies at the surface above the mine field, approximately 171 million tons. With the conventional mining method, 800,000 tons of coal gangue is discharged annually to the surface, occupying a great amount of land.
Based on the working face in the above coal seam, considering the height of each layer and boundary effects, the numerical model established in this study is 300 m long and 111 m high. The model and mesh of the example are shown in Figure 6. Two measuring points are prescribed in the model at a distance of 1.5 m below the roof and 100 m from the left and right boundaries of the coal seam. The material parameters of the mining strata are shown in Table 1, where the interface parameters are set according to stratum parameters, including stiffness parameters and strength parameters. The stiffness is 10 times that of the rock layer, and the strength is 3 times that of the rock layer. Gmsh software is used to discretize the model into 13,849 nodes and 27,255 elements, assuring the accuracy and efficiency.
(a)
(b)
3.2. Simulation Process
The numerical simulation includes three stages. The first stage is the elastic stage, in which boundary normal constraints are applied to the left and right sides and the bottom of the material. Then, the in situ stress exerted by the overlying strata is initialized at the top of the model. The block element and interfaces are both linear elastic materials. In the second stage, the material model adopts the Mohr–Coulomb model, whereas the virtual interfaces adopt the brittle fracture model. In this elastic-plastic stage, the static stress field under the influence of gravity is solved by using the virtual mass method, and the stress balance is obtained. The third stage is the coal seam mining and gangue filling stage, considering the coupled particle-block contact effects. The goaf is created by excavating the coal seam from left to right at an interval of 5 m for 40 iterations. In the next section, if it is not specified, the No. 9 coal seam will be mined. The coal gangue is filled into the goaf after each excavation to replace part of the stope. Therefore, coal pillars and the gangue filling body support the overlying strata and help to reduce ground subsidence and mitigate rockburst. In the numerical studies, we also consider different filling ratios of coal gangue and compare fluidization mining to conventional caving mining. The gangue parameters are shown in Table 2. After excavating the working face, the coal pillars in certain sections are set as empty models, and then, the goaf is backfilled with coal gangue particles with different filling ratios, supporting the overlying strata, as shown in Figure 7.
(a) First excavation
(b) The first coal in situ replacement
(c) Second excavation
(d) All coal in situ replacement
4. Results and Discussion
4.1. Deformation of the Roof
The strata overlying the backfill exhibit compression deformation in the long term. However, because of the frictional forces among the gangue particles inside the filled goaf, the grain-scale deformation of the filling body gradually becomes stable. Initially, there is no supporting effect on the roof, and the roof exerts no pressure on the gangue filling body in the beginning. When the roof subsides to the equivalent mining height, part of the roof is supported by the filling body and coal pillar, and some pores among the particles in the gangue filling body are compressed, restricting the settlement of the roof, as illustrated in Figure 8.
The supporting performance of different filling materials depends on (i) the deformation moduli of the backfill material and (ii) the filling ratio. Considering different elastic moduli of the backfill material and filling ratios, the long-term roof subsidence is illustrated in Figure 9.
(a)
(b)
According to Figure 9, the elastic moduli of the filling materials have relatively small influences on the roof subsidence, while the filling ratios have much greater influences. In each curve, there is a platform part at the bottom where the coal gangue filling material has reached the maximum compression and the equivalent mining height has reached.
4.2. Ground Subsidence
The ground subsidence depends on the filling ratio and mining depth. The filling ratios of 0%, 70%, 80%, and 90% are considered for investigating their influences on ground subsidence, in which a 0% filling ratio is close to that of traditional caving mining. The long-term ground subsidence is shown in Figure 10(a), indicating that higher filling ratios help to control the ground subsidence. Considering monitoring point 1 (see Figure 6), the evolution of its displacement is shown in Figure 10(b). Figure 10(b) shows that the stabilization of the roof for all filling ratios will take more than 10 days. Moreover, for the case with a filling ratio of 0%, the displacement curve is not very smooth, mainly caused by serious fragmentation of the bent roof.
(a)
(b)
To investigate the influence of excavation depth on ground subsidence, excavations of coal seams No. 5 and No. 8 in Table 1 are considered. The ground deformation (vertical and horizontal displacement) curves are shown in Figure 11. The results indicate that the ground displacements decrease with increasing mining depth.
(a)
(b)
4.3. Deformation and Failure Pattern of the Coal Seam
Considering a filling ratio of 90%, the evolutions of the vertical stress component at monitoring points 1 and 2 (see Figure 6) are shown in Figure 12. The results indicate that at both points, the peak stresses of fluidization mining are smaller than those of caving mining, which helps to reduce rockburst events.
(a)
(b)
Considering a filling ratio of 90%, the contour plot of is shown in Figure 13(a), indicating that the roof breaks like a deep beam under self-weight. For different filling ratios, the distributions of of the coal pillar are shown in Figure 13(b), indicating that the supporting stress decreases with increasing filling ratio.
(a)
(b)
The rupture degree is defined as the ratio of the area of activated interface elements to the area of all interface elements implemented in numerical calculations. For different filling ratios, the evolution of the rupture degree is shown in Figure 14. The rupture degree decreases with increasing filling rate.
5. Conclusion
Coupled particle-block elements built in the framework of CDEM are used for simulating the in situ fluidization mining processes of a deep coal seam where the coal gangue is filled back. The responses of the strata are analyzed and compared to those of traditional caving mining. Some conclusions are drawn: (i) to control the roof and ground subsidence, the filling ratio of coal gangue has a greater influence than its mechanical properties. When the filling rate increases from 70% to 90%, the roof settlement decreases by 1 m, and the ground subsidence decreases by 0.12 m. (ii) In situ fluidization mining greatly helps to control the fragmentation and deformation of strata. When the filling rate increases from 70% to 90%, the peak stress decreases by 2.5 MPa, and the rupture degree decreases by 1%.
Data Availability
The 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
The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (NSFC) (52178324) and National Key Research and Development Project of China, The Ministry of Science and Technology of China (2018YFC1505504).