Nonlinear Creep Model of Deep Gangue Backfilling Material and Time-Dependent Characteristics of Roof Deformation in Backfilling Mining

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 20MPa. 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.


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 manifesta-tion 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][12][13][14][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][17][18][19]. However, many roof models of solid backfilling mining are based on elastic foundation assumptions, ignoring the viscoelasticity of the backfilling materials [20][21][22][23][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][27][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 longterm deformation.

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.

Testing
Result. The relationships between strain and time, deformation rate, and time of gangue under four different axial pressures are shown in Figure 2.
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.

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 highprecision creep model. The fractional calculus theory has obvious advantages in the accuracy and prediction of the creep model [30][31][32].
where 0 and t are the limits of integration at the left and right sides of D; ΓðαÞ is the gamma function with argument α, which is defined by The fractional derivative is defined by where α ≥ 0, n-1 ≤ α ≤ n, and n are positive integers.

Abel
Dashpot. AD is a component derived from the integer-order creep element, as shown in Figure 2.
2 Geofluids AD's stress and strain have the following expressions: where α is the fractional order; t is the time, h; η c is the coefficient of viscosity; σ is the stress, MPa; σ c 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 grad-ual 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.  Figure 4.
The constitutive relationship of each element in the fractional Maxwell model is According to the series relationship,

Geofluids
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 σ 1 ðtÞ = σ 0 , 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 σ 1 ðtÞ = σ 0 , 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

Fractional Poynting-Thomson Model Parameter
Identification. The expression of generalized Hookes law in three-dimensional form is as follows: where σ m is the first invariant of the stress tensor, MPa; ε m is the first invariant of strain tensor; S ij is the deviatoric stress  4 Geofluids tensor, MPa; e ij is the deviatoric strain tensor; K is the bulk modulus, GPa; and G 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 σ m δ ij is the spherical stress tensor in MPa and ε m δ ij 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 R 2 are above 0.99.

Roof Mechanics Model of Backfilling Mining
Space Based on Fractional Viscoelastic Foundation

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 qðxÞ and the lower part supported by the foundation reaction force pðxÞ, and the deflection of the beam section is expressed as wðxÞ as shown in Figure 6. The working strike direction is established to be the positive x direction, the direction of inclination to be the positive y direction, and the vertical direction to be the positive z direction.

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 pðxÞ is given by where k 0 is the elastic foundation coefficient which can be obtained by the following formula: E 0 is the elastic modulus of the cushion, GPa; h 0 is the thickness of the cushion, m. The relationship among the corner θ, the bending moment M, the shear force Q, and the deflection wðxÞ of the beam section is , 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 d 1 , d 2 , d 3 , and d 4 are the parameters.

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 w x, s ð Þ= e α 0 x d 1 sin

Geofluids
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 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.

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%.

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, Δw 1 which is smaller than Δw 2 . 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.

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, Δw 1 is smaller than Δw 2 . 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).

Conclusions
The viscoelastic characteristics of gangue during creep deformation were analyzed. A fractional order viscoelastic 8 Geofluids 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.