4.1 Derivation of the constitutive equation
Numerous studies have shown that the creep process of rock is also a process of damage accumulation, and the viscoplastic deformation of rock exhibits nonlinear characteristics under longterm loading (X. Li et al., 2022; Ma et al., 2021; X. Yang & Jiang, 2022). However, the current classical element model is mostly used to describe the linear constitutive relationship. Although it can characterize decay creep and steadystate creep, it cannot reflect the accelerated creep process of rock. Therefore, in this paper, the damage factors in the creep process are considered, and the nonlinear damage viscous element with damage variable Ds is introduced to describe the nonlinear viscoplastic deformation of the red sandstone.
The element is essentially a nonlinear Newton body. According to Newton's law of viscosity, its constitutive equation is
(2)
where
is the creep rate and
η
1(
D
s) is the damage viscosity coefficient. The damage viscosity coefficient can be expressed as follows:
(3)
where
η
1 is the initial viscosity coefficient at the viscoplastic stage and
D
s is the damage variable. Accordingly, the damage variable has the following relationship with time (W. Yang et al.,
2014):
(4)
where
α is the coefficient related to rock material properties;
t denotes time; and 0 <
D
s < 1. With Equations (
3) and (
4) being substituted into (
2), the following Equation (
5) can be obtained:
(5)
A physical model for describing the nonlinear creep of the red sandstone can be constructed by combining the results of the loading‒unloading creep test with the nonlinear damage viscous element. The NCD model includes an instantaneous elastic element, an instantaneous plastic element, a viscoelastic element, and a viscoplastic element, all of which are connected in series as shown in Figure 5.
In Figure
5, the instantaneous elastic element is used to describe the instantaneous elastic stage in the creep test curve. The instantaneous plastic element is used to describe the instantaneous plastic stage. The viscoelastic element is used to describe the decay creep stage. The viscoplastic element is used to describe the steadystate creep and accelerated creep stages. The constitutive relations of each element are described below.
For the threedimensional elastic element, its stress tensor
σij can be decomposed into a spherical stress tensor
δijσ
m and a deviatoric stress tensor
sij, which change the volume and shape of the object, respectively. The instantaneous elastic element strain
is expressed as follows:
(6)
where
K and
G denote the bulk modulus and the shear modulus, respectively.
According to the loading‒unloading creep test curve of the red sandstone, obvious instantaneous plastic deformation occurs in the sandstone sample at all stress levels. Therefore, by assumption, the threshold
of the plastic slider is small, and the element will deform under any stress level. In this case, the constitutive equation of the instantaneous plastic element can be obtained by analogy with the constitutive equation of the elastic element. Thus, the instantaneous plastic element strain
can be expressed as follows:
(7)
where
and
are the deviatoric and spherical stress tensors of the plastic slider, respectively; and
is the mean normal stress of the slider. When
i =
j = 1,
(8)
The spherical stress tensor mainly causes the volume change of the sample. However, in the viscoelastic stage, the sample volume almost remains constant. Accordingly, the creep deformation caused by the spherical stress tensor can be neglected. The deformation of the viscoelastic element is mainly contributed by the deviatoric stress tensor. The viscoelastic element strain
is expressed as follows:
(9)
The plastic flow law should be considered to solve the constitutive equation of the viscoplastic element under threedimensional stress. Therefore, the constitutive equation of the viscoplastic element can be expressed as follows:
(10)
where
denotes the strain rate of the viscoplastic element;
F refers to the yield function;
Q represents the plastic potential function;
F
0 denotes the initial value of the yield function; The P function has the following properties:
(11)
For rock materials,
x is usually 1. When
F ≥ 0, it can be known from relevant flow laws that
F =
Q. The viscoplastic element strain
can be obtained by integrating both sides of Equation (
10).
(12)
The initial value of yield function
F
0 is set to 1, and according to the Mises yield criterion
(13)
where
J
2 is the second invariant of the deviatoric stress tensor, MPa.
Substituting Equation (
13) into Equation (
12) yields the following:
(14)
Since each element is connected in series, the total strain variable of the model should be the sum of the strain variables of each element. The constitutive equation of the NCD model can be obtained by combining Equations (
6), (
7), (
9), and (
14). The total strain
εij is expressed as follows:
(15)
(16)
In the triaxial compression test, the threedimensional normal stress is the threedimensional principal stress, which is
σ
1,
σ
2, and
σ
3, and
σ
2 =
σ
3. Accordingly, it can be obtained when
i =
j = 1:
(17)
By substituting Equations (
8) and (
17) into Equations (
15) and (
16), the axial constitutive equations of the NCD model can be obtained as follows:
(18)
(19)
4.2 Determination of the NCD model parameters through the inversion method
To further determine the unknown parameters in the NCD model, the model parameters must be inverted by using the triaxial creep test results. At present, for nonlinear function parameter inversion problems, the leastsquares method is most commonly used. Since many unknown parameters are in the model, MATLAB is used to determine the NCD model parameters so as to avoid setting their initial values. The inversion results of the parameters at all stress levels are shown in Table 2, and the fitting curves are shown in Figure 6.
Table 2. Parameters of the NCD model.
Deviatoric stress (MPa) 
K1 (GPa) 
G1 (GPa) 
K2 (GPa) 
G2 (GPa) 
G3 (GPa) 
η1 (GPa·h) 
η2 (GPa·h) 
α (10−3) 
R2 
5 
11.14 
33.26 
11.67 
34.39 
39.32 

122.3 

0.9822 
10 
7.96 
19.73 
8.29 
20.51 
85.28 

506.9 

0.9583 
15 
6.87 
6.08 
6.84 
3.91 
17.41 

225.2 

0.9857 
20 
6.36 
8.16 
6.45 
8.41 
55.11 

461.3 

0.9755 
25 
7.83 
5.80 
7.83 
6.21 
65.74 
300.0 
300.0 
5.420 
0.9933 
27 
5.46 
6.05 
5.46 
6.05 
36.38 
350.0 
350.0 
1.160 
0.9856 
30 
5.23 
5.23 
5.23 
5.22 
18.60 
337.9 
193.1 
0.001 
0.9959 
35 
50.50 
1.19 
50.50 
50.50 
5.56 
6867.0 
42.5 
676.200 
0.9986 
Although this model contains seven elements, compared with other models, it has fewer parameters, thus being more convenient to select, and better describing the properties of rock at various creep stages. According to Table 2 and Figure 6, all correlation coefficients of the creep test curves fitted by the NCD model are above 0.95, and the fitted curves are highly consistent with the test data, indicating that the NCD model can better reflect the creep characteristics of the red sandstone.