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

Authors: Peng Huang; Jixiong Zhang; Qi Zhang; Ntigurirwa Jean Damascene; Yuming Guo
Journal: Geofluids
 (2020)

---

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


---

## Body

## 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.Figure 1
Compacted steel barrel and dimensions.
(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 Figure2.Figure 2
Gangue creep properties.
(a)
Relationships between strain and time(b)
Relationships between deformation rate and timeAs shown in Figure2, 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.


## 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.Figure 1
Compacted steel barrel and dimensions.
(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 Figure2.Figure 2
Gangue creep properties.
(a)
Relationships between strain and time(b)
Relationships between deformation rate and timeAs shown in Figure2, 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

Thenth derivative of a function ft is expressed as Dαft=dαft/dtα. 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 ft of order α is defined as [33]
(1)Dt−αft=1Γα∫0tt−ξα−1fξdξReα>0,t>0,where 0 andt are the limits of integration at the left and right sides of D; Γα is the gamma function with argument α, which is defined by
(2)Γα=∫0∞e−ttα−1dtReα>0.The fractional derivative is defined by(3)Dαft=dnD−n−αftdtn,whereα≥0, n‐1≤α≤n, and n are positive integers.


#### 3.1.2. Abel Dashpot

AD is a component derived from the integer-order creep element, as shown in Figure2.AD’s stress and strain have the following expressions:(4)σt=ηcdαεtdtα0≤α≤1,α∈R+,t>0,(5)εt=σcηtαΓα+10≤α≤1,α∈R+,t>0,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.Figure 3
Effect of different orders ofα on creep properties.It can be seen from Figure3 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 Figure4.Figure 4
Fractional Poynting-Thomson model.The constitutive relationship of each element in the fractional Maxwell model is(6)σ1=E1ε1t,σ1=ηαDαε2t.According to the series relationship,(7)εt=ε1t+ε2t.Laplace transformation is carried out on formula (6) and substituted into formula (7):
(8)E1ηαsαε¯s=ηαsασ1¯s+E1σ1¯s.Carrying out the Laplace inverse transformation on formula (8), the constitutive relationship of Maxwell model can be obtained as
(9)E1ηαDαεt=ηαDασ1t+E1σ1t.Whenσ1t=σ0, Laplace transformation and inverse transformation were performed on equation (9), and the creep relationship of Maxwell model can be obtained as
(10)εt=σ0E1+σ0ηαtαΓα+1.The constitutive relationship of each element in the fractional Poynting-Thomson model is(11)Dασ1t+E1ηασ1t=E1Dαε1t,σ2t=E2ε2t.From the parallel relationship and the Laplace transformation of equation (11),
(12)ηαsασ¯s+E1σ¯s=E1ηαsαε¯s+E1E2ε¯s+E2ηαsαε¯s.Carrying out the Laplace inverse transformation on formula (12), the constitutive relationship of the fractional Poynting-Thomson model can be obtained as
(13)σt+ηαE1Dασt=E2εt+ηα+E2ηαE1Dαεt.Whenσ1t=σ0, Laplace transformation was performed on equation (13) to obtain
(14)ε¯s=E1σ0E1E2+E1+E2ηαsαs+ηασ0E1E2s1−α+ηαE1+E2s.After proving the formula (14), formula (15) was obtained:
(15)ε¯t=E1σ0E1+E2−1η−αs−1E1E2/E1+E2ηα+sα+σ0E1+E2−1sα−1E1E2/ηαE1+E2+sα.The Laplace variation expression of the two-parameter Mittag-Leffler function is(16)∫0∞e−sttαk+β−1Eα,βk±atαdt=k!sα−βsα∓ak+1Rep>a1/α.Carrying out Laplace inverse transformation on equation (15) from equation (16), the creep relationship of the fractional Poynting-Thomson model can be obtained as
(17)εt=E1σ0E1+E2ηαtαEα,α+1−E1E2E1+E2ηαtα+σ0E1+E2Eα,1−E1E2E1+E2ηαtα.


### 3.2. Fractional Poynting-Thomson Model Parameter Identification

The expression of generalized Hookes law in three-dimensional form is as follows:(18)εm=σm3K,eij=Sij2G,where σm is the first invariant of the stress tensor, MPa; εm is the first invariant of strain tensor; Sij is the deviatoric stress tensor, MPa; eij 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:(19)σij=σmδij+Sij,εij=εmδij+eij,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(20)εije=Sij2G+σmδij3K.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(21)εt=G1Sij2G1+G2ηαtαEα,α+1−G1G2G1+G2ηαtα+Sij2G1+G2Eα,1−G1G2G1+G2ηαtα,with the confining pressure and axial pressure having the following relationship:(22)σ1=λσ2=λσ3.Therefore, the three-dimensional creep equation of gangue under constant pressure in the steel drum is expressed as follows:(23)εt=G1σ11−λ3G1+G2ηαtαEα,α+1−G1G2G1+G2ηαtα+σ11−λ3G1+G2Eα,1−G1G2G1+G2ηαtα.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.Table 1
Parameters of fractional Poynting-Thomson model.
σ1 (MPa)G1 (MPa)G2 (MPa)η (MPa·min)αR25.03.10.5644.670.210.9910.03.970.7498.960.190.9915.04.861.1368.210.110.9920.05.751.25100.000.110.99In 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 Figure5. 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 R2 are above 0.99.Figure 5
Fractional Poynting-Thomson model.


## 3.1. Fractional Viscoelastic Creep Model

### 3.1.1. Basic Concept of Fractional Calculus

Thenth derivative of a function ft is expressed as Dαft=dαft/dtα. 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 ft of order α is defined as [33]
(1)Dt−αft=1Γα∫0tt−ξα−1fξdξReα>0,t>0,where 0 andt are the limits of integration at the left and right sides of D; Γα is the gamma function with argument α, which is defined by
(2)Γα=∫0∞e−ttα−1dtReα>0.The fractional derivative is defined by(3)Dαft=dnD−n−αftdtn,whereα≥0, n‐1≤α≤n, and n are positive integers.


### 3.1.2. Abel Dashpot

AD is a component derived from the integer-order creep element, as shown in Figure2.AD’s stress and strain have the following expressions:(4)σt=ηcdαεtdtα0≤α≤1,α∈R+,t>0,(5)εt=σcηtαΓα+10≤α≤1,α∈R+,t>0,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.Figure 3
Effect of different orders ofα on creep properties.It can be seen from Figure3 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 Figure4.Figure 4
Fractional Poynting-Thomson model.The constitutive relationship of each element in the fractional Maxwell model is(6)σ1=E1ε1t,σ1=ηαDαε2t.According to the series relationship,(7)εt=ε1t+ε2t.Laplace transformation is carried out on formula (6) and substituted into formula (7):
(8)E1ηαsαε¯s=ηαsασ1¯s+E1σ1¯s.Carrying out the Laplace inverse transformation on formula (8), the constitutive relationship of Maxwell model can be obtained as
(9)E1ηαDαεt=ηαDασ1t+E1σ1t.Whenσ1t=σ0, Laplace transformation and inverse transformation were performed on equation (9), and the creep relationship of Maxwell model can be obtained as
(10)εt=σ0E1+σ0ηαtαΓα+1.The constitutive relationship of each element in the fractional Poynting-Thomson model is(11)Dασ1t+E1ηασ1t=E1Dαε1t,σ2t=E2ε2t.From the parallel relationship and the Laplace transformation of equation (11),
(12)ηαsασ¯s+E1σ¯s=E1ηαsαε¯s+E1E2ε¯s+E2ηαsαε¯s.Carrying out the Laplace inverse transformation on formula (12), the constitutive relationship of the fractional Poynting-Thomson model can be obtained as
(13)σt+ηαE1Dασt=E2εt+ηα+E2ηαE1Dαεt.Whenσ1t=σ0, Laplace transformation was performed on equation (13) to obtain
(14)ε¯s=E1σ0E1E2+E1+E2ηαsαs+ηασ0E1E2s1−α+ηαE1+E2s.After proving the formula (14), formula (15) was obtained:
(15)ε¯t=E1σ0E1+E2−1η−αs−1E1E2/E1+E2ηα+sα+σ0E1+E2−1sα−1E1E2/ηαE1+E2+sα.The Laplace variation expression of the two-parameter Mittag-Leffler function is(16)∫0∞e−sttαk+β−1Eα,βk±atαdt=k!sα−βsα∓ak+1Rep>a1/α.Carrying out Laplace inverse transformation on equation (15) from equation (16), the creep relationship of the fractional Poynting-Thomson model can be obtained as
(17)εt=E1σ0E1+E2ηαtαEα,α+1−E1E2E1+E2ηαtα+σ0E1+E2Eα,1−E1E2E1+E2ηαtα.


## 3.1.1. Basic Concept of Fractional Calculus

Thenth derivative of a function ft is expressed as Dαft=dαft/dtα. 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 ft of order α is defined as [33]
(1)Dt−αft=1Γα∫0tt−ξα−1fξdξReα>0,t>0,where 0 andt are the limits of integration at the left and right sides of D; Γα is the gamma function with argument α, which is defined by
(2)Γα=∫0∞e−ttα−1dtReα>0.The fractional derivative is defined by(3)Dαft=dnD−n−αftdtn,whereα≥0, n‐1≤α≤n, and n are positive integers.


## 3.1.2. Abel Dashpot

AD is a component derived from the integer-order creep element, as shown in Figure2.AD’s stress and strain have the following expressions:(4)σt=ηcdαεtdtα0≤α≤1,α∈R+,t>0,(5)εt=σcηtαΓα+10≤α≤1,α∈R+,t>0,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.Figure 3
Effect of different orders ofα on creep properties.It can be seen from Figure3 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 Figure4.Figure 4
Fractional Poynting-Thomson model.The constitutive relationship of each element in the fractional Maxwell model is(6)σ1=E1ε1t,σ1=ηαDαε2t.According to the series relationship,(7)εt=ε1t+ε2t.Laplace transformation is carried out on formula (6) and substituted into formula (7):
(8)E1ηαsαε¯s=ηαsασ1¯s+E1σ1¯s.Carrying out the Laplace inverse transformation on formula (8), the constitutive relationship of Maxwell model can be obtained as
(9)E1ηαDαεt=ηαDασ1t+E1σ1t.Whenσ1t=σ0, Laplace transformation and inverse transformation were performed on equation (9), and the creep relationship of Maxwell model can be obtained as
(10)εt=σ0E1+σ0ηαtαΓα+1.The constitutive relationship of each element in the fractional Poynting-Thomson model is(11)Dασ1t+E1ηασ1t=E1Dαε1t,σ2t=E2ε2t.From the parallel relationship and the Laplace transformation of equation (11),
(12)ηαsασ¯s+E1σ¯s=E1ηαsαε¯s+E1E2ε¯s+E2ηαsαε¯s.Carrying out the Laplace inverse transformation on formula (12), the constitutive relationship of the fractional Poynting-Thomson model can be obtained as
(13)σt+ηαE1Dασt=E2εt+ηα+E2ηαE1Dαεt.Whenσ1t=σ0, Laplace transformation was performed on equation (13) to obtain
(14)ε¯s=E1σ0E1E2+E1+E2ηαsαs+ηασ0E1E2s1−α+ηαE1+E2s.After proving the formula (14), formula (15) was obtained:
(15)ε¯t=E1σ0E1+E2−1η−αs−1E1E2/E1+E2ηα+sα+σ0E1+E2−1sα−1E1E2/ηαE1+E2+sα.The Laplace variation expression of the two-parameter Mittag-Leffler function is(16)∫0∞e−sttαk+β−1Eα,βk±atαdt=k!sα−βsα∓ak+1Rep>a1/α.Carrying out Laplace inverse transformation on equation (15) from equation (16), the creep relationship of the fractional Poynting-Thomson model can be obtained as
(17)εt=E1σ0E1+E2ηαtαEα,α+1−E1E2E1+E2ηαtα+σ0E1+E2Eα,1−E1E2E1+E2ηαtα.


## 3.2. Fractional Poynting-Thomson Model Parameter Identification

The expression of generalized Hookes law in three-dimensional form is as follows:(18)εm=σm3K,eij=Sij2G,where σm is the first invariant of the stress tensor, MPa; εm is the first invariant of strain tensor; Sij is the deviatoric stress tensor, MPa; eij 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:(19)σij=σmδij+Sij,εij=εmδij+eij,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(20)εije=Sij2G+σmδij3K.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(21)εt=G1Sij2G1+G2ηαtαEα,α+1−G1G2G1+G2ηαtα+Sij2G1+G2Eα,1−G1G2G1+G2ηαtα,with the confining pressure and axial pressure having the following relationship:(22)σ1=λσ2=λσ3.Therefore, the three-dimensional creep equation of gangue under constant pressure in the steel drum is expressed as follows:(23)εt=G1σ11−λ3G1+G2ηαtαEα,α+1−G1G2G1+G2ηαtα+σ11−λ3G1+G2Eα,1−G1G2G1+G2ηαtα.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.Table 1
Parameters of fractional Poynting-Thomson model.
σ1 (MPa)G1 (MPa)G2 (MPa)η (MPa·min)αR25.03.10.5644.670.210.9910.03.970.7498.960.190.9915.04.861.1368.210.110.9920.05.751.25100.000.110.99In 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 Figure5. 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 R2 are above 0.99.Figure 5
Fractional Poynting-Thomson model.


## 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 loadqx and the lower part supported by the foundation reaction force px, and the deflection of the beam section is expressed as wx 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.Figure 6
Fractional viscoelastic and elastic foundation beam.


### 4.2. Solution of Integer Order Solution of Elastic Foundation

Suppose the basic differential equation of elastic foundation beam deflection curve is(24)EId4wxdx4=qx−px.When the foundation is an elastic foundation, the expression ofpx is given by
(25)px=k0wx,where k0 is the elastic foundation coefficient which can be obtained by the following formula:
(26)k0=E0h0.E0 is the elastic modulus of the cushion, GPa; h0 is the thickness of the cushion, m.The relationship among the cornerθ, the bending moment M, the shear force Q, and the deflection wx of the beam section is
(27)θx=dwxdx,Mx=−EId2wxdx2,Qx=−EId3wxdx3.Taking characteristic coefficients asα0=k/4EI4, 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):
(28)wx=eα0xd1sinα0x+d2cosα0x+e−α0xd3sinα0x+d4cosα0x+q0k0,where thed1, d2, d3, and d4 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
(29)ηαE1Dα+1px,t=E2+ηα+E2ηαE1Dαwx,t.Laplace transformation of equation (24) can be expressed as
(30)EId4w¯x,sdx4+k¯0s⋅w¯x,s=q0s.Laplace transformation of equation (13) can obtain the Laplace relationship between foundation reaction force and deformation as follows:
(31)p¯x,s=E1E2+E1+E2sαηαE1+sαηαh0w¯x,s.From equation (29) and equation (31), the corresponding relationship between elastic foundation and viscoelastic foundation is derived and expressed as
(32)σt1t2+Dα=εtE2ηατα1+ηαE2+ηαE1Dα,k¯0s=E1E2+E1+E2sαηαE1+sαηαh0.From equation (28), when the foundation is an elastic material, the Laplace expression of the bending deflection equation of the thin plate is
(33)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinαx+d4cosα0x+q0k¯0⋅s.Substituting formula (32) into (33), the Laplace deflection equation of the viscoelastic foundation sheet becomes
(34)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0E1E2+E1+E2sαηα/E1+sαηαh0⋅s.By arranging the formula (34) to its final form, we can get
(35)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−α0xd3sinα0x+d4cosα0x+q0h0s−1E1+sαηα⋅E1+E2−1η−αE1E2/E1+E2ηα+sα.Separating the variables from equation (35) and extracting the common factor, we get
(36)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0h0E1+E2ηαE1s−1E1E2/E1+E2ηα+sα+ηαsα−1E1E2/E1+E2ηα+sα.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
(37)wx,t=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0h0E1+E2ηαE1tαEα,α+1E1E2E1+E2ηα⋅tα+ηαEα,1−E1E2E1+E2ηα⋅tα.


### 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 loadq0 above the roof mainly comes from the weight of potential caving strata a, q0=2.0MPa. The creep time t is 24 hours. Considering the calculation time and accuracy, the first three terms of infinite series m and n in equation (37) and the first ten terms of infinite series k in double-parameter Mittag-Leffler function (16) are calculated. The boundary conditions of the beam are as follows:
(38)wx,tx=0x=L=0,θx,tx=0x=L=0.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.Table 2
Parameters of fractional viscoelastic foundation beam.
Foundation modeld1d2d3d4Elastic (integer order)0.0008610.000165-0.126192-0.125165Fractional order viscoelastic0.0011710.001911-0.142520-0.139436The roof subsidence of elastic foundation and fractional foundation can be obtained as shown in Figure7.Figure 7
Roof subsidence of elastic foundation and fractional foundation.It can be seen from Figure7 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.


## 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 loadqx and the lower part supported by the foundation reaction force px, and the deflection of the beam section is expressed as wx 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.Figure 6
Fractional viscoelastic and elastic foundation beam.


## 4.2. Solution of Integer Order Solution of Elastic Foundation

Suppose the basic differential equation of elastic foundation beam deflection curve is(24)EId4wxdx4=qx−px.When the foundation is an elastic foundation, the expression ofpx is given by
(25)px=k0wx,where k0 is the elastic foundation coefficient which can be obtained by the following formula:
(26)k0=E0h0.E0 is the elastic modulus of the cushion, GPa; h0 is the thickness of the cushion, m.The relationship among the cornerθ, the bending moment M, the shear force Q, and the deflection wx of the beam section is
(27)θx=dwxdx,Mx=−EId2wxdx2,Qx=−EId3wxdx3.Taking characteristic coefficients asα0=k/4EI4, 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):
(28)wx=eα0xd1sinα0x+d2cosα0x+e−α0xd3sinα0x+d4cosα0x+q0k0,where thed1, d2, d3, and d4 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
(29)ηαE1Dα+1px,t=E2+ηα+E2ηαE1Dαwx,t.Laplace transformation of equation (24) can be expressed as
(30)EId4w¯x,sdx4+k¯0s⋅w¯x,s=q0s.Laplace transformation of equation (13) can obtain the Laplace relationship between foundation reaction force and deformation as follows:
(31)p¯x,s=E1E2+E1+E2sαηαE1+sαηαh0w¯x,s.From equation (29) and equation (31), the corresponding relationship between elastic foundation and viscoelastic foundation is derived and expressed as
(32)σt1t2+Dα=εtE2ηατα1+ηαE2+ηαE1Dα,k¯0s=E1E2+E1+E2sαηαE1+sαηαh0.From equation (28), when the foundation is an elastic material, the Laplace expression of the bending deflection equation of the thin plate is
(33)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinαx+d4cosα0x+q0k¯0⋅s.Substituting formula (32) into (33), the Laplace deflection equation of the viscoelastic foundation sheet becomes
(34)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0E1E2+E1+E2sαηα/E1+sαηαh0⋅s.By arranging the formula (34) to its final form, we can get
(35)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−α0xd3sinα0x+d4cosα0x+q0h0s−1E1+sαηα⋅E1+E2−1η−αE1E2/E1+E2ηα+sα.Separating the variables from equation (35) and extracting the common factor, we get
(36)w¯x,s=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0h0E1+E2ηαE1s−1E1E2/E1+E2ηα+sα+ηαsα−1E1E2/E1+E2ηα+sα.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
(37)wx,t=eα0xd1sinα0x+d2cosα0x+e−αxd3sinα0x+d4cosα0x+q0h0E1+E2ηαE1tαEα,α+1E1E2E1+E2ηα⋅tα+ηαEα,1−E1E2E1+E2ηα⋅tα.


## 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 loadq0 above the roof mainly comes from the weight of potential caving strata a, q0=2.0MPa. The creep time t is 24 hours. Considering the calculation time and accuracy, the first three terms of infinite series m and n in equation (37) and the first ten terms of infinite series k in double-parameter Mittag-Leffler function (16) are calculated. The boundary conditions of the beam are as follows:
(38)wx,tx=0x=L=0,θx,tx=0x=L=0.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.Table 2
Parameters of fractional viscoelastic foundation beam.
Foundation modeld1d2d3d4Elastic (integer order)0.0008610.000165-0.126192-0.125165Fractional order viscoelastic0.0011710.001911-0.142520-0.139436The roof subsidence of elastic foundation and fractional foundation can be obtained as shown in Figure7.Figure 7
Roof subsidence of elastic foundation and fractional foundation.It can be seen from Figure7 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 Figure8. And the relationship of maximum roof subsidence and time is shown in Figure 9.Figure 8
Roof subsidence with different time.Figure 9
Relationship of maximum roof subsidence and time.It can be seen from Figures7 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 Figure10. And the relationship of maximum roof subsidence and viscosity coefficient is shown in Figure 11.Figure 10
Roof subsidence with different time.Figure 11
Relationship of maximum roof subsidence and time.It can be seen from Figures10 and 11 that as the viscosity coefficient increases, the roof deformation decreases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, Δw1 which is smaller than Δw2. 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 Figure12. And the relationship of maximum roof subsidence and fractional order is shown in Figure 13.Figure 12
Roof subsidence with different fractional orders.Figure 13
Relationship of maximum roof subsidence and fractional order.It can be seen from Figure12 that as the fractional order increases, the roof deformation increases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, Δw1 is smaller than Δw2. 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).


## 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 Figure8. And the relationship of maximum roof subsidence and time is shown in Figure 9.Figure 8
Roof subsidence with different time.Figure 9
Relationship of maximum roof subsidence and time.It can be seen from Figures7 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 Figure10. And the relationship of maximum roof subsidence and viscosity coefficient is shown in Figure 11.Figure 10
Roof subsidence with different time.Figure 11
Relationship of maximum roof subsidence and time.It can be seen from Figures10 and 11 that as the viscosity coefficient increases, the roof deformation decreases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, Δw1 which is smaller than Δw2. 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 Figure12. And the relationship of maximum roof subsidence and fractional order is shown in Figure 13.Figure 12
Roof subsidence with different fractional orders.Figure 13
Relationship of maximum roof subsidence and fractional order.It can be seen from Figure12 that as the fractional order increases, the roof deformation increases accordingly. However, the increase in roof deformation gradually becomes smaller, for example, Δw1 is smaller than Δw2. 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



---

*Source: 8816871-2020-11-18.xml*