Fractional order calculus is an extension of integer order calculus. Since the 1990s, the theory and methods of fractional order calculus have been widely applied to various fields of the natural and social sciences. In the area of viscous fluid mechanics, the introduction of fractional order calculus theory allows for more realistic theoretical models to be developed in the study of mechanical physical problems associated with real fluids, leading to accurate conclusions. The Riemann–Liouville (R-L) type fractional order calculus is commonly used in the theory of rock rheology studies (Zhou et al., 2018), but its fractional order derivatives are hyper-singular and limited in applications in engineering and technology and in physical modeling. In this paper, we use a theory of fractional order derivatives with weak singular properties proposed by the Italian geophysicist Caputo. The Caputo fractional order derivative solves the fractional order initial value problem in the definition of R-L type fractional order calculus and has been widely used in the modeling process of many practical application problems (Liu et al., 2021). The definition of Caputo fractional order derivative is D t β 0 C f t = 0 D t β − n D n f t = 1 Γ n − β ∫ 0 t t − ξ n − β − 1 f n ξ d ξ , β > 0 (2) where n is the smallest integer greater than or equal to β ; ξ is an integral variable of 0 , t ; f n ξ is the n th order derivative of function f ξ ; Γ is Gamma function, defined as Γ z = ∫ 0 ∞ t z − 1 e − t d t = 2 ∫ 0 ∞ t 2 z − 1 e − t 2 d t (3) Γ 1 + z = z Γ z z ∈ N * , Re z > 0 . (4)

The Caputo fractional operator is shown in Eq. 5: D t − β 0 C f t = 0 D t − β f t = 1 Γ β ∫ 0 t t − ξ β − 1 f ξ d ξ , β > 0 (5)

The Laplace transform formula for the Caputo fractional order derivative is L D t β 0 C f t = s β U s − ∑ j = 0 n − 1 u j 0 s β − j − 1 , n − 1 < u ≤ n (6) where U s is the Laplace transform operator of the function f t .

The theory of fractional order calculus is applied to the traditional Newtonian dashpot to construct a fractional dashpot, called the Abel dashpot. As shown in Figure 4:

Define the stress-strain equation for the Abel dashpot as σ t = η 0 d β ε t d t β (7)

Where η is the viscosity coefficient; σ t is the axial stress and ε t is the axial strain.

When β = 0 , η 0 = E , σ = E ε , representing linear elastomer, namely, Hooke elastomer; When β = 1 , σ = η 0 d ε / d t , corresponding to Newtonian dashpot and satisfying the ideal fluid. So, the physical meaning of Abel dashpot can be defined as a fluid element between Hooke elastomer and Newtonian body.

In the study of rock creep mechanics, σ t = σ = c o n s t , Eq. 7 is integrated with Caputo fractional operator, and the creep equation of fractional Abel dashpot can be expressed as ε t = σ η 0 ⋅ t β Γ 1 + β , m ≤ β ≤ m + 1 (8)

For Eq. 8, if σ / η 0 is a constant, select 0 < β ≤ 1 and β > 1 to draw the figure of strain time as follows:

As shown in Figures 2, 3, the growth rate of ε is variate non-linear follow β . Combining with Figure 1, the constitutive relationship curves of Abel dashpot can be used to describe the typical creep properties of the rock. When the stress level is below the long-term strength of the rock, the decay creep stage can be characterized by Abel daspot with 0 < β < 1 . When the stress level exceed the long-term strength of the rock, the accelerating creep stage can be characterized by Abel despot with β > 1 . Therefore, the Abel dashpot body is introduced by paralleling with the friction plate to construct a fractional order viscoplastic body, as shown in Figure 4. The stress relationship of fractional order plastic body as: σ = σ s + σ A b e l (9)

Substituting Eq. 7 into Eq. 9, the fractional order viscoplastic body creep constitutive relationship as: ε v p t = σ − σ s η 2 ⋅ t k Γ 1 + k , k > 1 (10) where η 2 is the viscosity coefficient of the viscoplastic body, and σ s is the long-term strength of the rock.

Based on fractional order calculus theory and Boltzmann superposition principle, a six-element non-linear elastic-viscoplastic creep model is proposed as shown in Figure 4 in this paper. The model consists of a Hooke elastomer, an Abel dashpot body, a Kelvin body, and a fractional order viscoplastic body in series. In this case, the instantaneous creep of the rock is characterized by Hooke elastomer. The non-linear decay creep is characterized by the Abel dashpot body. The steady-state creep is approximated by a constant strain with time, as t → ∞ , the slope of the strain-time curve k = d ε t d t = c o n s t , and the constitutive relationship for this stage is described by the conventional linear Kelvin body. Finally, the accelerating creep stage is described by the fractional order viscoplastic body.

According to the series-parallel law for the combined element model, when the stress is a constant, the stress-strain relationship for the six-element non-linear elastic-viscoplastic creep model as the following equation shows: σ = σ e = σ v e 0 = σ v e l = σ v p ε t = ε e t + ε v e 0 t + ε v e l t + ε v p t (11) where σ e , σ v e 0 , σ v e l , σ v p is the stress of Hooke elastomer, Abel dashpot body, Kelvin body and fractional order viscoplastic body, respectively. ε t e , ε t v e 0 , ε t v e l , ε t v p is the strain of Hooke elastomer, Abel dashpot body, Kelvin body and fractional order viscoplastic body, respectively.

(i) When σ < σ s , there is only transient elastic creep, decay creep and steady-state creep in the rock, no accelerating creep has occurred, and only Hooke elastomer, Abel dashpot, kelvin body at work, the stress-strain relationship of the model equation is

The creep equation of the Hooke elastomer and Abel dashpot body is ε e t = σ E 0 ε v e 0 t = σ η 0 t p Γ 1 + p (13)

For kelvin body, according to the parallel law, there is σ = σ E 1 + σ η 1 ε = ε E 1 + ε η 1 σ E 1 = E 1 ε σ η 1 = η 1 ε ˙ (14)

So, the creep constitutive equation of kelvin body is σ = E 1 ε + η 1 d ε d t (15)

Solve Eq. 15, obtain the creep equation of the linear kelvin body is ε v e l t = σ E 1 1 − exp − E 1 η 1 t (16) where σ is strain, E 1 , η 1 is elastic modulus and viscosity coefficient, respectively.

Substituting Eqs 13, 16 into Eq. 12, access the equation of ε t as ε t = σ E 0 + σ η 0 t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t , σ < σ s , 0 < p < 1 (17)

(ii) When σ ≥ σ s , the rock undergoes accelerating creep and the fractional order viscoplastic body is added to the work, the stress-strain relationship of the model is given by

Combining Eq. 10, Eq. 13 and Eq.16 and substituting them into Eq. 18, access the equation of ε t as ε t = σ E 0 + σ η 0 t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t + σ − σ s η 2 t k Γ 1 + k , σ ≥ σ s , k > 1 (19)

In summary, the non-linear creep equation for rock in a one-dimensional stress state is as follows: ε t = σ E 0 + σ η 0 t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t , σ < σ s , 0 < p < 1 σ E 0 + σ η 0 t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t + σ − σ s η 2 t k Γ 1 + k , σ ≥ σ s , k > 1 (20)

In geotechnical engineering, rocks are often in a complex three-dimensional stress state. Therefore, in order to reflect the creep properties of rocks in geotechnical engineering more accurately, the three-dimensional creep constitutive equation should be adopted. According to the theory of elasticity, the internal stress tensor σ i j of rock can be decomposed into a spherical stress tensor σ m and deviatoric stress tensor s i j under the condition of the three-dimensional stress. Similarly, the total strain tensor ε i j of rock can be decomposed into a spherical strain tensor ε m and a deviatoric strain tensor e i j , and their constitutive relationship can be expressed, respectively, as follow: σ i j = s i j + σ i j δ i j ε i j = e i j + ε m δ i j (21) where σ i j is the Kronecker delta. The relationship between different stress tensors and strain tensors is as follows: σ m = 1 3 σ 1 + σ 2 + σ 3 = 1 3 σ k k ε m = 1 3 ε 1 + ε 2 + ε 3 = 1 3 ε k k (22)

Under the condition of the three dimensional stress, assuming the total strain of the non-linear creep model is ε i j t . ε i j 0 t , ε i j v e 0 t , ε i j v e l t , ε i j v p t is the strain of Hooke elastomer, Abel dashpot body, Kelvin body and fractional order viscoplastic body, respectively. Based on the theoretical superposition principle of the component combination model, the relationship of strain is ε i j t = ε i j e t + ε i j v e 0 t + ε i j v e l t + ε i j v p t (23)

In rock creep tests, when the first order stress level applied is less than the long-term strength of the rock, the rock develops transient strain rapidly, and the constitutive relationship at this stage can be described by a Hooke elastomer.

For the Hooke elastomer, the elastic constitutive relation can be expressed by Hooke’s Law as s i j = 2 G e i j σ i j = 3 K ε i j (24) where G , K is the shear modulus and bulk modulus, respectively. The relationship between the shear modulus G , bulk modulus K , elastic modulus E , and Poisson`s ratio μ of soil is G = E 2 1 + μ K = E 3 1 − 2 μ (25)

Hence, the strain of the Hooke elastomer can be written as ε i j e t = 1 2 G 0 S i j + 1 3 K 0 σ m δ i j (26)

When the stress keep loading, the rock enters a non-linear decay creep phase, which is described by the Abel dashpot in this paper, and the three-dimensional creep constitutive equation as: ε i j v e 0 = 1 2 η 0 S i j t p Γ 1 + p (27)

Then, the rock will undergo steady-state creep, use the conventional linear Kelvin body to describe this phase of the rock creep process in this paper. Assuming that the volume change is elastic and the rheological properties are mainly in terms of shear deformation (Sun, 2007b), the three-dimensional creep constitutive equation is ε i j v e 1 = 1 2 G 1 S i j 1 − exp − G 1 η 1 t (28) where G 1 and η 1 are the shear modulus and viscosity coefficient of the Kelvin body, respectively.

When the stress deviator tensor s ≥ σ s , significant plastic deformation occurs, and the rock enters a phase of accelerating creep until it breaks down.

In a three-dimensional stress state, when stress exceeds the viscoplastic yield surface, a viscoplastic strain will be generated. Based on Perzyna’s limit stress flow law (Perzyna, 1966; Aydan, 2016), the three-dimensional creep equation for fractional order viscoplastic body in the accelerating creep phase can be obtained as ε i j v p t = 1 η 2 Φ F F 0 ∂ Q ∂ σ i j t k Γ 1 + k , k > 1 (29) where Φ F F 0 is the power function, Φ F F 0 = F F 0 m . Q is the plastic potential function. • is switch function, expressed as: 〈 Φ F F 0 〉 = 0 , F < 0 F F 0 , F ≥ 0 (30) F is the rock yield function. F 0 is the initial value of the rock yield function, generally taken as F 0 = 1 (Al-Rub et al., 2013). The exponent m is a constant and generally taken as m = 1 .

Combining Eqs 26–30 into Eq. 23, obtain the three-dimensional creep equation as follow ε i j t = 1 2 G 0 S i j + 1 3 K 0 σ m δ i j + 1 2 η 0 S i j t p Γ 1 + p + 1 2 G 1 S i j 1 − exp − G 1 η 1 t F < 0 1 2 G 0 S i j + 1 3 K 0 σ m δ i j + 1 2 η 0 S i j t p Γ 1 + p + 1 2 G 1 S i j 1 − exp − G 1 η 1 t + 1 η 2 F F 0 m ∂ Q ∂ σ i j t k Γ 1 + k F ≥ 0 (31)

The classical rock strength criterion includes Mohr-coulomb criterion, Tresca criterion, Von-Mises criterion and Drucker-Prager criterion, et al. However, the Mohr-coulomb and Tresca criterion considers only the maximum and minimum principal stress of the three principal stresses and does not consider the effect of intermediate principal stresses in the material. The Von-Mises criterion does not consider the effect of hydrostatic pressure on yielding and damage. The Drucker-Prager yield criterion improves the corner singularity problem of the Mohr-Coulomb criterion, and is also suitable for describing the yield behavior of rock materials. Therefore, the Drucker-Prager criterion is chosen as the yield criterion for rock creep analysis in this paper. In creep deformation, the creep yield deformation of rock materials mainly results from the deviator stress tensor, and the spherical stress tensor has little effect on yield deformation (Zheng and Kong, 2006), which is defined by F = J 2 − σ s / 3 (32) where J 2 is the second deviatoric stress tensor invariant. The tested material is suitable for the associated flow rule when F = Q (Moghadam et al., 2013).

In the true triaxial stress environment, the rock is stressed from three directions, and there is σ 1 > σ 2 > σ 3 σ m = σ 1 + σ 2 + σ 3 3 S 11 = σ 1 − σ m = 2 σ 1 − σ 2 − σ 3 3 (33)

Considering the triaxial creep experiment in the normal triaxial stress status, i.e., σ1 >σ2 = σ3, we can obtain σ m = σ 1 + 2 σ 3 3 S 11 = 2 σ 1 − σ 3 3 J 2 = 1 3 σ 1 − σ 3 F ∂ F ∂ σ 11 = σ 1 − σ 3 − σ s 3 (34)

Substituting Eqs 32–34 into Eq. 31, the full-stage creep strain under the three-dimensional stress state is obtained as follows: ε 11 t = σ 1 − σ 3 3 G 0 + σ 1 + 2 σ 3 9 K 0 + 1 η 0 σ 1 − σ 3 3 t p Γ 1 + p + σ 1 − σ 3 3 G 1 1 − exp − G 1 η 1 t σ 1 − σ 3 < σ s σ 1 − σ 3 3 G 0 + σ 1 + 2 σ 3 9 K 0 + 1 η 0 σ 1 − σ 3 3 t p Γ 1 + p + σ 1 − σ 3 3 G 1 1 − exp − G 1 η 1 t + 1 η 2 σ 1 − σ 3 − σ s 3 t k Γ 1 + k σ 1 − σ 3 ≥ σ s (35) Where ε 11 t represents the vertical strain of the sample under constant stress.

The establishment method of the empirical model for layered rock creep is based on the model of isotropic materials. The creep equation of layered rock is derived by introducing the influence of laminae on the mechanical properties of rock creep. For example, Tang (Tang et al., 2018) gave the relationship between the elastic modulus, creep rate and bedding angle at the same moisture content by uniaxial compression creep test as E = A e B β 1 η = D e m β (36) where A , B , C , D , m is a fitting coefficient, β is the bedding angle.

The paper introduces functional expressions of elastic modulus and creep rate, and proposes a one-dimensional creep constitutive equation for layered shale based on the Burgers model. However, the creep model for the accelerating creep stage was not given in that paper. Wang (Wang, 2020) obtained the anisotropic characteristics of long-term strength, peak strength and elastic modulus of sandstones with different bedding angles by triaxial creep tests.

In this paper, based on the research, the non-linear elastic-viscoplastic whole process creep equation for layered rock is derived analogously from the isotropic rock material creep model constitutive equation. The one-dimensional creep equation of the elastic-viscoplastic creep of layered rock is ε t = σ E 0 β + σ η 0 β t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t , σ < σ s , 0 < p < 1 σ E 0 β + σ η 0 β t p Γ 1 + p + σ E 1 1 − exp − E 1 η 1 t + σ − σ s β η 2 t k Γ 1 + k , σ ≥ σ s , k > 1 (37) where β is bedding angle, and the diagram of bedding angle is shown in Figure 5. E 0 β , η 0 β , σ s β is the functions of the effect of changes in bedding angle on the elastic modulus, viscosity coefficient and long-term strength of rock, respectively.

Accordingly, with reference to the derivation of the three-dimensional creep equation for isotropic rock materials, the three-dimensional creep equation of the elastic-viscoplastic creep of layered rock is ε 11 t = σ 1 − σ 3 3 G 0 β + σ 1 + 2 σ 3 9 K 0 β + 1 η 0 β σ 1 − σ 3 3 t p Γ 1 + p + σ 1 − σ 3 3 G 1 1 − exp − G 1 η 1 t σ 1 − σ 3 < σ s β σ 1 − σ 3 3 G 0 β + σ 1 + 2 σ 3 9 K 0 β + 1 η 0 β σ 1 − σ 3 3 t p Γ 1 + p + σ 1 − σ 3 3 G 1 1 − exp − G 1 η 1 t + 1 η 2 σ 1 − σ 3 − σ s β 3 t k Γ 1 + k σ 1 − σ 3 ≥ σ s β (38)

The plasticity theory analytical method based on the assumption of constant Poisson’s ratio (Dathe et al., 2001). It is assumed that Poisson’s ratio does not change with time and stress, and is equivalent to the value of elastic stage, μ σ , t = μ . Based on the creep constitutive equation established under one-dimensional condition, the equation can be extended from one-dimensional stress state to three-dimensional. Let the creep compliance J t in Eq. 20 is given as: J t = 1 E 0 + 1 η 0 t p Γ 1 + p + 1 E 1 1 − exp − E 1 η 1 t , σ < σ s 1 E 0 + 1 η 0 t p Γ 1 + p + 1 E 1 1 − exp − E 1 η 1 t + 1 − σ s / σ η 2 t k Γ 1 + k , σ ≥ σ s (39)

Then, Eq. 20 can be described as ε = J t σ (40)

For isotropic rock, the creep compliance substitution method can be used to obtain the basic form of the three-dimensional creep equation of rock as follows (Li et al., 2021): ε = J t A σ (41) where A is Poisson’s ratio matrix for isotropic material; ε and σ are strain tensor and stress tensor, respectively.

For Eq. 41, scholars (Aravas et al., 1995; Kou et al., 2023) have carried out a detailed solution, which is not repeated in this paper. Based on the creep constitutive model established in this paper, the three-dimensional creep constitutive equation for layered rock is derived as ε 11 t = 1 E 0 + 1 η 0 t p Γ 1 + p + 1 E 1 1 − exp − E 1 η 1 t + 1 − σ s σ η 2 t k Γ 1 + k × − μ ′ n sin 4 ⁡ θ + cos 4 ⁡ θ + sin 2 ⁡ 2 θ 2 σ x + sin 4 ⁡ θ + cos 4 ⁡ θ n + sin 2 ⁡ 2 θ 4 1 + 1 n ] σ y − μ sin 2 ⁡ θ μ ′ n cos 2 ⁡ θ σ z } (42) where μ , μ ′ is Poisson’s ratio parallel to and perpendicular to foliation plane, respectively. θ is the bedding angle, defined as the angle with the horizontal plane. Defining n = E θ = 90 ° / E θ = 0 ° , E θ = 90 ° 、 E θ = 0 ° is elastic modulus perpendicular to and parallel to foliation plane, individually. σ x , σ y , σ z are the three axes of positive pressure in the direction of the overall orthogonal coordinate axes, respectively.

The uniaxial compression non-linear creep model is shown in Eq. 20, and parameters of the model to be determined are E 0 , η 0 , p , E 1 , σ s , η 2 , k .

(i) The parameters of E 0 , σ s

In the transient elastic creep stage, E 0 could be calculated from E = σ / ε . σ s could be obtained from experimental data.

(ii) The parameters of η 0 , p , E 1 , η 1

Selecting the corresponding creep test data of σ σ ≥ σ s to perform step-by-step fitting, the operation process as follows: the first step is to take the test data from the decay creep and steady creep stages and follow the operation of (i) and (ii) in turn to obtain E 0 , η 0 , p , E 1 , η 1 , σ s , the second step is to take the test data of σ σ ≥ σ s to perform the non-linear fitting, then η 2 , k obtained.

The triaxial compression non-linear creep model is shown in Eq. 35, and the model parameters to be determined are G 0 , K 0 , η 0 , G 1 , η 1 , p , σ s , η 2 , k . Let the constant as M = σ 1 − σ 3 3 G 0 + σ 1 + 2 σ 3 9 K 0 N = σ 1 − σ 3 3 R = σ 1 − σ 3 − σ s 3 (43) then Eq. 35 can be described as ε 11 t = M + 1 η 0 N t p Γ 1 + p + N G 1 1 − exp − G 1 η 1 t σ 1 − σ 3 < σ s M + 1 η 0 N t p Γ 1 + p + N G 1 1 − exp − G 1 η 1 t + 1 η 2 R t k Γ 1 + k σ 1 − σ 3 ≥ σ s (44)

The process of parameter identification is:

(i) The parameters of G 0 , K 0 , η 0 , G 1 , η 1 , p

Based on the elastic modulus E and Poisson’s ratio μ from conventional triaxial compression tests on rocks under the same circumferential pressure, G 0 and K 0 are obtained through Eq. 25. Then, combining G 0 , K 0 into Eq. 35 to perform the non-linear fitting with the creep experimental data from the stress of σ σ < σ s .

(ii) The parameters of σ s , η 2 , k

The parameter of σ s could be calculated from triaxial creep experimental data. Selecting the full process of high-stress creep experimental data and fitting it to obtain G 0 , K 0 , η 0 , G 1 , η 1 , p firstly, then referring to the step-wise fitting method in Section 4.1.1, we could obtain η 2 , k .

The experimental data used for parameter identification in the one-dimensional creep equation is chosen from Wang (Wang et al., 2020). The test was carried out using the MTS815.02 Multifunctional Servo Test System for uniaxial compression creep testing of sandstone with graded loading, and the long-term strength given in the paper is 70 MPa. In this paper, data from the water content test set ω = 0 % is taken for parameter identification of the creep model. The graded loading scheme of the uniaxial compression creep test is 60 M P a , 70 M P a , 80 M P a , 90 M P a , 100 M P a , respectively. The contrastive analysis of the creep calculation curve and experimental data are illustrated in Figure 6.

Showing in Figure 7, the proposed creep model can accurately describe the characteristics of three phases of rock creep, the rationality is verified. The parameters of the creep model listed in Table 2.

The experimental data used for parameter identification in the three-dimensional creep equation is chosen from Ye (Ye et al., 2022). This test was carried out for triaxial compression creep tests at different water contents ω at an confining pressure of 1 MPa. The graded loading scheme of the triaxial compression creep test is shown in Table 3 and the long-term strength of rock is 2.48 M P a , 1.28 M P a , 0.88 M P a , respectively.

Using the method of the previous Section 4.1.2 for model parameters identification. When the loading stress level is less than the yield strength, the first equation of Eq. 35 is used for nonlinear fitting. The fit of the experimental data to the creep equation was obtained as shown in Figures 8–10.

When the loading stress level is greater than the yield strength, a step-wise non-linear fitting is made by the second equation of Eq. 35. The fitting curve of the moisture content under low-stress level is shown in Figure 8, and the fitting curve of the moisture content of ω = 8.47 % , ω = 12.38 % under high-stress level is shown in Figures 11, 12. The parameters of the creep model listed in Table 4.

Analyzing Figures 7, 10, the creep model not only describes the viscoelastic behavior in the decay creep and steady-state creep phases of the rock in uniaxial creep experiment and triaxial creep experiment well under low stress, but also has a good representation of the non-linear mechanical behavior in the accelerating creep phase of the rock under high-stress conditions. Furthermore, as seen in Tables 2, 4, the correlation between experimental data and fitting curves is high. The reasonableness and applicability of the fractional order creep model proposed in this paper are verified.

(i) Creep Equation of Empirical Model

The experimental data used for parameter identification in the one-dimensional creep equation is chosen from Tang (Tang et al., 2018). Selecting the data of 60 0 ( ω = 0 % ), and fitting by Eq. 37. The fitting curves are shown in Figure 11.

The parameters of the creep model listed in Table 5.

The experimental data used for parameter identification in the three-dimensional creep equation is chosen from Wang (J. Wang et al., 2020). The author took the triaxial compression creep tests on five sets of specimens with bedding angles of 0 0 , 30 0 , 45 0 , 60 0 and 90 0 under the confining pressure condition of 5 MPa. In this paper, selecting five sets of experimental data from the accelerating creep phase with different bedding angles to be fitted to identify the model parameters. The fitting curve is shown in Figures 12, 13.

The parameters of the creep model are listed in Table 6.

Based on the parameters from Table 6, giving the plotting of G 0 and 1 / η 0 with different bedding angles separately, and comparing with the experimental data in the original paper. We find that the parameters have a similar trend to the parameters from the experiment, verifying the influence of the bedding angle on the creep properties of the rock. As shown in Figures 14, 15.

(ii) Creep Equation of Plasticity Theory

Experimental data were used from the literature (Kou et al., 2023). In the literature, step-wise loading triaxial creep tests of phyllite specimens with three kinds of bedding angles (0°, 45° and 90°) are carried out with the confining pressure of 10 MPa. In this paper, the last stage accelerating creep test data of 0° and 90° layered rock is taken to validate the plasticity theory creep equation and compared with the empirical model creep equation simultaneously. The elastic mechanical parameters of the rock are shown in Table 7. For the bedding angle of 0°, the elastic modulus, poission ratio and long-term strength of rock is 26.25 G P a , 0.29,91 M P a , respectively. For the bedding angle of 90°, the elastic modulus, poission ratio and long-term strength of rock is 29.31 G P a , 0.36,91 M P a , respectively.

For Eq. 42, let F θ = − μ ′ n sin 4 ⁡ θ + cos 4 ⁡ θ + sin 2 ⁡ 2 θ 2 σ x + sin 4 ⁡ θ + cos 4 ⁡ θ n + sin 2 ⁡ 2 θ 4 1 + 1 n σ y − μ sin 2 ⁡ θ μ ′ n cos 2 ⁡ θ σ z (45)

Then Eq. 42 is expressed as ε 11 t = F θ 1 E 0 + 1 η 0 t p Γ 1 + p + 1 E 1 1 − exp − E 1 η 1 t + 1 − σ s σ η 2 t k Γ 1 + k } (46)

Substituting the elastic parameters from Table 9 into Eq. 45, then F θ = 76.6305 , θ = 0 ° 85.6391 , θ = 90 ° (47)

Substituting Eq. 47 into Eq. 46, fitting the experimental data as shown in Figure 16. For the creep equation of the empirical model, the same set of experimental data was fitted by Eq. 38 using the method described in Section 4.1.2.

The creep constitutive equations for layered rock derived by the above two methods were fitting by the same set of experimental data, and a curve fit of the experimental data to the creep equation was obtained as shown in Figure 16.

The parameters of the two models are shown in Table 7.

As seen in Figure 16 and Table 7, the fitting curves of the creep constitutive equation for layered rock derived by the two methods are almost identical and have the same fitting correlation coefficients.

In summary, the non-linear creep model derived in this paper can not only better characterize the creep properties of layered rock under low and high-stress conditions, but also better reflect the influence of different bedding angles on the creep mechanical properties of rocks. In addition, the empirical creep constitutive equation derived in this paper and the plastic theory creep constitutive equation has almost the same fitting results for the same set of experimental data, which further verifies the accuracy and applicability of the creep constitutive model for layered rock established in this paper.