Constitutive Model of Solid Backfill Materials and Numerical Simulation of Overburden Movement and Deformation in Backfill Mining

Solid backﬁll mining is an eﬃcient and environmental-friendly coal mining technology, which can eﬀectively solve the problems of coal gangue pollution, water resource loss, and surface subsidence. Based on the mechanical behavior of backﬁll materials in the compaction process, volume strain was used to express the deformation modulus, and a constitutive model of backﬁll materials was proposed in this study. The ABAQUS UMATwas used to develop the numerical calculation subroutine of the model. Finally, the rationality of the model was veriﬁed that simulated stress-strain curves of the backﬁll materials during the compaction process agree well with experiments. Based on the proposed constitutive model, the inﬂuence of three factors (the initial compaction rate of the ﬁlling body, the mining height, and the mining depth) on the key strata and surface subsidence was analyzed systematically. The results show that the initial compaction rate and the height of coal seams have signiﬁcant inﬂuences on surface subsidence. When the thickness of topsoil is only changed and the structural composition and lithology of overburden are not changed, the mining depth has little inﬂuence on surface subsidence, but a signiﬁcant inﬂuence on surface subsidence range. The inﬂuence of mining height and mining depth on the deformation of key strata of overburden and surface subsidence is approximately linear, while the inﬂuence of the initial compaction rate is nonlinear.


Introduction
e rapid development of China's economy relies heavily on coal resources. With the gradual reduction of resources, the intensity of coal mining is gradually increasing. In traditional coal mining, coal gangue will be generated during the roadway excavation and coal washing, accounting for about 15% of the total weight. e accumulation of coal gangues has caused a series of environmental problems [1][2][3], such as land occupation, water pollution, and air pollution after spontaneous combustion [4][5][6][7]. In addition, coal mining can also cause the fracture of key strata [8][9][10][11], resulting in the loss of groundwater resources [12][13][14] and surface subsidence [15,16]. Solid backfill mining is an environmental coal mining technology, in which coal gangue is used as the main filling body to replace the underground coal resources [17,18]. is method can be used to consume the existing coal gangue, alleviate the environmental problems caused by the accumulation of coal gangue, reduce the disturbance of coal mining on overburden, and prevent the fracture of key strata and greatly reduce the surface subsidence [19][20][21].
In solid backfill mining, backfill materials are mainly composed of coal gangue, fly ash, and other industrial wastes. e mechanical properties of backfill materials determine the control ability of backfill mining on the movement and deformation of overburden. At present, the studies on mechanical behavior of backfill materials have been systematically studied, including the axial stress-strain curve of the compacted cylinder [22,23], the creep relaxation [24], the influence of particle size ratio on the mechanical behavior [1], and the influence of confining pressure and lateral pressure on mechanical properties [25].
Under the axial compression in cylinder, the stress-strain curve of the backfill materials presents typical nonlinear characteristics, which can be accurately fitted by polynomial negative exponential functions [26].
Due to the resistance of the filling body to overburden, the mining strength of solid backfill mining is small. e finite element method based on continuum mechanics can be used to numerically simulate the movement and deformation of overburden rock, so as to evaluate the influence of different backfill mining schemes on key strata and surface subsidence. rough the numerical simulation, Huang Yanli et al. analyzed the control effect of different compaction rates of the filling body on the strata movement and surface subsidence during the solid backfill mining [27]. Zha Jianfeng et al. used numerical simulation method to study the influence of gangue backfill as double yield model and elastic model on strata movement. However, the backfill materials are considered as elastic materials in the abovementioned studies, and the compaction process of gangue filling into goaf under the action of overburden cannot be correctly reflected [28]. Li Meng proposed a backfill material model considering pressure coefficient and performed the numerical simulation of overburden movement and deformation characteristics based on FLAC3D software [29]. However, in this constitutive model of backfill materials, the influence of lateral stress was not considered in the axial deformation modulus. Based on the negative exponential function curve of backfill material, Zhang et al. proposed a negative exponential function model of backfill mining and predicted the strain and deformation law of overburden during backfill mining by the numerical simulation in ABAQUS software [30]. However, the effect of lateral stress in the compaction experiment was considered in this model, and the compaction effect of the compaction machine on backfill materials after backfill mining was not considered in the simulation. As a result, the mechanical behavior of backfill materials cannot be accurately described.
According to the continuum theory, experimental results, and mechanical characteristics of solid backfill materials, a constitutive model of backfill materials was proposed based on the relationship between volume strain and deformation modulus. e platform of ABAQUS UMAT was used, then the numerical development of the model was performed, and the influence of the initial compaction rate of backfill materials, mining height, and mining depth on the movement and deformation law of overburden was systematically analyzed.

Mechanical Properties of Solid
Backfill Materials e physical and mechanical properties of solid backfill materials are the important basis for the analysis of overburden movement in solid backfill mining. Solid backfill materials are generally composed of broken coal gangue, fly ash, loess, and other industrial wastes in the form of granular materials. e mechanical property test of solid backfill material was carried out with the compaction experimental device, as shown in Figure 1. During the experiment, the backfill materials were placed inside the cylinder, and then the piston was pressed by the mechanical testing machine to compact the backfill materials. During the compaction process, the compression force and displacement of the machine were collected, and the loading rate was quasilinear loading. After the experiment, the axial stress and strain of backfill materials were calculated through the pressure and displacement of the testing machine. Figure 2 shows the stress-strain curves of the four kinds of broken rock during the compaction process. It can be seen that the stress-strain curve of the backfill materials changes in a nonlinear way during the compaction process. At the initial compressive deformation, the resistance of backfill material is low because there are numerous of gaps in the backfill materials. e stress increases slightly at the initial compressive deformation but its growing rate rises gradually at subsequent deformation. With the increase of strain, the increasing range of stress gradually increases. By using regression analysis, some scholars proposed that the stress-strain relationship of broken rock is an exponential function [31]: where σ and ε are the axial stress and strain of fractured rock mass; A and B are regression coefficients.

Constitutive Model of Solid Backfill
Material. e backfill materials of solid backfill mining are relatively loose, and there are a lot of voids in the material. erefore, the solid backfill material is not an ideal elastic body, but a porous discontinuous medium with elasticity, plasticity, and viscosity. When the finite element method is used to calculate the surface subsidence law after backfill mining, the stressstrain relationship of backfill materials during the loading process should be established to effectively reveal the resistance of backfill materials to the deformation of overburden. Firstly, based on the physical and mechanical properties of backfill materials during the backfill mining, the following assumptions are put forward for backfill materials: (1) Each component in the backfill materials is evenly distributed, and the backfill materials are isotropic materials. Based on the abovementioned hypothesis, the constitutive model of backfill materials is the constitutive relation of elastic materials: where σ, ε, and D are the stress tensor, strain tensor, and fourth-order elastic tensor of backfill materials, respectively. Besides, D is related to deformation modulus E and Poisson's ratio μ. As the stress and deformation of backfill materials are nonlinear in the process of compaction, it is necessary to find out the relationship between deformation modulus E, Poisson's ratio μ, and strain, so as to further determine the constitutive model of backfill materials. e Poisson's ratio μ is difficult to be determined because the backfill material is loose, so Poisson's ratio μ is assumed to be a constant. According to the generalized Hooke law, the stress-strain relationship of backfill material in the compaction experiment is as follows: where σ z and ε z are the stress and strain in the compression direction, respectively; εv is the volume strain, and ε v � ε x + ε y + ε z . In the compaction process, because the rigidity of the cylinder is large, ε x and ε y are far less than ε z . erefore, ε v ≈ ε z , and the above formula can be further simplified as follows: where m is a constant related to Poisson's ratio. Combined with equation (1), equation (5) can be obtained: By introducing equation (5) into equation (1), the constitutive relation of backfill materials can be obtained.

Numerical Simulation.
Based on ABAQUS UMAT, the subroutine of solid backfill materials is developed. e calculation flow ( Figure 3) is as follows: (1) At the beginning of the incremental step, the basic parameters of the material are given, including parameters A, B, and μ and initial volume strain ε vini . After backfill materials are filled into the goaf, the compaction machine of filling hydraulic support tamps the filling body to ensure that the backfill material has a certain initial compaction rate. e compaction rate is the ratio of the compacted volume to the initial volume, and the relationship between the compaction rate and the volume strain is as follows:

Shock and Vibration
where V ys is the volume of backfill materials after compaction, V s is the initial volume, and K y sin i is the initial compaction rate. (2) Because ABAQUS calls the same subroutine every time, the volume strain of the model is different in the first load step and the subsequent load step. en, the volume strain should be judged before starting the program. If ε x + ε y + ε z � 0, the volume strain is the initial volume strain ε ini ; if ε x + ε y + ε z ≠ 0, the volume strain is the volume strain stored in the previous loading step for updating the deformation modulus. (3) e deformation modulus is updated by equation (5), and then the Jacobian matrix is updated, and the current volume strain is stored. (4) After the end of the incremental step, ABAQUS main program can automatically perform the convergence judgment. If it converges, the next incremental step is entered; if not, the incremental step is reduced until it converges.

Model Validation.
e compaction process of backfill material is spatially axisymmetric, and the central axis of the cylinder is the axis of symmetry. erefore, the finite element model is established by using the axisymmetric element in ABAQUS, as shown in Figure 4. In the experiment, the inner wall of the cylinder is coated with lubricating oil. erefore, friction-free contact is established between the cylinder and the filling body and between the piston and the filling body to simulate the interaction between the backfill materials and the cylinder. e bottom of the model restricts the vertical displacement and limits the horizontal displacement of the symmetry axis. e loading mode is displacement loading with the size of U z � 60.75mm, i.e., ε z � 0.45. e elastic modulus and Poisson's ratio of steel used for cylinder and piston are 220 GPa and 0.3, respectively. e parameters of backfill materials are shown in Table 1. Figures 5 and 6, respectively, show the stress-strain curves and elastic modulus-volume strain curves of fractured shale, mudstone, coal, and sandstone. It can be seen that the experimental values and numerical simulation values of the stress-strain curves and elastic modulus-volume strain curves of the four broken rock masses are consistent. It indicates that the backfill material model developed by UMAT is reliable and can be UMAT start Equation (5)

Numerical Calculation Model of Overburden Movement and Deformation Law in Solid
Backfill Mining e engineering geological conditions of the third mining area in No. 12 coal mine of Pingdingshan Coal Group Co., Ltd., was taken as the background, and the large-scale general finite element calculation software ABAQUS was used as the analysis and calculation platform. According to the geological conditions and the symmetry characteristics of the working face, the numerical mechanical model of the plane strain was established, as shown in Figure 7. e model size is 1000 m × 380 m; the horizontal displacement is constrained on both sides of the model, and the vertical displacement is constrained at the bottom. e rock layer exerts a body force of ρg, where ρ is the density and g is the acceleration of gravity.
e Mohr-Coulomb model of ABAQUS is used for the parameters of each stratum, and the physical and mechanical parameters are shown in Table 2. Considering the factors of calculation accuracy and calculation time, the mesh near the coal seam is refined. e model is divided into 4085 grids, as shown in Figure 7. e mechanical properties of the filling body play an important role in controlling the overburden movement of the working face. According to the stress-strain test results of the solid backfill material, broken rock #4 is selected as the backfill material, and the material parameters are shown in Table 2.
To study the law of overburden deformation and surface subsidence under the condition of backfill mining, the influence of the initial compaction rate of the filling body, mining height, and mining depth on overburden deformation and surface subsidence is mainly considered in the model; the numerical calculation scheme is shown in Table 3.  Shock and Vibration the right of the middle of the working face, and the vertical dimension is 380 m from the bottom of the model to the surface. It can be seen that with the decrease of the initial compaction rate, the small displacement contour (0.05) of overburden gradually moves to the surface, and the scope of overburden between adjacent isolines gradually increases. Figure 9(a) shows the subsidence curve of the main key strata. When the initial compaction rate of filling decreases from 0.95 to 0.85, the deformation of the key strata of overburden gradually decreases, and the reduction amplitude is significant. Figures 9(b) and 9(c) show the surface subsidence curve and the curve of the maximum surface subsidence varying with the initial compaction rate of the filling body. It can be seen that, with the increase of the initial compaction rate, the surface subsidence decreases obviously. Similar to the deformation law of the key strata of overburden, when the initial compaction ratio of filling decreases from 0.95 to 0.85, the surface subsidence gradually decreases, and the reduction amplitude is obvious. As shown in Figures 9(a) and 9(b), the surface subsidence curve is similar to that of the main key strata, but the surface subsidence value under the same compaction rate of the filling body is smaller than that of the main key strata, because the stress redistribution occurs in the overburden after mining, and the main key strata are close to the stope, which will cause great disturbance after mining. In addition, the main key strata and surface subsidence have a nonlinear change law with the initial compaction rate, because the initial compaction rate reflects the initial volume strain of the backfill material, which has a nonlinear relationship with the deformation modulus of the backfill material.

Results and Discussion
ere is a complex influence of mining depth on the surface subsidence and overburden deformation, and the structural composition and lithologic characteristics of    Shock and Vibration 600 m, the w kmax of the key strata of overburden increases from 366.9 mm to 455.5 mm, with an increase of 24.1%. is increase range is obviously smaller than that of mining depth (71.4%). In fact, with the increase of mining depth, the structure and lithologic composition of overburden have a great influence on the deformation of key strata of overburden, which is the main factor of overburden deformation resistance. Figures 11(b)   subsidence basin L increases from 1184 m to 1970 m, with an increase of 66.4%. It indicates that mining depth has a great influence on the surface subsidence range. Figure 12 shows the contour map of vertical displacement of overburden at the different mining heights of the filling bod (2.5-5.0 m). In this figure, the horizontal dimension is 0 m-350 m, that is, 350 m from the middle of the working face to the right, and the vertical dimension is 380 m, that is, 380 m from the bottom of the model to the surface. Figure 13(a) shows the subsidence curve of the main key strata at different mining heights and the variation curve of the subsidence peak value of the main key strata with the mining height. It can be seen that with the increase of mining height, the contour lines of the same size gradually move to the upper part of the stope, and the deformation of the key strata of the overburden increases gradually. When the mining height h increases from 2.5 m to 5.0 m, the subsidence of the key strata of overburden increases from 267.4 mm to 530.2 mm. It indicates that the mining height has a significant influence on the deformation of the key strata of the overburden. Figure 13( Figure 9: (a) Subsidence curve of main key strata, (b) surface subsidence curve, and (c) variation curve of maximum subsidence under different initial compaction rates. and the variation curve of the peak value w max of surface subsidence with the mining height h. When the initial compaction rate of the filling body is 0.8, the w max increases with the increase of mining height. Figure 13(c) shows the variation curves of the maximum subsidence of the surface and the main key strata with the mining height, both of which show a linear change law with the mining height.
When the mining height h increases from 2.5 m to 6.0 m, the w max increases from 218.9 mm to 434.0 mm, indicating that the mining height has a significant impact on the surface subsidence in solid backfill mining. e surface subsidence is similar to the deformation law of the key strata of overburden, but its subsidence value is less than that of the main key strata.

Conclusion
(1) In this paper, the constitutive model of solid backfill materials based on elastic constitutive relation of solid material is proposed and used to accurately reveal the compaction mechanical characteristics of solid backfill materials. e ABAQUS UMAT subroutine, which is developed by the constitutive model, can be used to effectively predict the movement and deformation law of overburden in solid backfill mining; (2) In solid backfill mining, the initial compaction rate and the height of coal seams have a significant impact on the surface subsidence. When the thickness of the topsoil is only changed and the structural composition and lithology of the overburden are not changed, the mining depth has little impact on the surface subsidence, but has a significant impact on the surface subsidence range; (3) In solid backfill mining, the influence of mining height and mining depth on the deformation of key strata in the overburden and surface subsidence is approximately linear, but the influence of the initial compaction rate is nonlinear.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.