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 long-term 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 steady-state 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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0004](https://onlinelibrary.wiley.com/cms/asset/0311c426-0d87-41fd-8bdb-906908b31e90/dug212033-math-0004.png)
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.
Nonlinear creep damage element model.
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 steady-state creep and accelerated creep stages. The constitutive relations of each element are described below.
For the three-dimensional 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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0008](https://onlinelibrary.wiley.com/cms/asset/09b66d71-d934-4660-9152-4ed3bb0b8c89/dug212033-math-0008.png)
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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0010](https://onlinelibrary.wiley.com/cms/asset/24518ee2-cc57-40f3-bf46-2403f151b847/dug212033-math-0010.png)
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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0011](https://onlinelibrary.wiley.com/cms/asset/84bdd766-15f8-4c64-a93f-6040785b7e27/dug212033-math-0011.png)
can be expressed as follows:
(7)
where
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0013](https://onlinelibrary.wiley.com/cms/asset/416df42e-5bcc-47f5-bf70-33a07d07c124/dug212033-math-0013.png)
and
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0014](https://onlinelibrary.wiley.com/cms/asset/911e3913-a421-4eba-8b54-cc05fc8dbe92/dug212033-math-0014.png)
are the deviatoric and spherical stress tensors of the plastic slider, respectively; and
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0015](https://onlinelibrary.wiley.com/cms/asset/dcca51a3-e284-4f5b-8215-75a8afd93376/dug212033-math-0015.png)
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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0017](https://onlinelibrary.wiley.com/cms/asset/087d4319-8e5e-45a3-b74a-68455ebfd245/dug212033-math-0017.png)
is expressed as follows:
(9)
The plastic flow law should be considered to solve the constitutive equation of the viscoplastic element under three-dimensional stress. Therefore, the constitutive equation of the viscoplastic element can be expressed as follows:
(10)
where
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0020](https://onlinelibrary.wiley.com/cms/asset/387f5e0a-3653-4eae-8b51-098d9e6e6cc3/dug212033-math-0020.png)
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
![urn:x-wiley:20970668:media:dug212033:dug212033-math-0022](https://onlinelibrary.wiley.com/cms/asset/7262e8d4-4934-46c0-bd70-ef2b812cbb6a/dug212033-math-0022.png)
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 three-dimensional normal stress is the three-dimensional 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 least-squares 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 |
Fitting curves of the nonlinear creep damage model.
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.