# The 3D Nonstationary Creep Constitutive Model of Sandstone Based on Fractional Order

Authors: Shuguang Zhang; Wenbo Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Based on nonlinear creep characteristics of rock, triaxial creep tests were conducted on sandstone by using an MTS812.02 test machine. The creep properties of rocks under different conditions were analyzed. Establishing a creep model that describes parameter degradation is necessary to address the deterioration of rock creep parameters over time. A nonstationary Abel pot was constructed on the basis of laws of damage accumulation and deterioration. Then, the nonstationary creep model was established according to fractional order theory. The parameters of the triaxial creep test curve under different confining pressures were identified on the basis of Levenberg–Marquardt algorithm in Hengda coal mine in Fuxin, China. By comparing the model curve with the experimental curve, we found that the nonlinear model not only accurately reflected the creep characteristics of the decay and steady creep stages but also overcame the limitation of traditional creep model, that is, difficulty in describing the accelerated creep. Finally, the influence of fractional order and nonstationary parameters on rock creep was obtained by parameter sensitivity analysis. Overall, the model had high fitting degree as well as good prediction and analysis for 3D creep test data.


---

## Body

## 1. Introduction

Deep rock mass is generally in a state of high geostress, temperature, and permeability flow field. Compared with geotechnical materials under low stress conditions, rocks under high stress state exhibit various mechanical strength characteristics [1–3]. Rheological properties are observed, and the state of rock failure changes from being brittle to ductile as the confining pressure increases. Furthermore, the mechanical strength parameters have varying degrees of deterioration as stress and time increase [4].By adopting methods to improve the basic components of rheological models, scholars have achieved remarkable results in describing the rheological properties of rocks. However, when the rock is in a complex geological environment, the creep indicates a nonlinear degradation of the rock structure, as well as nonlinear changes in its physical and mechanical parameters due to changes in the stress state and time [5, 6]. Fahimifar et al. [7] addressed the limitations of the Sterpi model in predicting steady-state creep under low stress conditions by connecting a series of plastic elements, further satisfying the Mohr–Coulomb criterion. Pramthawee et al. [8] proposed an approach to extend a modified creep model to the time-dependent analysis of high rockfill dams. The prediction of dam deformation can be considerably enhanced by the proposed model. In another study, by introducing a time index, a nonlinear viscoelastic triaxial creep model is proposed to improve the prediction of the long-term deformation behavior of coal [9]. Liu et al. [10] adopted the Kachanov creep damage rate and introduced a damage variable into a rheological differential constitutive equation to derive an analytical integral solution for the creep damage equation of the Bingham model. Moreover, a new four-element creep model based on the variable-order fractional derivative and continuous damage mechanics was proposed by Tang et al. [11].The present study established a nonstationary Abel stick based on damage accumulation; this tool denormalizes other creep parameters. Then, a nonstationary creep model was established by applying fractional order theory. The sandstone of Fuxin Hengda Coal Mine was subjected to a triaxial creep test. Through a combination of this approach with the Levenberg–Marquardt algorithm, the test curves under different confining pressures were fitted. Furthermore, the correctness and rationality of the nonstationary fractional creep model were verified.


## 2. Establishment of Sandstone Fractional 3D Creep Constitutive Model

### 2.1. Establishment of Deteriorated Abel Dashpot Creep Model

The variation law of fractional order functions could be described in many forms. Riemann–Liouville calculus method was commonly used to define fractional orders. For any real numberγ, a Riemann–Liouville fractional order function exists [11] as follows: (1)
d
-
γ
f
t
d
t
-
γ
=
1
Γ
γ
∫
0
t
t
-
ξ
γ
-
1
f
ξ
d
ξ
,where γ is the order of the fractional order.Equation (1) is the γ-th order Riemann–Liouville fractional integral of the functionf (t).The Gamma function Γ(γ) is as follows:(2)
Γ
γ
=
∫
0
∞
t
γ
-
1
e
-
t
d
tGiven the fact that conventional Newtonian does not accurately describe the accelerated creep properties of the rock, the integer-order clay dashpot was transformed into a fractional clay dashpot. The Abel dashpot is shown in Figure1. The constitutive equation is as follows:(3)
σ
t
=
η
γ
d
γ
ε
t
d
t
γ
,
0
≤
γ
≤
1
,whereη is a viscosity coefficient.Figure 1
Abel dashpot.When the stress is constant, (3) is integrated to obtain the creep constitutive equation of the Abel dashpot [12].(4)
ε
=
σ
η
γ
t
γ
Γ
1
+
γ
,whereσ is the stress level value loaded by the creep test.In the rheological process, the individual creep parameters were not fixed. The rheological test of rock is a process in which rock damage accumulates continuously and mechanical parameters deteriorate. For the dashpot, the viscosity coefficientη that controls the rate of creep deformation is affected by damage accumulation and deterioration. Generally, the damage variation process of rock could be described by introducing damage variableD. To effectively describe the regularity of viscoplastic creep strain and time variation, we define the viscosity coefficient of the nonlinear Abel dashpot as the amount that changes with the accumulation of damage.(5)
η
γ
t
,
D
=
η
γ
1
-
D
,where η
γ is the viscosity coefficient of the Abel dashpot after deterioration.According to a previous study [13], the damage-evolution equation of rock satisfied the exponential function form as follows:(6)
D
=
1
-
exp
⁡
-
α
t
.Equations (5) and (6) were substituted into (3) to obtain the viscous creep coefficient equation of the Abel dashpot.(7)
σ
=
η
γ
exp
⁡
-
α
t
d
γ
ε
t
d
t
γ
.When the stress was constant, (7) was integrated to obtain the creep constitutive equation of the Abel dashpot after deterioration.(8)
ε
=
σ
η
γ
t
γ
∑
k
=
0
∞
α
t
k
Γ
1
+
γ
+
k
0
<
γ
≤
1
.Based on the coupling of external load on the rock mass, the internal mineral structure changed and recombined with time, causing the weakening of the physical and mechanical properties. Thus, the parameters were gradually degraded through time.When only the influence of the load action time is considered, the expression of the creep parameter is as follows:(9)
X
t
,
D
=
X
1
-
D
,whereX represents the creep parameter in the rock creep model.Then, (6) was substituted into (9) to obtain the nonstationary creep parameter after degradation.(10)
X
t
,
D
=
X
exp
⁡
-
α
t
.


### 2.2. Establishment of Rock Fractional 1D Creep Model

In general, the creep of rock could be divided into three stages, namely, decay creep phase, stable creep phase, and accelerated creep phase. Nishihara model (Figure2(a)) could effectively describe the creep properties of the first two stages of rock, but describing the characteristics of the accelerated creep is difficult. Fractional calculus is the theory of studying the differential and integral of arbitrary order, indicating the generalization of integer calculus to the arbitrary order. The linear Newtonian pot can be transformed into a nonlinear Newtonian pot to accurately simulate the creep process of the actual rock.Figure 2
Schematic of rheological constitutive model.
(a)
Nishihara model
(b)
Fractional creep modelIn this study, all the other creep parameters were converted into nonstationary parameters by (10). Then, the original linear components in the Nishihara model were replaced to obtain the nonstationary state. The fractional creep model is presented in Figure 2(b).The 1D state: the total strainε of the fractional model according to the rheological model theory would satisfy the following conditions [14, 15]:(11)
ε
=
ε
e
+
ε
v
e
+
ε
v
p
,whereε
e is the elastic strain,ε
ve is the viscoelastic strain, andε
vp is the viscoplastic strain.For elastic Hook body components, the parameterE
0 is degraded at different stress levels. Based on (10), the elastic strainε
e of the Hook element is obtained as follows: (12)
ε
e
=
σ
E
0
=
σ
E
0
exp
⁡
-
α
0
σ
,whereE
0 is the elastic modulus and α
0 is the coefficient related to stone properties.According to a previous study [16] and combined with nonstationary stick expression, the constitutive relation of the viscoelastic body is as follows:(13)
ε
v
e
=
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
e
-
α
2
t
/
η
1
γ
e
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
,where η
1
γ is the viscosity coefficient of Abel dashpot,E
1 is the elastic modulus of a viscoelastic model, and α
1 and α
2 are the coefficients related to stone properties.Whenσ < σ
s, the friction block does not work; thus, the whole viscoplastic element cannot work.(14)
ε
v
p
=
0Whenσ ≥ σ
s, the friction block begins to work and the nonlinear pot is deformed. When this variable is combined with the expression of nonstationary dashpot, the viscoplastic rheological equation is as follows:(15)
ε
v
p
=
σ
-
σ
s
η
2
γ
t
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
0
<
γ
≤
1
,where η
2
γ is the viscosity coefficient of the Abel dashpot and α
3 is the fractional parameter of the dashpot.In summary, (12)–(15) are substituted into (11), and the constitutive equation of the nonstationary fractional 1D creep model is obtained.Whenσ < σ
s
,
(16)
ε
=
σ
E
0
exp
⁡
-
α
0
σ
+
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(17)
ε
=
σ
E
0
exp
⁡
-
α
0
σ
+
σ
-
σ
s
η
2
γ
t
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
+
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
e
x
p
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k


### 2.3.3D Creep Model Derivation

The surrounding rock is generally under three-way stress, with few cases of one-way stress. Applying the 1D constitutive model to practical engineering scenarios is difficult. Therefore, the 1D creep model needs to be converted into a 3D model. However, all the components were obtained from the creep model derived in the 1D state. The issue of whether or not they can be used in the 3D state needs to be addressed. As previously proposed, the creep constitutive equation of rock is consistent in the 3D and 1D states [17]. Therefore, the 3D model of the nonstationary creep of rock can be derived by analogy.In the 3D state, the total strainε
ij of the creep model according to the rheological model theory satisfies the following condition:(18)
ε
i
j
=
ε
i
j
e
+
ε
i
j
v
e
+
ε
i
j
v
p
,where ε
e
i
j is the elastic deformation tensor, ε
v
e
i
j is the viscoelastic deformation tensor, and ε
v
p
i
j is the viscoplastic deformation tensor.In the 3D state, the stress tensorσ
ij can be decomposed into the spherical stress tensorσ
m and deviatoric stress tensorS
ij. Similarly, the strain tensorε
ij can be decomposed into the spherical strain tensorε
m and partial strain tensore
ij.(19)
σ
i
j
=
S
i
j
+
δ
i
j
σ
m
ε
i
j
=
e
i
j
+
δ
i
j
ε
m
,whereδ
ij is the Kronecker function, the spherical stress tensorσ
m = (σ
1 + 2σ
3) / 3, and the spherical strain tensorε
m = (ε
1 + 2ε
3) / 3;S
ij =σ
1-σ
3.The 3D stress condition of the elastic element tensorS
ij, spherical stress tensorσ
m, spherical strain tensorε
m, and partial strain tensore
ij can be obtained as follows [10, 18]:(20)
S
i
j
=
2
G
σ
-
σ
A
t
e
i
j
σ
m
=
3
K
σ
-
σ
A
t
ε
m
.Then, with the pseudo-triaxialσ
1 >σ
2 =σ
3, the elastic strain in the 3D state of the rock is as follows:(21)
ε
i
j
e
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
.The viscoelastic strain of the rock in the 3D state can also be possibly obtained:(22)
ε
v
e
=
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Although the analogy method can suitably derive the elastic and viscoelastic deformation of rock, the viscoplastic deformation also involves the yield functionF and plastic potential functionQ. Simply replacing the stressσ in the original 1D model with the stress tensorS
ij is impossible. To determine the viscoplastic strain in a 3D model, a yield function is required [19, 20]. Therefore, the constitutive relationship of the nonstationary Abel dashpot in the 3D state is as follows:(23)
ε
i
j
v
p
=
t
γ
η
2
γ
Φ
F
F
0
n
∂
Q
∂
σ
i
j
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
,whereF
0 is the initial reference value of the yield function of the rock andn is the material constant, typically n = 1.The initial yield function value of the rock is assumed to beF
0 = 1, and according to the flow law in plastic theory, (23) can be changed as follows:(24)
ε
i
j
v
p
=
F
t
γ
η
2
γ
∂
F
∂
σ
i
j
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
.Generally, F is based on the following Drucker–Plage yield function:(25)
F
=
J
2
-
σ
s
3
=
σ
1
-
σ
3
-
σ
s
3
,whereJ
2 is the stress deviator second invariant.Finally, (20)–(25) are substituted with (18) to obtain the nonstationary fractional creep constitutive equation of sandstone under 3D state.Whenσ < σ
s
,
(26)
ε
i
j
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(27)
ε
i
j
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
+
t
γ
η
2
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
σ
1
-
σ
3
-
σ
s
3
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
,whereG
1 andG
2 are the viscoelastic shear moduli,η
1 andη
2 are the shear viscous coefficients, andK is the bulk modulus.Based on (26) and (27) and rheological model theory, a radial creep constitutive model of the rock could be obtained.Whenσ < σ
s
,
(28)
ε
i
j
=
-
σ
1
-
σ
3
6
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
-
σ
1
-
σ
3
2
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(29)
ε
i
j
=
-
σ
1
-
σ
3
6
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
-
t
γ
η
2
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
σ
1
-
σ
3
-
σ
s
6
-
σ
1
-
σ
3
2
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k


## 2.1. Establishment of Deteriorated Abel Dashpot Creep Model

The variation law of fractional order functions could be described in many forms. Riemann–Liouville calculus method was commonly used to define fractional orders. For any real numberγ, a Riemann–Liouville fractional order function exists [11] as follows: (1)
d
-
γ
f
t
d
t
-
γ
=
1
Γ
γ
∫
0
t
t
-
ξ
γ
-
1
f
ξ
d
ξ
,where γ is the order of the fractional order.Equation (1) is the γ-th order Riemann–Liouville fractional integral of the functionf (t).The Gamma function Γ(γ) is as follows:(2)
Γ
γ
=
∫
0
∞
t
γ
-
1
e
-
t
d
tGiven the fact that conventional Newtonian does not accurately describe the accelerated creep properties of the rock, the integer-order clay dashpot was transformed into a fractional clay dashpot. The Abel dashpot is shown in Figure1. The constitutive equation is as follows:(3)
σ
t
=
η
γ
d
γ
ε
t
d
t
γ
,
0
≤
γ
≤
1
,whereη is a viscosity coefficient.Figure 1
Abel dashpot.When the stress is constant, (3) is integrated to obtain the creep constitutive equation of the Abel dashpot [12].(4)
ε
=
σ
η
γ
t
γ
Γ
1
+
γ
,whereσ is the stress level value loaded by the creep test.In the rheological process, the individual creep parameters were not fixed. The rheological test of rock is a process in which rock damage accumulates continuously and mechanical parameters deteriorate. For the dashpot, the viscosity coefficientη that controls the rate of creep deformation is affected by damage accumulation and deterioration. Generally, the damage variation process of rock could be described by introducing damage variableD. To effectively describe the regularity of viscoplastic creep strain and time variation, we define the viscosity coefficient of the nonlinear Abel dashpot as the amount that changes with the accumulation of damage.(5)
η
γ
t
,
D
=
η
γ
1
-
D
,where η
γ is the viscosity coefficient of the Abel dashpot after deterioration.According to a previous study [13], the damage-evolution equation of rock satisfied the exponential function form as follows:(6)
D
=
1
-
exp
⁡
-
α
t
.Equations (5) and (6) were substituted into (3) to obtain the viscous creep coefficient equation of the Abel dashpot.(7)
σ
=
η
γ
exp
⁡
-
α
t
d
γ
ε
t
d
t
γ
.When the stress was constant, (7) was integrated to obtain the creep constitutive equation of the Abel dashpot after deterioration.(8)
ε
=
σ
η
γ
t
γ
∑
k
=
0
∞
α
t
k
Γ
1
+
γ
+
k
0
<
γ
≤
1
.Based on the coupling of external load on the rock mass, the internal mineral structure changed and recombined with time, causing the weakening of the physical and mechanical properties. Thus, the parameters were gradually degraded through time.When only the influence of the load action time is considered, the expression of the creep parameter is as follows:(9)
X
t
,
D
=
X
1
-
D
,whereX represents the creep parameter in the rock creep model.Then, (6) was substituted into (9) to obtain the nonstationary creep parameter after degradation.(10)
X
t
,
D
=
X
exp
⁡
-
α
t
.


## 2.2. Establishment of Rock Fractional 1D Creep Model

In general, the creep of rock could be divided into three stages, namely, decay creep phase, stable creep phase, and accelerated creep phase. Nishihara model (Figure2(a)) could effectively describe the creep properties of the first two stages of rock, but describing the characteristics of the accelerated creep is difficult. Fractional calculus is the theory of studying the differential and integral of arbitrary order, indicating the generalization of integer calculus to the arbitrary order. The linear Newtonian pot can be transformed into a nonlinear Newtonian pot to accurately simulate the creep process of the actual rock.Figure 2
Schematic of rheological constitutive model.
(a)
Nishihara model
(b)
Fractional creep modelIn this study, all the other creep parameters were converted into nonstationary parameters by (10). Then, the original linear components in the Nishihara model were replaced to obtain the nonstationary state. The fractional creep model is presented in Figure 2(b).The 1D state: the total strainε of the fractional model according to the rheological model theory would satisfy the following conditions [14, 15]:(11)
ε
=
ε
e
+
ε
v
e
+
ε
v
p
,whereε
e is the elastic strain,ε
ve is the viscoelastic strain, andε
vp is the viscoplastic strain.For elastic Hook body components, the parameterE
0 is degraded at different stress levels. Based on (10), the elastic strainε
e of the Hook element is obtained as follows: (12)
ε
e
=
σ
E
0
=
σ
E
0
exp
⁡
-
α
0
σ
,whereE
0 is the elastic modulus and α
0 is the coefficient related to stone properties.According to a previous study [16] and combined with nonstationary stick expression, the constitutive relation of the viscoelastic body is as follows:(13)
ε
v
e
=
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
e
-
α
2
t
/
η
1
γ
e
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
,where η
1
γ is the viscosity coefficient of Abel dashpot,E
1 is the elastic modulus of a viscoelastic model, and α
1 and α
2 are the coefficients related to stone properties.Whenσ < σ
s, the friction block does not work; thus, the whole viscoplastic element cannot work.(14)
ε
v
p
=
0Whenσ ≥ σ
s, the friction block begins to work and the nonlinear pot is deformed. When this variable is combined with the expression of nonstationary dashpot, the viscoplastic rheological equation is as follows:(15)
ε
v
p
=
σ
-
σ
s
η
2
γ
t
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
0
<
γ
≤
1
,where η
2
γ is the viscosity coefficient of the Abel dashpot and α
3 is the fractional parameter of the dashpot.In summary, (12)–(15) are substituted into (11), and the constitutive equation of the nonstationary fractional 1D creep model is obtained.Whenσ < σ
s
,
(16)
ε
=
σ
E
0
exp
⁡
-
α
0
σ
+
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(17)
ε
=
σ
E
0
exp
⁡
-
α
0
σ
+
σ
-
σ
s
η
2
γ
t
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
+
σ
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
E
1
e
x
p
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k


## 2.3.3D Creep Model Derivation

The surrounding rock is generally under three-way stress, with few cases of one-way stress. Applying the 1D constitutive model to practical engineering scenarios is difficult. Therefore, the 1D creep model needs to be converted into a 3D model. However, all the components were obtained from the creep model derived in the 1D state. The issue of whether or not they can be used in the 3D state needs to be addressed. As previously proposed, the creep constitutive equation of rock is consistent in the 3D and 1D states [17]. Therefore, the 3D model of the nonstationary creep of rock can be derived by analogy.In the 3D state, the total strainε
ij of the creep model according to the rheological model theory satisfies the following condition:(18)
ε
i
j
=
ε
i
j
e
+
ε
i
j
v
e
+
ε
i
j
v
p
,where ε
e
i
j is the elastic deformation tensor, ε
v
e
i
j is the viscoelastic deformation tensor, and ε
v
p
i
j is the viscoplastic deformation tensor.In the 3D state, the stress tensorσ
ij can be decomposed into the spherical stress tensorσ
m and deviatoric stress tensorS
ij. Similarly, the strain tensorε
ij can be decomposed into the spherical strain tensorε
m and partial strain tensore
ij.(19)
σ
i
j
=
S
i
j
+
δ
i
j
σ
m
ε
i
j
=
e
i
j
+
δ
i
j
ε
m
,whereδ
ij is the Kronecker function, the spherical stress tensorσ
m = (σ
1 + 2σ
3) / 3, and the spherical strain tensorε
m = (ε
1 + 2ε
3) / 3;S
ij =σ
1-σ
3.The 3D stress condition of the elastic element tensorS
ij, spherical stress tensorσ
m, spherical strain tensorε
m, and partial strain tensore
ij can be obtained as follows [10, 18]:(20)
S
i
j
=
2
G
σ
-
σ
A
t
e
i
j
σ
m
=
3
K
σ
-
σ
A
t
ε
m
.Then, with the pseudo-triaxialσ
1 >σ
2 =σ
3, the elastic strain in the 3D state of the rock is as follows:(21)
ε
i
j
e
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
.The viscoelastic strain of the rock in the 3D state can also be possibly obtained:(22)
ε
v
e
=
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Although the analogy method can suitably derive the elastic and viscoelastic deformation of rock, the viscoplastic deformation also involves the yield functionF and plastic potential functionQ. Simply replacing the stressσ in the original 1D model with the stress tensorS
ij is impossible. To determine the viscoplastic strain in a 3D model, a yield function is required [19, 20]. Therefore, the constitutive relationship of the nonstationary Abel dashpot in the 3D state is as follows:(23)
ε
i
j
v
p
=
t
γ
η
2
γ
Φ
F
F
0
n
∂
Q
∂
σ
i
j
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
,whereF
0 is the initial reference value of the yield function of the rock andn is the material constant, typically n = 1.The initial yield function value of the rock is assumed to beF
0 = 1, and according to the flow law in plastic theory, (23) can be changed as follows:(24)
ε
i
j
v
p
=
F
t
γ
η
2
γ
∂
F
∂
σ
i
j
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
.Generally, F is based on the following Drucker–Plage yield function:(25)
F
=
J
2
-
σ
s
3
=
σ
1
-
σ
3
-
σ
s
3
,whereJ
2 is the stress deviator second invariant.Finally, (20)–(25) are substituted with (18) to obtain the nonstationary fractional creep constitutive equation of sandstone under 3D state.Whenσ < σ
s
,
(26)
ε
i
j
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(27)
ε
i
j
=
σ
1
-
σ
3
3
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
+
t
γ
η
2
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
σ
1
-
σ
3
-
σ
s
3
σ
1
-
σ
3
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
,whereG
1 andG
2 are the viscoelastic shear moduli,η
1 andη
2 are the shear viscous coefficients, andK is the bulk modulus.Based on (26) and (27) and rheological model theory, a radial creep constitutive model of the rock could be obtained.Whenσ < σ
s
,
(28)
ε
i
j
=
-
σ
1
-
σ
3
6
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
-
σ
1
-
σ
3
2
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k
.Whenσ ≥ σ
s
,
(29)
ε
i
j
=
-
σ
1
-
σ
3
6
G
0
exp
⁡
-
α
0
σ
1
-
σ
3
+
σ
1
+
2
σ
3
9
K
exp
⁡
-
α
0
σ
1
-
σ
3
-
t
γ
η
2
γ
∑
k
=
0
∞
α
3
t
k
Γ
1
+
γ
+
k
σ
1
-
σ
3
-
σ
s
6
-
σ
1
-
σ
3
2
η
1
γ
exp
⁡
-
α
1
t
∑
k
=
0
∞
-
G
1
exp
⁡
-
α
2
t
/
η
1
γ
exp
⁡
-
α
1
t
t
γ
1
+
k
γ
1
+
k
Γ
γ
1
+
k


## 3. Analysis of Results of Sandstone Triaxial Creep Tests

### 3.1. Test Plan

To reduce the differences in the samples and ensure the comparability of the tests, we derived the selected rock blocks from the same section. The rock mass, which was grayish white, was primarily composed of quartz, albite, dolomite, biotite, feldspar, and kaolinite. The samples were processed in a standard cylinder with diameter of 50 mm and height of 100 mm. The sandstone samples were depicted in Figure3(c). Postprocessing test samples were screened and examined visually in accordance with the abovementioned accuracy requirements. Samples with defects in appearance and evident differences were removed. Rock samples with similar wave velocities were selected as test pieces for the triaxial tests using the acoustic wave detection system of the MTS815.02 rock test machine.Figure 3
Rock Mechanics Test System of MTS815.02.
(a)
Test system.(
1
) Hydraulic system. (
2
) Axial and confining pressure control systems. (
3
) Loading room
(b)
Sample installation.(
4
) Sample. (
5
) Axial extensometer. (
6
) Radial extensometer. (
7
) Bearing plate and loading platform
(c)
Sandstone samplesThe main equipment used in this study was a multifunctional electrohydraulic servo-controlled rigid testing machine that was specially developed for analyzing rock and concrete (MTS Corporation, USA). The test system consisted of a loading part, a test part, and a control part. The device had three independent closed-loop servo control functions, each for controlling the axial pressure, confining pressure, and water pressure. The equipment test data were objective and reliable. The double-mean axial extensometer and circumferential extensometer achieved high test accuracy, and multiple servo-controlled methods were available. The major technical parameters of the test machine were as follows: stiffness of 7.0 × 109 N/m, maximum axial pressure of 1600 KN, maximum confining pressure of 70 MPa, and maximum pore-water pressure of 70 MPa. The test equipment is shown in Figure 3.In this study, the load control method was adopted for loading, and the loading rate was 200 N/s. Results of the triaxial compression test showed that the average compressive strength was 89.2 MPa at a confining pressure of 10 MPa. Furthermore, the average compressive strength at a confining pressure of 15 MPa was 114.7 MPa. For the sample with a confining pressure of 10 MPa, the initial load was 60% of the compressive strength. To facilitate loading, we considered the first-stage load as an integer and then loaded. The first-stage load was 55 MPa. The elastic modulus of the rock was 29.89 GPa, and the Poisson’s ratio was 0.255. For the sample with a confining pressure of 15 MPa, the initial load was 60% of the compressive strength. The first-stage load was 70 MPa. Here, the elastic modulus of the rock was 34.28 GPa, and the Poisson’s ratio was 0.236. The graded loads were increased by 5 MPa in turn and were performed in three different time groups: 10, 60, and 100 hours. If the rock sample broke during the time, the creep test was terminated; otherwise, the rock sample would continue to be tested in the next period of the creep test. When installing the sample, we first applied Vaseline on the two sample-contacting surfaces of the test piece to reduce the friction between the rigid block and end faces of the test piece. This method minimized the end effect. During the process of installing the sample, the lateral displacement sensor remained horizontal with the bottom surface of the test piece. The contact point of the axial displacement sensor was symmetrical and horizontal. Finally, the installed test piece was placed in the center of the servo tester pressure plate.


### 3.2. Axial Creep Analysis

For the samples without creep damage during the aforementioned creep test, the mechanical parameters were measured by the single-test piece method. The test data were automatically collected by the testing machine and converted into the corresponding stress–strain test output in the data acquisition system. The axial creep and time curves under different axial stresses are plotted in Figure4.Figure 4
Creep of axial strain over time under different axial stress states.
(a)
σ
3 = 10 MPa
(b)
σ
3 = 15 MPaAs shown in Figure4, the creep-time variation regularity of sandstone under different confining pressures was basically the same. At lower stress levels, the rock experienced only attenuated creep and stable creep, and the creep time was longer than 100 hours without damage. At the medium stress levels, the rock experienced a long period of attenuated and stable creep. The rate of the stable creep was greater than zero. When the internal defect development and damage accumulated to a certain extent, the creep rate suddenly increased, marking the end of the stable creep. The rock then entered the accelerated creep stage. Creep damage occurred after a short period under high stress level. The rock at the high stress level was destroyed only after 60 hours. At the high stress level, the rock underwent all three stages of creep, and the accelerated creep destabilized awhile after it began. In comparison, under the last load, the rock experienced attenuated creep and stable creep within an extremely short period, and the rock entered the accelerated creep stage, thereby producing damage quickly. After the creep had lasted for 10 hours, the rock under the highest level of stress experienced visible instability failure.Meanwhile, the duration of decay creep was gradually reduced as the stress level increased under the same confining pressure. For example, in the confining pressure of 10 MPa, the stress was 45 MPa at 32.62 hours after initiation, and the stress was 85 MPa at 0.67 hours. This phenomenon occurred because under low levels of stress, the deformation of the rock during the decay creep phase was mainly caused by the compaction of the pores and the mutual slip between the particles. Defects such as microcracks had not yet been fully formed. However, as the stress level increased, creeps deformed due to the development of fissures and cavities. Therefore, after undergoing a short phase of pore compaction and mutual slippage between the particles, defects such as cracks began to fully develop.


### 3.3. Radial Creep Analysis

The radial creep and time curves under different axial stresses are plotted in Figure5.Figure 5
Creep of radial strain over time under different radial-stress states.
(a)
σ
3 = 10 MPa
(b)
σ
3 = 15 MPaAs shown in Figure5, the variation law and characteristics of the radial creep curve are basically consistent with those of the axial creep curve. Only the radial deformation was numerically negative. Whether it was creep deformation, the instantaneous deformation or the instantaneous deformation increment gradually increased as the stress level increased. Relative to axial deformation, the rate of increase of radial transient deformation was slightly larger.


## 3.1. Test Plan

To reduce the differences in the samples and ensure the comparability of the tests, we derived the selected rock blocks from the same section. The rock mass, which was grayish white, was primarily composed of quartz, albite, dolomite, biotite, feldspar, and kaolinite. The samples were processed in a standard cylinder with diameter of 50 mm and height of 100 mm. The sandstone samples were depicted in Figure3(c). Postprocessing test samples were screened and examined visually in accordance with the abovementioned accuracy requirements. Samples with defects in appearance and evident differences were removed. Rock samples with similar wave velocities were selected as test pieces for the triaxial tests using the acoustic wave detection system of the MTS815.02 rock test machine.Figure 3
Rock Mechanics Test System of MTS815.02.
(a)
Test system.(
1
) Hydraulic system. (
2
) Axial and confining pressure control systems. (
3
) Loading room
(b)
Sample installation.(
4
) Sample. (
5
) Axial extensometer. (
6
) Radial extensometer. (
7
) Bearing plate and loading platform
(c)
Sandstone samplesThe main equipment used in this study was a multifunctional electrohydraulic servo-controlled rigid testing machine that was specially developed for analyzing rock and concrete (MTS Corporation, USA). The test system consisted of a loading part, a test part, and a control part. The device had three independent closed-loop servo control functions, each for controlling the axial pressure, confining pressure, and water pressure. The equipment test data were objective and reliable. The double-mean axial extensometer and circumferential extensometer achieved high test accuracy, and multiple servo-controlled methods were available. The major technical parameters of the test machine were as follows: stiffness of 7.0 × 109 N/m, maximum axial pressure of 1600 KN, maximum confining pressure of 70 MPa, and maximum pore-water pressure of 70 MPa. The test equipment is shown in Figure 3.In this study, the load control method was adopted for loading, and the loading rate was 200 N/s. Results of the triaxial compression test showed that the average compressive strength was 89.2 MPa at a confining pressure of 10 MPa. Furthermore, the average compressive strength at a confining pressure of 15 MPa was 114.7 MPa. For the sample with a confining pressure of 10 MPa, the initial load was 60% of the compressive strength. To facilitate loading, we considered the first-stage load as an integer and then loaded. The first-stage load was 55 MPa. The elastic modulus of the rock was 29.89 GPa, and the Poisson’s ratio was 0.255. For the sample with a confining pressure of 15 MPa, the initial load was 60% of the compressive strength. The first-stage load was 70 MPa. Here, the elastic modulus of the rock was 34.28 GPa, and the Poisson’s ratio was 0.236. The graded loads were increased by 5 MPa in turn and were performed in three different time groups: 10, 60, and 100 hours. If the rock sample broke during the time, the creep test was terminated; otherwise, the rock sample would continue to be tested in the next period of the creep test. When installing the sample, we first applied Vaseline on the two sample-contacting surfaces of the test piece to reduce the friction between the rigid block and end faces of the test piece. This method minimized the end effect. During the process of installing the sample, the lateral displacement sensor remained horizontal with the bottom surface of the test piece. The contact point of the axial displacement sensor was symmetrical and horizontal. Finally, the installed test piece was placed in the center of the servo tester pressure plate.


## 3.2. Axial Creep Analysis

For the samples without creep damage during the aforementioned creep test, the mechanical parameters were measured by the single-test piece method. The test data were automatically collected by the testing machine and converted into the corresponding stress–strain test output in the data acquisition system. The axial creep and time curves under different axial stresses are plotted in Figure4.Figure 4
Creep of axial strain over time under different axial stress states.
(a)
σ
3 = 10 MPa
(b)
σ
3 = 15 MPaAs shown in Figure4, the creep-time variation regularity of sandstone under different confining pressures was basically the same. At lower stress levels, the rock experienced only attenuated creep and stable creep, and the creep time was longer than 100 hours without damage. At the medium stress levels, the rock experienced a long period of attenuated and stable creep. The rate of the stable creep was greater than zero. When the internal defect development and damage accumulated to a certain extent, the creep rate suddenly increased, marking the end of the stable creep. The rock then entered the accelerated creep stage. Creep damage occurred after a short period under high stress level. The rock at the high stress level was destroyed only after 60 hours. At the high stress level, the rock underwent all three stages of creep, and the accelerated creep destabilized awhile after it began. In comparison, under the last load, the rock experienced attenuated creep and stable creep within an extremely short period, and the rock entered the accelerated creep stage, thereby producing damage quickly. After the creep had lasted for 10 hours, the rock under the highest level of stress experienced visible instability failure.Meanwhile, the duration of decay creep was gradually reduced as the stress level increased under the same confining pressure. For example, in the confining pressure of 10 MPa, the stress was 45 MPa at 32.62 hours after initiation, and the stress was 85 MPa at 0.67 hours. This phenomenon occurred because under low levels of stress, the deformation of the rock during the decay creep phase was mainly caused by the compaction of the pores and the mutual slip between the particles. Defects such as microcracks had not yet been fully formed. However, as the stress level increased, creeps deformed due to the development of fissures and cavities. Therefore, after undergoing a short phase of pore compaction and mutual slippage between the particles, defects such as cracks began to fully develop.


## 3.3. Radial Creep Analysis

The radial creep and time curves under different axial stresses are plotted in Figure5.Figure 5
Creep of radial strain over time under different radial-stress states.
(a)
σ
3 = 10 MPa
(b)
σ
3 = 15 MPaAs shown in Figure5, the variation law and characteristics of the radial creep curve are basically consistent with those of the axial creep curve. Only the radial deformation was numerically negative. Whether it was creep deformation, the instantaneous deformation or the instantaneous deformation increment gradually increased as the stress level increased. Relative to axial deformation, the rate of increase of radial transient deformation was slightly larger.


## 4. Nonstationary Fractional Creep Model Verification

Based on the definition of long-term rock strength by Kravcov et al. [21], the long-term strength of rock can be derived from the isochronous stress–strain curve. According to axial isochronous stress–strain curve of rock under the action of 15 MPa confining pressure, the curve exhibited apparent diverging characteristics at a stress level of 80 MPa and higher. Before this stress level was applied, the isochronous stress–strain curves were linear. The curves also began to deflect as they started to diverge. With the deflection, the creep deformation of the sandstone changed from linear to nonlinear. Therefore, the yield strength of sandstone was 80 MPa (Figure 6(a)). Similarly, under the confining pressure of 10 MPa, the long-term strength of the rock was 65 MPa (Figure 6(b)).Figure 6
Isochronous stress–strain curve.
(a)
σ
3 = 15 MPa
(b)
σ
3 = 10 MPaBy using the least squares method, we generated the parameters of the nonstationary axial creep model under different conditions as shown in Tables1 and 2. The comparison between the creep test data and nonstationary axial creep constitutive model is presented in Figures 7(a) and 7(b). Similarly, the experimental data are fitted by the radial creep model, and the radial model fitting parameters are presented in Tables 3 and 4. The comparison between the creep test data and nonstationary radial creep constitutive model is shown in Figures 7(c) and 7(d).Table 1
Fitting of axial model parameters (σ
3 = 15 MPa).
σ
1/MPa
70
75
80
85
90
95
G
0
/GPa
451.45
440.53
439.323
390.77
363.75
319.74
K/GPa
2106.81
2055.83
2052.19
1826.46
-1157.16
-912.03
G
1/GPa
1001.22
1226.95
1041.27
781.55
506.36
394.84
η
1/GPa⋅h
406.01
337.27
260.93
244.81
194.63
178.04
η
2/GPa⋅h
-
-
1150.49
1035.20
970.80
787.89
γ
0.268
0.294
0.308
0.320
0.346
0.362
α
0
13.62
12.322
12.308
11.29
10.62
9.71
α
1
0.317
0.357
0.358
0.467
0.512
0.589
α
2
-
-
-2.077
-2.038
-1.829
-1.726
α
3
-
-
-0.519
-0.775
-0.648
-0.624
R
2
0.998
0.997
0.997
0.986
0.991
0.978Table 2
Fitting of axial model parameters (σ
3 = 10 MPa).
σ
1/MPa
55
60
65
70
75
G
0
/GPa
140.67
148.99
155.58
138.50
124.21
K/GPa
213.21
255.02
297.41
-253.83
-210.82
G
1/GPa
904.10
1128.25
14301.75
1242.03
490.22
η
1/GPa⋅h
277.43
368.43
460.72
247.27
139.964
η
2/GPa⋅h
-
-
1763.72
1381.34
1247.29
γ
0.265
0.290
0.303
0.314
0.338
α
0
13.47
12.15
12.10
11.06
10.38
α
1
0.084
0.103
0.108
0.146
0.173
α
2
-
-
-0.629
-0.639
-0.618
α
3
-
-
-0.056
-0.113
-0.112
R
2
0.999
0.998
0.999
0.991
0.987Table 3
Fitting of radial model parameters (σ
3 = 15 MPa).
σ
1/MPa
70
75
80
85
90
95
G
0
/GPa
34.98
119.9
92.62
88.44
78.32
66.44
K/GPa
113.04
361.08
235.56
225
-199.56
-169.2
G
1/GPa
329.36
379.96
770.04
782
594.32
450.8
η
1/GPa⋅h
133.2
101.52
259.92
327.6
311.04
295.92
η
2/GPa⋅h
-
-
1860.48
2940.12
2208.96
2485.56
γ
0.465
0.45
0.201
0.098
0.11
0.067
α
0
11.07
12.24
5.33
2.13
1.79
1.53
α
1
0.267
0.259
0.116
0.057
0.063
0.039
α
2
-
-
-0.155
-0.076
-0.085
-0.052
α
3
-
-
-0.811
-0.648
-0.546
-0.469
R
2
0.982
0.984
0.995
0.997
0.999
0.999Table 4
Fitting of radial model parameters (σ
3 = 10 MPa).
σ
1/MPa
55
60
65
70
75
G
0
/GPa
78.32
66.44
55.66
47.74
34.54
K/GPa
199.56
169.2
141.48
121.32
87.72
G
1/GPa
594.32
450.8
321.08
195.96
45.08
η
1/GPa⋅h
311.04
295.92
189.36
158.4
97.92
η
2/GPa⋅h
-
-
1004.76
875.28
10.44
γ
0.111
0.068
0.139
0.078
0.445
α
0
1.8
1.55
1.89
1.48
0.95
α
1
0.064
0.039
0.08
0.045
0.255
α
2
-
-
-0.155
-0.076
-0.085
α
3
-
-
-0.288
-0.648
-0.546
R
2
0.983
0.985
0.996
0.998
0.997Figure 7
Comparison of fitted creep model curves with experimental data.
(a)
σ
3 = 15 MPa (axial)
(b)
σ
3 = 10 MPa (axial)
(c)
σ
3 = 10 MPa (radial)
(d)
σ
3 = 15 MPa (radial)As shown in Figure7, the actual sandstone creep test data and fitted values from the model are similar. The correlation coefficient was above 0.98, indicating an ideal fitness of the model. This model overcame the limitation of the traditional visco-elastoplastic model to describe the accelerated creep stage, showing the deformation characteristics during all three stages of sandstone creep. Therefore, the nonstationary nonlinear constitutive model is a correct and rational approach to improve the traditional visco-elastoplastic model.To prove that the model established in this study could effectively describe the accelerated creep characteristics compared with the traditional viscoelastic-plastic model, we compared the traditional model curves with the model in this study and the experimental data (Figure8).Figure 8
Model comparison.
(a)
Stress level 95 MPa (axial)
(b)
Stress level 85 MPa (axial)
(c)
Stress level 85 MPa (radial)
(d)
Stress level 95 MPa (radial)As shown in Figure8, the theoretical value of the creep model agrees well with the experimental results. It not only fully reflected the instantaneous elastic deformation of red sandstone after loading but also reflected the attenuation creep and stable creep. In particular, this model overcame the difficulty of the traditional Nishihara model in describing the accelerated creep phase and successfully fitted the deformation of the sandstone accelerated creep stage. The fitting result is superior to that of the traditional Nishihara model.


## 5. Conclusions

Based on the deterioration of rock creep parameters, the parameters satisfied the negative exponential change law, and the nonstationary Abel clay pot under the condition of damage accumulation was established. Finally, a nonstationary fractional 3D creep model was developed.Under low stress levels, rock deformation during the decay creep phase was mainly caused by the compaction of the pores and the mutual slip between the particles. Defects such as microcracks had not yet been fully formed. However, as the stress level increased, creeps deformed mainly because of the development of fissures and cavities. Therefore, after undergoing a short phase of pore compaction and mutual slippage between the particles, defects such as cracks began to fully develop.The new model proposed in this study overcame the disadvantage of the traditional visco-elastoplastic model in describing the accelerated creep stage. This model also showed the deformation characteristics during the three stages of sandstone creep. Thus, the nonstationary nonlinear constitutive model is a correct and rational approach to improve the traditional visco-elastoplastic model.The influence of the change in mechanical parameters with time on the rock structure degradation was investigated. By including the rock parameter degradation, the model is closer to the true state in which the surrounding rock creep occurred. This method is a novel approach to generate rock creep models and assist with practical engineering applications.



---

*Source: 6031842-2019-05-14.xml*