Abstract
With the gradual increase in mining depth of coal resource exploitation, deep backfilling mining has effectively solved the impact of strong deep mine pressure and strong mining disturbances. However, after deep backfilling mining, the backfilling material is subjected to high stress for a long time, and its viscoelasticity has a significant impact on the roof control effect. This paper uses a large-scale bulk confinement test device to analyze the viscoelastic properties of gangue, establishes a high-precision fractional viscoelastic creep model, and identifies the gangue parameters. The established fractional viscoelastic model was used as the foundation model of the beam, and the roof model based on the fractional viscoelastic foundation was solved. The top deformation characteristics of elastic foundation and fractional foundation were compared and analyzed, and the time effect, viscoelastic effect, and order effect of the fractional order viscoelastic foundation beam were discussed. The results show that the viscosity of gangue increased under the action of deep high stress. As time increased, the roof deformation also increased. In order to more effectively control the long-term deformation of the roof, the viscosity coefficient of the backfilling material should be greater than 20 MPa. This research provides relevant guidance for the requirements of backfilling materials for deep backfilling mining and the prediction of long-term dynamic deformation of the roof in underground excavations.
1. Introduction
As the mining depth increases for coal resources exploitation, deep backfilling mining technology has achieved considerable development in China and has achieved significant economic, social, and environmental benefits in many mining areas [1]. Due to the fundamental changes of the environment in deep mining, the surrounding rock environment of the gangue backfill body will face the conditions of high ground stress, high water pressure, and strong mining disturbances, which will inevitably affect the internal void structure and mechanical properties of gangue and long-term deformation characteristics [2, 3]. Under the above-mentioned complicated deep environmental conditions, the mechanical properties of the rock are significantly different from those under the shallow in situ stress state. The specific manifestation is that the rock mass under high stress in the deep has continuous strong rheological properties [4]. In addition, deep rock mass is prone to brittle-ductile transformation, and its failure shows strong plastic deformation characteristics [5]. The gangue backfill body is used as a medium to prevent the deformation of the rock formation [6]. The mechanical properties, especially the long-term rheological and viscoplastic properties under the above-mentioned high stress-water environmental conditions, are very important for the roof control and backfilling effect.
Scholars worldwide have carried out many researches on backfilling material properties and solid backfilling mining. Zhang et al. [7] studied the distribution and characteristics of granular materials, cement materials and higher water content materials, and a guide for the selection of backfilling materials. The finite compression experiment was conducted to analyze the deformation and crushing characteristics of the gangue bulk materials under different water immersion heights by Li et al. [8]. Li et al. [9] tested compaction properties of crushed gangue backfilling materials with different particle sizes and also analyzed the influences of particle size of the backfilling materials on surface subsidence. Huang et al. [10] analyzed the mesostructure, stress variation, energy dissipation, and backfilling effects of the five common solid backfilling materials. In addition, the time-dependent characteristics of gangue and other backfilling materials are also a hot subject of research. Researches showed that gangue creep has viscoelastic properties [11–15]. Moreover, the viscoelastic properties of gangue also have an important impact on the deformation of the roof of the backfilled space in backfilling mining [16–19]. However, many roof models of solid backfilling mining are based on elastic foundation assumptions, ignoring the viscoelasticity of the backfilling materials [20–24].
Fractional calculus is a branch of mathematics. As an excellent mathematical tool, it has been successfully applied in fields such as rheology, electromagnetics, signal processing, and chaotic systems [25]. Abel dashpot based on fractional derivative, showing the characteristics of Spring and Newtonian dashpot, with viscoelastic characteristics. Some classic viscoelastic models, such as the Maxwell model, Kelvin-Voigt model, and general Kelvin-Voigt model, are improved by using the fractional Abel dashpot [26–28]. Those fractional viscoelastic models have much better advantages for describing the material’s rheological mechanical properties than the integer-order calculus [29].
This paper uses a large-scale confined compression test device to study the creep characteristics of gangue and establishes a high-precision fractional viscoelastic creep model of gangue. The fractional viscoelastic creep model is used as the foundation of the beam, and the fractional viscoelastic foundation roof model is also established. The deformation, viscoelastic effect, and time effect of the backfilling mining roof are analyzed. This research can be used to predict the dynamic subsidence of the deep backfilling mining roof and provides relevant theoretical guidance for rock strata long-term deformation.
2. Gangue Creep Test
2.1. Sample and Testing Procedure
The particle size range of gangue is 0-30 mm, from the Pingdingshan coal mines, which matches the actual gradation of gangue in the goaf. The main lithology of gangue materials is mudstone with an average density of 1450 kg·m-3 and a moisture content of 0.42%. The X-ray diffraction (XRD) results show that mudstone mainly consists of 37.36% quartz, 34.23% kaolinite, 11.22% illite, 2.11% chlorite, 1.15% feldspar, and other minerals. The YAS-5000 electrohydraulic servocontrolled rock mechanic test system was employed for the backfilling material compaction tests. The compacted steel barrel is a cylindrical compacted steel barrel designed by the authors, including the loading platen, handle, bolt and base, as shown in Figure 1. The steel barrel’s inner diameter is 250 mm, height is 305 mm, and yield limit is 235 MPa. The uniaxial creep tests were performed on gangue materials by using the steel circular cylinder, for constant loads of 5 MPa, 10 MPa, 15 MPs, and 20 MPa stress in the axial direction for 60 mins each, and with a loading rate of 1 kN/s.
(a) Compacted steel barrel
(b) Dimensions
2.2. Testing Result
The relationships between strain and time, deformation rate, and time of gangue under four different axial pressures are shown in Figure 2.
(a) Relationships between strain and time
(b) Relationships between deformation rate and time
As shown in Figure 2, the deformation of gangue is not linear with time in lateral confined creep tests. With the increase of axial pressure, the creep strain of gangue reduces gradually, such as Δε1, Δε2, and Δε3. The deformation of gangue can be divided into two stages: deceleration deformation and stable deformation, which clearly defines the viscoelastic properties of gangue. Therefore, when predicting the long-term deformation of gangue, the viscoelastic characteristics should not be ignored.
3. Gangue Material Creep Model and Parameter Identification
In order to study the viscoelastic characteristics of the gangue creep process, it is necessary to establish a high-precision creep model. The fractional calculus theory has obvious advantages in the accuracy and prediction of the creep model [30–32].
3.1. Fractional Viscoelastic Creep Model
3.1.1. Basic Concept of Fractional Calculus
The th derivative of a function is expressed as . If α is a fraction, this expression becomes a fractional derivative. Fractional calculus is usually expressed in terms of the Riemann-Liouville definition. The Riemann-Liouville fractional integration of function of order α is defined as [33]
where 0 and are the limits of integration at the left and right sides of ; is the gamma function with argument α, which is defined by
The fractional derivative is defined by
where , , and are positive integers.
3.1.2. Abel Dashpot
AD is a component derived from the integer-order creep element, as shown in Figure 2.
AD’s stress and strain have the following expressions:
where is the fractional order; is the time, h; is the coefficient of viscosity; is the stress, MPa; is the constant stress, MPa; and is the strain.
According to formula (5), given the relevant mechanical parameters of the material, the creep characteristic curves under different orders are obtained, as shown in Figure 3.
It can be seen from Figure 3 that under the condition of constant stress, as α increases, the corresponding strain rate increases; when α is a fractional order, the strain curve increases slowly, rather than linearly like Newtonian fluid increase. However, the elastic body strain remains unchanged; the creep of the material shows a nonlinear gradual process, indicating that the properties of the material are between fluid and solid; and the fractional viscous pot can well reflect the viscoelastic properties of the material.
3.1.3. Fractional Poynting-Thomson Model
According to the viscoelastic relationship of the gangue creep experiment, the fractional Poynting-Thomson model is selected as the viscoelastic foundation model of gangue. The model consists of a fractional Maxwell model with elastic elements in parallel, as shown in Figure 4.
The constitutive relationship of each element in the fractional Maxwell model is
According to the series relationship,
Laplace transformation is carried out on formula (6) and substituted into formula (7):
Carrying out the Laplace inverse transformation on formula (8), the constitutive relationship of Maxwell model can be obtained as
When , Laplace transformation and inverse transformation were performed on equation (9), and the creep relationship of Maxwell model can be obtained as
The constitutive relationship of each element in the fractional Poynting-Thomson model is
From the parallel relationship and the Laplace transformation of equation (11),
Carrying out the Laplace inverse transformation on formula (12), the constitutive relationship of the fractional Poynting-Thomson model can be obtained as
When , Laplace transformation was performed on equation (13) to obtain
After proving the formula (14), formula (15) was obtained:
The Laplace variation expression of the two-parameter Mittag-Leffler function is
Carrying out Laplace inverse transformation on equation (15) from equation (16), the creep relationship of the fractional Poynting-Thomson model can be obtained as
3.2. Fractional Poynting-Thomson Model Parameter Identification
The expression of generalized Hookes law in three-dimensional form is as follows: where is the first invariant of the stress tensor, MPa; is the first invariant of strain tensor; is the deviatoric stress tensor, MPa; is the deviatoric strain tensor; is the bulk modulus, GPa; and is the shear modulus, GPa.
Usually, when the rock is in the three-dimensional stress state, the spherical stress tensor will cause the volume deformation of the element without changing its shape, while the deviatoric stress tensor causes element shape change without volume change. The expressions of stress state and strain state are as follows: where is the spherical stress tensor in MPa and is the spherical strain tensor.
Therefore, the strain expression of elastomer is
In the elastic state, assuming that the spherical strain tensor and deviatoric strain tensor are related to the deviatoric stress tensor, respectively, the three-dimensional expression of the fractional Poynting-Thomson model is
with the confining pressure and axial pressure having the following relationship:
Therefore, the three-dimensional creep equation of gangue under constant pressure in the steel drum is expressed as follows:
There are usually two methods to determine the parameters of the creep model. One is the graphical method, according to the relationship between the creep curve and the physical meaning [34], and another method is the optimization analysis method, such as the nonlinear least square method, which is widely used [35, 36]. In this paper, the widely used Levenberg–Marquardt nonlinear least square method is adopted for parameter identification, as shown in Table 1.
In order to verify the validity of the creep model parameters, the comparison between the creep test data of gangue and the model curve is shown in Figure 5. It can be seen from the figure that the fractional Poynting-Thomson model is in good agreement with the experimental data, and the accuracy and are above 0.99.
4. Roof Mechanics Model of Backfilling Mining Space Based on Fractional Viscoelastic Foundation
4.1. Mechanical Model of Fractional Viscoelastic Foundation Beam
The basic roof in solid backfilling coal mining can be regarded as a fixed beam on both sides supported by the upper part of the overlying rock load and the lower part supported by the foundation reaction force , and the deflection of the beam section is expressed as as shown in Figure 6. The working strike direction is established to be the positive direction, the direction of inclination to be the positive direction, and the vertical direction to be the positive direction.
4.2. Solution of Integer Order Solution of Elastic Foundation
Suppose the basic differential equation of elastic foundation beam deflection curve is
When the foundation is an elastic foundation, the expression of is given by where is the elastic foundation coefficient which can be obtained by the following formula: is the elastic modulus of the cushion, GPa; is the thickness of the cushion, m.
The relationship among the corner , the bending moment , the shear force , and the deflection of the beam section is
Taking characteristic coefficients as , when the gangue backfilling material is regarded as an elastic material, the deflection curve equation of the roof beam can be obtained by solving equation (24):
where the , , , and are the parameters.
4.3. Solution of Fractional Viscoelastic Foundation
From the constitutive relationship of the fractional Poynting-Thomson model (13), the relationship between the stress and deformation of the viscoelastic foundation is given by
Laplace transformation of equation (24) can be expressed as
Laplace transformation of equation (13) can obtain the Laplace relationship between foundation reaction force and deformation as follows:
From equation (29) and equation (31), the corresponding relationship between elastic foundation and viscoelastic foundation is derived and expressed as
From equation (28), when the foundation is an elastic material, the Laplace expression of the bending deflection equation of the thin plate is
Substituting formula (32) into (33), the Laplace deflection equation of the viscoelastic foundation sheet becomes
By arranging the formula (34) to its final form, we can get
Separating the variables from equation (35) and extracting the common factor, we get
Carrying out the Laplace inverse transformation of equation (36) and introducing the two-parameter Mittag-Leffler function, the deflection equation of the viscoelastic foundation sheet can be obtained as
4.4. Case Study
According to the geological conditions of the example coal mine panel, it has a length of 280 m and a width of 90 m. The average mining height of the coal seam is 2.7 m. The roof of the coal seam is sandstone with 4.5 m thickness, elastic modulus of 5.8 GPa, and Poisson’s ratio of 0.27. The load above the roof mainly comes from the weight of potential caving strata a, . The creep time is 24 hours. Considering the calculation time and accuracy, the first three terms of infinite series and in equation (37) and the first ten terms of infinite series in double-parameter Mittag-Leffler function (16) are calculated. The boundary conditions of the beam are as follows:
By substituting boundary conditions (formulas in (38)) and the above parameters and creep parameters (from Table 1) into formulas (28) and (37), the parameters are as shown in Table 2.
The roof subsidence of elastic foundation and fractional foundation can be obtained as shown in Figure 7.
It can be seen from Figure 7 that the elastic foundation and the fractional viscoelastic foundation beam have the same deformation trend, and the maximum subsidence is the same. However, the local deformation of the fractional viscoelastic foundation is smaller, just because the fractional viscoelastic foundation has the higher global correlation. The roof subsidence calculated by the classic elastic foundation beam has been confirmed by many studies. It can be seen from the calculation of the final subsidence that the fractional viscoelastic foundation beam model is also reliable.
5. Discussion
5.1. Time Characteristics of Fractional Viscoelastic Foundation
Referring to the above for other parameters, the roof subsidence with different times can be obtained as shown in Figure 8. And the relationship of maximum roof subsidence and time is shown in Figure 9.
It can be seen from Figures 7 and 8 that as time increases, the roof deformation also increases accordingly. However, the increase in roof deformation gradually becomes smaller. The amount of roof subsidence increased from 154 mm in 7 d to 304 mm in 10 years, which shows an increase of 49.3%. Compared with the elastic foundation and fractional viscoelastic foundation (10 years later), the roof deformation increased by 56.6%.
5.2. Viscosity Characteristics of Fractional Viscoelastic Foundation
The roof subsidence with different viscosity coefficients can be obtained as shown in Figure 10. And the relationship of maximum roof subsidence and viscosity coefficient is shown in Figure 11.
It can be seen from Figures 10 and 11 that as the viscosity coefficient increases, the roof deformation decreases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, which is smaller than . The amount of roof subsidence reduced from 637 mm with a viscosity coefficient with 5 MPa to 43.1 mm with a viscosity coefficient with 100 MPa, a proportional decrease of 93.2%. Compared with the elastic foundation and fractional viscoelastic foundation (10 years later), the roof deformation increased by 56.6% as in Section 5.1. It is also worth noting that the roof subsidence decreases sharply with the increase in coefficient of viscosity. When the viscosity coefficient is greater than 20 MPa-min, the effect of viscosity of the backfilling material on the roof control is limited. In other words, the viscosity coefficient of the backfilling material should reach 20 MPa-min in order to better control the sinking of the roof.
5.3. Fractional Order Characteristics of Fractional Viscoelastic Foundation
The roof subsidence with different fractional orders can be obtained as shown in Figure 12. And the relationship of maximum roof subsidence and fractional order is shown in Figure 13.
It can be seen from Figure 12 that as the fractional order increases, the roof deformation increases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, is smaller than . It can be also seen from Figure 13 that the relationship of maximum roof subsidence and fractional order is a nonlinear relationship. The fractional order can be obtained by fitting the experimental data, and the calculated roof subsidence is more accurate due to the fact that the fractional order is variable (moves in a changing trend).
6. Conclusions
The viscoelastic characteristics of gangue during creep deformation were analyzed. A fractional order viscoelastic foundation beam model was established; the differential equation of the fractional order viscoelastic foundation beam was obtained by Laplace transform; the deflection of classical elastic foundation and the fractional order viscoelastic foundation beam were compared and analyzed; and the time, viscosity, and order pair were discussed. The drawn conclusions are the following: (1)The viscoelastic properties of gangue were analyzed through the graded creep test experiment. The deformation of gangue was divided into two stages of deceleration deformation and stable deformation. A high-precision fractional Poynting-Thomson creep model was established, and the parameters of gangue fractional Poynting-Thomson creep model were identified(2)The viscoelastic foundation beam model was established by the gangue fractional Poynting-Thomson creep model, the deflection equation of the roof was obtained by Laplace transform, and the parameters of the viscoelastic foundation beam were obtained according to the boundary conditions. Compared with the classical elastic foundation beam model, the fractional viscoelastic foundation beam result is the most reliable(3)According to the established viscoelastic foundation beam, the influence of time, viscosity, and order parameters on roof deformation was discussed. The fractional viscoelastic foundation beam can predict roof deformation according to time; appropriately increasing the viscosity of the backfilling material can be effective in the control of roof subsidence; and the nonlinear fractional order can improve the accuracy of the creep model, thereby improving the accuracy of the solution for the viscoelastic foundation beam
Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was funded by the National Key R&D Program of China (grant number 2018YFC0604704), the National Science Fund for Distinguished Young Scholars (grant number 51725403), the China Postdoctoral Science Foundation (grant number 2019M661990), and the Fundamental Research Funds for the Central Universities (grant number 2020QN05).