A Nonlinear Creep Damage Model Considering the Effect of Dry-Wet Cycles of Rocks on Reservoir Bank Slopes

Abstract

Water has a crucial effect on the time-dependent behavior of rocks. The long-term cyclical fluctuations of reservoir water level lead to dry–wet (DW) cycles of rocks on reservoir bank slopes, making this influential factor more complex. To deeply understand the time-dependent behavior of rocks under DW cycles, argillite from the reservoir bank slope of Longtan Hydropower Station was used to perform a series of triaxial creep tests. Subsequently, based on analysis of creep test results after different DW cycles, a damage nonlinear Burgers viscoelastic-plastic (DNBVP) model considering the effect of saturation–dehydration cycles was proposed by introducing a nonlinear viscoplastic body and a damage variable describing DW cycles. Then, the three-dimensional creep equations of the new model were derived and its creep parameters were identified. Comparison between the theoretical curves and the test results shows that the theoretical curves of the DNBVP model were able to describe rock creep tests results after different DW cycles. Furthermore, by comparing classical creep models with the proposed model, it was found that the DNBVP model can accurately reflect the nonlinear characteristics of rocks at the accelerated creep stage. Finally, the sensitivity of the DNBVP model was analyzed and discussed, and three-dimensional central difference expressions necessary for secondary development of the new model were also derived in detail. The proposed new model with secondary development may provide a basis for improving the geotechnical design of reservoir bank slopes and the control of reservoir bank landslides.

1. Introduction

The time-dependent behavior of rocks is closely related to the safety of geotechnical engineering [1,2]. Numerical studies of engineering practices and other studies show that rock instability and failure have direct relationships with time [3,4], and variation in the mechanical properties of rocks over time can easily lead to natural disasters, especially landslides [5,6,7,8,9].

Water is considered a crucial factor in the time-dependent behavior of rocks [10,11,12,13]. With the construction and operation of numerous hydropower projects in China’s western regions, a series of high and steep reservoir bank slopes have appeared [14,15]. With reservoir impoundment, the reservoir water level undergoes dramatic cyclic fluctuations every year, resulting in dry–wet (DW) cycles affecting rocks in the hydro-fluctuation belts of reservoir bank slopes. The effect is a cumulative process, with the damage and deterioration of rocks accelerating over time [16,17], leading to the development of secondary geological disasters; therefore, reservoir areas tend to become severely afflicted areas, where landslide disasters occur frequently [18,19]. Thus, it is of great significance to fully understand the effect of DW cycles on the time-dependent behavior of rocks.

At present, studies on the mechanical behavior of rocks under the effect of DW cycles have focused mainly on uniaxial compressive strength [20], shear strength [15], dynamic compressive strength [21], tensile strength [17], slake durability [22], fracture toughness [23] and triaxial compression strength [24,25]. Previous results have shown that the mechanical properties of rocks deteriorate to some degree after DW cycles. For example, Özbek [20] performed 50 DW cycles on ignimbrites and found that the uniaxial compressive strength decreased by 16.08%. Zhang et al. [15] conducted a laboratory study on the shear strength of argillaceous siltstones after different DW cycles (0, 10, 20, 40 and 70 cycles) and concluded that the cohesion of argillaceous siltstones was most susceptible to DW cycles, which decreased from 349 kPa in the initial state to 200 kPa after 70 DW cycles. Hua et al. [23] used Brazilian disc samples to investigate the influence of different DW cycles (zero, one, three, five and seven cycles) on the fracture toughness and crack propagation of rocks, and reported that the failure behavior of sandstone changed from fragility to ductility under a high number of wet–dry cycles. Torres-Suarez et al. [24] found that after several DW cycles, the elasticity modulus of most rock samples, as well as the shear modulus of almost all samples, experienced degradation. For example, the shear modulus was about 60% of the initial value after three weeks. Wang et al. [25] studied the mechanical properties of sandstone under DW cycles and proposed an improved Duncan–Chang model which could describe the deformation characteristics of rocks at the initial compression stage and post-peak residual stage during DW cycles. Moreover, they carried out secondary development using the FLAC3D platform and constructed a dynamic link calculation program for this model using Microsoft Visual C++. However, the time-dependent behavior of rocks has not been investigated sufficiently by previous studies.

The creep test is considered the preferred method for obtaining the characteristics of the time-dependent behavior of rocks [26,27]. For example, Xu et al. [27] proposed a nonlinear creep model based on creep tests of green schist samples from Jinping dam foundation (Hohai model) and carried out secondary development of the numerical program of this model. Furthermore, they made use of examples to verify the creep law of rock mass at the dam abutment of the Jinping Hydropower Station at normal storage levels. Grgic and Amitrano [10] and Liu et al. [28] tested the creep properties of rocks under different saturation conditions and analyzed the effect of water saturation on the creep behavior of rocks, concluding that the creep effect became more significant with the increase of saturation. Hashiba and Fukui [29] and Hashiba et al. [12] performed strength tests and creep tests under uniaxial compression on andesite in dry and wet conditions and pointed out that water weakened the rock strength and reduced the creep lifetime. Moreover, a simple equation expressing the test results was developed. Deng et al. [16] conducted a series of creep tests on typical red-bed soft rock samples, which were subjected to different DW cycles to simulate the cyclical change process of reservoir water levels from the Three-Gorges Reservoir Area. They found that in the first three DW cycles, the strength of red-bed soft rock decreased rapidly and then tended to be stable. As the number of DW cycles increased, the third creep stage of the samples exhibited a trend of gradual increase, as the sample’s plasticity was enhanced under the effect of the DW cycles. Tang et al. [13] conducted creep tests on sandstone with different soaking durations and found that, compared with the dry state, the minimum creep strain rate increased while the failure time decreased with the soaking duration increased from four days to six days. With the soaking duration increased to eight days, the above phenomenon exhibited no significant change, suggesting that the sandstone was fully saturated after six days of soaking. These studies have enriched our knowledge of the creep properties of rocks under the effect of water–rock interactions. However, there are still some aspects worth exploring: (i) Existing studies on creep tests fail to obtain a deep understanding of the creep properties of rocks after DW cycles. (ii) Existing rock creep models rarely consider the damage effect of DW cycles, and there is a lack of constitutive models describing the nonlinear creep properties of rocks after different DW cycles. (iii) For creep models considering the effect of DW cycles, it is necessary to derive equations necessary for secondary development.

To address the problems above, argillite obtained from the reservoir bank slope of Longtan Hydropower Station was used to perform a series of triaxial creep tests after different DW cycles in this study. Firstly, the experimental data obtained in the creep tests were analyzed. On this basis, a nonlinear damage model accounting for the effects of saturation–dehydration cycles was proposed. Furthermore, the new model was verified and parameter sensitivity analyses of the proposed model were performed. Finally, the three-dimensional central difference expressions necessary for secondary development of the new model was derived in detail.

2. Triaxial Creep Test to Investigate the Effect of Dry–Wet Cycling on Creep

2.1. Preparation of Test Rock Specimens

Argillite was obtained from the reservoir bank slope of Longtan Hydropower Station in China (see Figure 1). In the dam site area, argillite makes up about 30.8% of the strata of the slope, and argillite in the “hydro-fluctuation belt” (Figure 1c) is subject to obvious creep under the influence of reservoir water level fluctuations. The sampling site is shown in Figure 1d. After the drilling of the boreholes, argillite samples (predominantly gray) were obtained in situ. Then, rock samples were packed and transported to the laboratory and processed into cylinders with a diameter of 50 mm and height of 100 mm, according to the GB/T50266-2013 standard for test methods of engineering rock mass [30], as shown in Figure 1e.

Table 1. The conventional physical property parameters of samples.

ρ (g/cm3)	ρd (g/cm3)	ρs (g/cm3)	W (%)	Ws (%)	Porosity (%)
2.48	2.38	2.57	4.1	7.4	7.2

Notes: ρ = Natural density, ρ_d = Dry density, ρ_s = Saturated density, W = Natural water content, W_s = Saturated moisture content.

lists the basic mechanical parameters of the argillite samples.

2.2. Experimental Procedure

Landslides frequently occur when the reservoir bank with a rock mass is in saturated conditions. Therefore, conventional triaxial tests and triaxial creep tests on rock samples after DW cycles were all performed in a saturated state in this study. In order to achieve saturation, a vacuum saturator was used to saturate the rock samples, according to the procedure mentioned in [31].

In the present study, specimens underwent 0 (saturated state), 1, 5, 10, 15 and 20 DW cycles, respectively. The DW cycles included a saturation process and a dehydration process. The saturation process of rock samples was as follows: firstly, a vacuum air pump was used to seal and saturate the rock samples, then a negative pressure of 0.1 MPa was maintained for 12 h to ensure that the samples were fully saturated. For the dehydration process, according to the standard testing method for engineering rock masses (GB/T50266-2013) [30], the samples were placed in an oven for 12 h at 104 °C to remove water completely and then were cooled naturally. A triaxial compressive strength test (with confining pressure of 4 Mpa) and a triaxial rheological test (with confining pressure of 4 Mpa) were performed after each saturation–dehydration cycle.

For the triaxial creep tests, a rock temperature-stress-seepage coupling triaxial creep meter produced by Top Industrie was used. The test procedure of this device is described in detail in [32]. During the experiment, the confining pressure was maintained at 4 MPa, and axial loads of the rock samples were applied level by level [33,34], which were approximately 30%, 45%, 60%, 75% and 85%, respectively, of the triaxial compressive strength. The loading scheme is shown in

Table 2. Loading scheme for argillite creep test.

Lithology	Load Level (MPa)	Axial Load (MPa)
		n = 0	n = 1	n = 5	n = 10	n = 15	n = 20
Argillite	First level	13.20	12.00	7.80	6.40	4.80	5.40
	Second level	19.80	18.00	11.70	9.60	7.20	8.10
	Third level	29.70	27.00	17.55	14.40	10.80	10.80
	Fourth level	39.60	36.00	23.40	19.20	14.40	13.50
	Fifth level	49.50	45.00	29.25	24.00	18.00	15.30

.

2.3. Test Results

Figure 2 shows the results of creep tests after different saturation–dehydration cycles. Deviatoric stress in this figure is the difference in axial load and confining pressure.

The strain-time and stress-time curves of argillite samples are presented in Figure 2. It can be seen that: (i) At the same level of deviatoric stress, the strain after stabilization of the creep curve increases nonlinearly as the number of DW cycle increases. For example, after the first level of deviatoric stress is applied, the values of stabilized creep strain after 0, 1, 5, 10, 15 and 20 saturation–dehydration cycles are about 0.0387%, 0.0404%, 0.0413%, 0.0421%, 0.0433% and 0.0440%, respectively. (ii) Under low axial loads, the rheological curves only exhibit two stages (decelerated stage and stable stage). However, three rheological stages—decelerated stage, stable stage and accelerated stage—may appear only when the axial load approaches or reaches the critical failure value of rocks [33,35]. (iii) As shown in Figure 2a–c, when the number of DW cycle is small, e.g., n = 0, 1 or 5, obvious shear failure zones are observed after sample failure, and conjugate shear fractures occur, with approximately X-shaped sample fracture surfaces. However, when n = 10, 15 or 20, the macromorphology of the samples after creep failure is very complex. As shown in Figure 2d–f, several shear fractures and obvious axial and horizontal tensile cracks and flaws were observed. A reasonable explanation of the above phenomenon is that the mechanical properties of the samples are degraded as the number of DW cycles increases, which also results in propagation of numerous micro-fissures inside the samples. In short, the macromorphology of the samples after creep failure becomes more complicated as the number of DW cycles increases.

Figure 3 shows the axial strain/strain rate/time curves after one saturation–dehydration cycle when the fifth level of the axial load is applied. In Figure 3, $t_{\mathrm{S}}$ is the moment when the rock samples enter the stable stage from the decelerated stage, and ${\epsilon }_{\mathrm{S}}$ is the strain value at that moment. $t_{\mathrm{P}}$ is the moment when the rock samples enter the accelerated stage from the stable stage, and ${\epsilon }_{\mathrm{P}}$ is the strain value at the dividing point between the stable and accelerated stages. $t_{\mathrm{F}}$ is the moment when failure of the rock samples occurs, and ${\epsilon }_{\mathrm{F}}$ is the strain value at that moment. At the decelerated stage, the strain continuously accumulates and the accumulation gradually increases, but the strain increase rate reduces continuously, i.e., $\epsilon >0$, $\dot{\epsilon }>0$, $\ddot{\epsilon }<0$. At the stable stage, the strain continuously accumulates over time, but the strain increase rate is constant, i.e., $\epsilon >0$, $\dot{\epsilon }>0$, $\ddot{\epsilon }=0$. At the accelerated stage, the strain develops rapidly, and the strain increase rate shows a trend of increase, i.e., $\epsilon >0$, $\dot{\epsilon }>0$, $\ddot{\epsilon }>0$. From the moment of $t_{\mathrm{P}}$, micro-fissures and pores in the rock samples propagate, converge and become connected. And at the moment $t_{\mathrm{F}}$, complete failure of the rock samples occurs.

3. Establishment of Creep Model

The Burgers model, one of the most effective constitutive models, is widely adopted to describe the creep behavior of rocks [29,36]. This model is capable of accurately describing the creep law of rocks above long-term strength [37,38]. However, the elements in the Burgers model are linear, so it is not capable to reproduce the nonlinear viscoelastic-plastic creep properties of rocks reasonably [39,40]. It is noteworthy that the irreversible viscoplasticity is considered a critical factor resulting in creep failure of rocks [41,42]. The Burgers model does not have a yield threshold; thus, it cannot describe viscoplastic deformation when the applied stress exceeds the long-term strength [33]. Based on this understanding, we proposed a nonlinear rheological model that can reflect the viscoelastic-plastic properties of rocks was proposed based on the Burgers model, as shown in Figure 4. The proposed model consists of incorporating a nonlinear viscoplastic body (consisting of a nonlinear viscous element and a plastic element connected in parallel) in the Burgers model. The new model is called the nonlinear Burgers viscoelastic-plastic (NBVP) model.

As can be seen in Figure 3, the rheological curves of the rock samples at the accelerated stage increase nonlinearly. At this stage, considering the change of parameters of the creep model over time would be an ideal method for describing the creep properties of the rock samples [43,44]. Therefore, the viscosity coefficient ${\eta }_3$ in Figure 4 should decrease over time [34,45]. The viscosity coefficients in Figure 4 are expressed in section. At stages other than the accelerated stage, the viscosity coefficient ${\eta }_3$ remains constant. At the accelerated stage, the viscosity coefficient ${\eta }_3$ is subject to nonlinear processing, and the constitutive equation of the fourth part in Figure 4 is

$${\epsilon }_4=\frac{({\sigma }_0-{\sigma }_{\mathrm{S}})(\exp <\frac{t-t_{\mathrm{P}}}{t_0}{>}^m-1)}{{\eta }_3}$$

(1)

where $t_0$ is the unit reference time (the value of which is one), $t_0$ has the effect of maintaining dimensional uniformity of the time unit [46]; m is a parameter to be determined through rheological experiment; and $<>$ is the conditional judgment switching function, in the equation:

$$<\frac{t-t_{\mathrm{P}}}{t_0}>=\{\begin{array}{l}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}t\le t_{\mathrm{P}}\\ \frac{t-t_{\mathrm{P}}}{t_0}\hspace{1em}\hspace{1em}\hspace{1em}t>t_{\mathrm{P}}\end{array}$$

(2)

3.1. One-Dimensional Creep Equations of the NBVP Model

One-dimensional creep equations of the NBVP model, based on the Laplace transform and the inverse Laplace transform [33], are as follows:

(1)

When ${\sigma }_0<{\sigma }_{\mathrm{S}}$,

$$\epsilon (t)=\frac{{\sigma }_0}{E_1}+\frac{{\sigma }_0}{E_2}(1-e^{-\frac{E_2}{{\eta }_1}t})+\frac{{\sigma }_0}{{\eta }_2}t$$

(3)

(2)

When ${\sigma }_0\ge {\sigma }_{\mathrm{S}}$, $t\le t_{\mathrm{P}}$,

$$\epsilon (t)=\frac{{\sigma }_0}{E_1}+\frac{{\sigma }_0}{E_2}(1-e^{-\frac{E_2}{{\eta }_1}t})+\frac{{\sigma }_0}{{\eta }_2}t+\frac{{\sigma }_0-{\sigma }_{\mathrm{S}}}{{\eta }_3}t$$

(4)

(3)

When ${\sigma }_0\ge {\sigma }_{\mathrm{S}}$, $t>t_{\mathrm{P}}$,

$$\epsilon (t)=\frac{{\sigma }_0}{E_1}+\frac{{\sigma }_0}{E_2}(1-e^{-\frac{E_2}{{\eta }_1}t})+\frac{{\sigma }_0}{{\eta }_2}t+\frac{({\sigma }_0-{\sigma }_{\mathrm{S}})(\exp <\frac{t-t_{\mathrm{P}}}{t_0}{>}^m-1)}{{\eta }_3}$$

(5)

3.2. Three-Dimensional Creep Equations of the NBVP Model

The triaxial creep stress state is ${\sigma }_2={\sigma }_3$ which remains constant, and the confining pressure remains constant after each level of load ${\sigma }_1$ is applied. The rheological equation of the NBVP model under triaxial creep test is shown as follows:

$$\epsilon (t)\{\begin{array}{l}=\frac{{\sigma }_1-{\sigma }_3}{3G_1}+\frac{{\sigma }_1-{\sigma }_3}{3G_2}(1-e^{-\frac{G_2}{{\eta }_1}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2}t+\frac{{\sigma }_1+2{\sigma }_3}{9K} ({\sigma }_1-{\sigma }_3<{\sigma }_{\mathrm{S}})\\ =\frac{{\sigma }_1-{\sigma }_3}{3G_1}+\frac{{\sigma }_1-{\sigma }_3}{3G_2}(1-e^{-\frac{G_2}{{\eta }_1}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2}t+(\frac{{\sigma }_1-{\sigma }_3-{\sigma }_{\mathrm{S}}}{2{\eta }_3})t\\ +\frac{{\sigma }_1+2{\sigma }_3}{9K} ({\sigma }_1-{\sigma }_3\ge {\sigma }_{\mathrm{S}} , t\le t_{\mathrm{P}})\\ =\frac{{\sigma }_1-{\sigma }_3}{3G_1}+\frac{{\sigma }_1-{\sigma }_3}{3G_2}(1-e^{-\frac{G_2}{{\eta }_1}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2}t+(\frac{{\sigma }_1-{\sigma }_3-{\sigma }_{\mathrm{S}}}{2{\eta }_3})\cdot (\exp <\frac{t-t_{\mathrm{P}}}{t_0}{>}^m-1)\\ +\frac{{\sigma }_1+2{\sigma }_3}{9K} ({\sigma }_1-{\sigma }_3\ge {\sigma }_{\mathrm{S}} , t>t_{\mathrm{P}})\end{array}$$

(6)

where K is the bulk modulus of argillite; $G_1$ and $G_2$ are the shear modulus corresponding to $E_1$ and $E_2$ in the three-dimensional stress state [47].

3.3. Analysis of the Damage Law under the Effect of Saturation–Dehydration Cycles

Defects such as holes, micro-cracks and micro-fissures in rocks propagate inevitably under the effect of external factors such as saturation–dehydration cycles [25]. The damage of rocks also means the development of micro-defects inside rocks [17].

Based on the theory of damage mechanics [48,49], the damage factor $D$ was introduced. It is assumed that (i) the damage factor of the rock samples gradually increases as the number of saturation–dehydration cycles increases, and its value is less than one, and (ii) the change in the damage factor with the number of saturation–dehydration cycle is continuous. Then,

$$D_1(n)=\frac{G_1(0)-G_1(n)}{G_1(0)}$$

(7)

where $D_1(n)$ is the damage factor of the instantaneous shear modulus of elasticity; $G_1(0)$ is the instantaneous shear modulus of elasticity when the number of saturation–dehydration cycles is zero. $G_1(n)$ is the instantaneous shear modulus of elasticity when the number of saturation–dehydration cycles is n; when n is zero, the rock samples have no damage, $D_1(n)=0$, and as the number of saturation–dehydration cycles increases, $D_1(n)$ increases but is less than one.

Since the level-by-level loading method is used for triaxial creep tests in this study and the creep mechanical properties of the rock samples vary with the deviatoric stress [6], creep parameters corresponding to each deviatoric stress are also different. To determine the final creep parameters of the rock samples, the mean value method is widely adopted in previous research [26,50].

Table 3. Values of creep parameters of argillite after different saturation–dehydration cycles.

Lithology	$\mathit{n}$	$\overline{{\mathit{G}}_1}/\left(\mathbf{GPa}\right)$	D1 (n)	$\overline{{\mathit{G}}_2}/\left(\mathbf{GPa}\right)$	D2 (n)	$\overline{{\mathit{\eta }}_1}/(\mathbf{GPa}·\mathbf{h})$	D3 (n)	$\overline{{\mathit{\eta }}_2}/(\mathbf{GPa}·\mathbf{h})$	D4 (n)
Argillite	0	15.78	0.00	1.73	0.00	10,864.51	0.00	11.24	0.00
	1	13.37	0.15	1.37	0.21	8896.80	0.18	8.97	0.20
	5	6.41	0.59	0.75	0.57	4717.41	0.57	5.74	0.49
	10	2.96	0.81	0.49	0.72	3333.74	0.69	4.05	0.64
	15	0.79	0.94	0.33	0.81	2400.78	0.78	3.39	0.70
	20	0.72	0.95	0.24	0.86	1899.63	0.83	3.02	0.73

Note: $\overline{G_1}$, $\overline{G_2}$, $\overline{{\eta }_1}$ and $\overline{{\eta }_2}$ are respectively the mean values of the creep parameters $G_2$, $G_2$, ${\eta }_1$ and ${\eta }_2$ of the NBVP model; their damage factors are $D_1(n)$, $D_2(n)$, $D_3(n)$ and $D_4(n)$, respectively.

lists the mean values of creep parameters of the NBVP model for argillite under different levels of deviatoric stress after different DW cycles, which were obtained before the rock samples reach the yield strength.

According to Table 3, the mean values of $\overline{G_1}$, $\overline{G_2}$, $\overline{{\eta }_1}$ and $\overline{{\eta }_2}$ of the creep parameters of argillite decrease as the number of saturation–dehydration cycles n increases. Referring to Equation (7), the damage factor of each creep parameter can be determined by fitting the data in Table 3 using the least square method. The equation expressing the relationship between the damage factor of each creep parameter and the number of saturation–dehydration cycles n is written as

$$\{\begin{array}{l}{D}_1(n)=1-\exp (-0.173n)\\ D_2(n)=1-\exp (-0.135n)\\ D_3(n)=1-\exp (-0.123n)\\ D_4(n)=1-\exp (-0.094n)\end{array}$$

(8)

3.4. A Nonlinear Creep Damage Model for Rocks Considering the Effect of Saturation–Dehydration Cycles

Creep damage of rocks is mainly due to the change of creep parameters, such as the shear modulus of elasticity and the viscosity coefficient [34,51]. Therefore, the introduction of a variable parameter (damage factor) into the creep constitutive model to express the damage evolution process provides a new method for studying the creep damage properties of rocks [50,52].

Based on the theory of damage mechanics, Equation (8) is introduced into the NBVP rheological model to establish a damage nonlinear Burgers viscoelastic-plastic (DNBVP) model, which is capable of considering the effect of saturation–dehydration cycles (see Figure 5).

Combining Equations (6) and (8), the rheological equations of the DNBVP model under triaxial creep tests are given as follows:

(1)

When ${\sigma }_1-{\sigma }_3<{\sigma }_{\mathrm{S}}$,

$$\epsilon (t)=\frac{{\sigma }_1-{\sigma }_3}{3G_1(1-D_1(n))}+\frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_2(1-D_2(n))}(1-e^{-\frac{G_2(1-D_2(n))}{{\eta }_1(1-D_3(n))}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2(1-D_4(n))}t$$

(9)

(2)

When ${\sigma }_1-{\sigma }_3\ge {\sigma }_{\mathrm{S}}$ and $t\le t_{\mathrm{P}}$,

$$\begin{array}{l}\epsilon (t)=\frac{{\sigma }_1-{\sigma }_3}{3G_1(1-D_1(n))}+\frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_2(1-D_2(n))}(1-e^{-\frac{G_2(1-D_2(n))}{{\eta }_1(1-D_3(n))}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2(1-D_4(n))}t\\ +(\frac{{\sigma }_1-{\sigma }_3-{\sigma }_{\mathrm{S}}}{2{\eta }_3})t\end{array}$$

(10)

(3)

When ${\sigma }_1-{\sigma }_3\ge {\sigma }_{\mathrm{S}}$ and $t>t_{\mathrm{P}}$,

$$\begin{array}{l}\epsilon (t)=\frac{{\sigma }_1-{\sigma }_3}{3G_1(1-D_1(n))}+\frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_2(1-D_2(n))}(1-e^{-\frac{G_2(1-D_2(n))}{{\eta }_1(1-D_3(n))}t})+\frac{{\sigma }_1-{\sigma }_3}{{\eta }_2(1-D_4(n))}t\\ +(\frac{{\sigma }_1-{\sigma }_3-{\sigma }_{\mathrm{S}}}{2{\eta }_3})\cdot (\exp <\frac{t-t_{\mathrm{P}}}{t_0}{>}^m-1)\end{array}$$

(11)

According to the stress state in rheological tests, relevant parameters of the DNBVP model can be obtained by combining Equations (9)–(11) and using the Origin fitting tool.

4. Discussion

In this section, the DNBVP model verification and parameter sensitivity analysis were conducted. Moreover, the equations necessary for secondary development of this model in FLAC3D were deduced, and the limitations of the DNBVP model were discussed.

4.1. Verification of the DNBVP Model

The rheological curves in Figure 2 were processed using the Boltzmann superposition principle [53,54], and parameters of the DNBVP rheological model were identified according to the creep curves. The specific steps are shown as follows:

(1)

When the deviatoric stress is less than the rock yield strength, the model has four rheological parameters, namely, G₁, G₂, η₁ and η₂. First, according to the rheological test data, let t = 0, and the instantaneous shear modulus G₁ can be determined according to Equation (9). Then, calculate the other three rheological parameters, G₂, η₁ and η₂, by fitting the rheological test data.

(2)

When the deviatoric stress is greater than the rock yield strength, and $t\le t_{\mathrm{P}}$, there are five corresponding rheological parameters, namely, η₁, η₂, η₃, G₁ and G₂, which can be calculated by Equation (10) and the fitting of test curves.

(3)

When the deviatoric stress is greater than the rock yield strength, and $t>t_{\mathrm{P}}$, there are six corresponding rheological parameters, namely, η₁, η₂, η₃, G₁, G₂ and m. As can be seen in Equation (11), the six rheological parameters should be determined in order to solve the equation. The rheological curves are fitted according to Step (2) when the t value is below t_p, and then the corresponding parameters η₁, η₂, G₁ and G₂ can be calculated; parameters η₃ and m can be given by plugging the above parameters into Equation (11).

With the parameters obtained according to three steps above, the resulting data from the creep tests under the fifth level of load are shown in Figure 2—these were used as examples to verify the proposed DNBVP model. Figure 6 illustrate the comparison between the theoretical curves and the test data. The results show that the theoretical curves in Figure 6 can successfully describe the rheological test data of rocks after different saturation–dehydration cycles.

The classic Kelvin–Voigt model [55,56] and Burgers model [29,36,37,38,39,40] were compared with the DNBVP model by fitting the test data with n = 20 in Figure 6, respectively (see Figure 7). The results show that the Kelvin–Voigt model and the Burgers model can only describe the first two stages of creep curves, whereas the proposed DNBVP model can describe all creep stages of rocks, especially the accelerated creep stage. The correlation coefficients of the fitted curves obtained from the three models are 0.85, 0.90 and 0.99, respectively. Reasons for the above difference can be explained as follows: at the first two creep stages, the classic models had little difference from the proposed model in terms of fitted curves. However, after the rock samples enter the accelerated creep stage, the classic Kelvin–Voigt and Burgers models cannot describe the nonlinear properties of the rock samples because they lack corresponding nonlinear parameters, whereas the DNBVP model can accurately reflect the nonlinear properties of the rock samples at this stage because it introduced a nonlinear viscoplastic body.

4.2. Analysis of Sensitivity of the DNBVP Model

With the same number of saturation–dehydration cycles, the most critical parameter in the DNBVP model that reflects the characteristics of the rock creep curves is the nonlinear parameter m. Curves at the accelerated creep stage of rocks are complex [27,33]. To better understand the effect of the change of parameter m on the accelerated creep stage of rocks, we analyzed its sensitivity as follows.

Taking the test data in Figure 7 as an example, the obtained parameters of the DNBVP model are listed in

Table 4. Parameters of the DNBVP model (σ₃ = 19.3 MPa, σ₃ = 4 MPa, n = 20).

P (MPa)	G1 (GPa)	G2 (GPa)	η1 (GPa·h)	η2 (GPa·h)	η3 (GPa·h)	m
15.3	0.77	44.26	81.18	323.82	416.48	NA
15.3 (accelerated stage)	0.77	44.26	81.18	323.82	201.53	0.4

according to the model parameter identification method mentioned in Section 4.1. By keeping the other creep parameters constant in Table 4 and only assigning different values to the nonlinear parameter m, the tendency of change of the creep curve with the variation in m is obtained, as shown in Figure 8.

As can be seen in Figure 8, the value of creep strain and the creep rate strongly depend on the value of m. When m is greater than 0.1, the creep curve exhibits the accelerated creep stage, and the slope of the creep curve increases as the value of m increases, which means the creep deformation rate changes more rapidly. Figure 8 also shows that the time required for the rock samples to enter the accelerated creep stage reduces as the value of m increases. For example, the time required to enter the accelerated creep stage is about 42.5 h when m = 0.1 and is approximately 30 h when m = 3. Moreover, the duration of accelerated creep stage is also shortened as the value of m increases. For example, the duration of the accelerated creep stage is about 20 h when m = 0.4 and is about 1.5 h when m = 3.

From the above analysis, it can be found that, with the same number of saturation–dehydration cycles, the creep rate and the creep curve pattern of the DNBVP model are largely restricted by the parameter m, and the change of m fully reflects the nonlinear accelerated creep property of rocks. This suggests the flexibility of the nonlinear parameter m when describing the DNBVP model.

4.3. Implementation of the DNBVP Model in FLAC3D

FLAC3D is an important tool for analyzing geotechnical engineering issues. To facilitate the practical application of the DNBVP model, FLAC3D was used as the platform in this study to carry out secondary development of the DNBVP model according to the implementation process described in the user manual [57]. The secondary development of the constitutive model using FLAC3D is actually a process in which new stress is obtained according to the stress and strain increment within the last time-step by using the constitutive model equation established [52]. Therefore, it is crucial to derive three-dimensional central difference expressions necessary for secondary development of the constitutive model equation [47]. The three-dimensional central difference expressions of the DNBVP model are derived as follows:

It is specified that $S_{ij}$ is the deviatoric stress; $e_{ij}$ is the total deviatoric stress; $S_{ij}^N$ and $S_{ij}^O$ are the new and old stress deviators within a certain time increment, respectively; and $e_{ij}^N$ and $e_{ij}^O$ are respectively the new and old strain deviators within a certain time increment. $S_{\mathrm{S}}$ is the three-dimensional yield strength.

(1)

When ${(S_{ij})}_0<S_{\mathrm{S}}$, the three-dimensional central difference expression for the new stress deviator is

$$S_{ij}^N=\frac{1}{a}\left[b\Delta e_{ij}+b(A-B)e_{ij}^{2,O}+b{S}_{ij}^O\right]$$

(12)

where

$$a=1+b[\frac{1}{2G_0(1-D_1(n))}+\frac{1}{{\eta }_2(1-D_4(n))}], b=\frac{4A{\eta }_1(1-D_3(n))}{\Delta t}$$

(13)

$$A=1+\frac{G_2(1-D_2(n))\Delta t}{2{\eta }_1(1-D_3(n))}, B=1-\frac{G_2(1-D_2(n))\Delta t}{2{\eta }_1(1-D_3(n))}$$

(14)

$e_{ij}^{2,O}$ is the new strain deviator of the rheological element of the second part in Figure 5.

(2)

When ${(S_{ij})}_0\ge S_{\mathrm{S}}$ and $t\le t_{\mathrm{P}}$, the three-dimensional central difference expression for the new stress deviator is

$$S_{ij}^N=\frac{1}{a}\left[b\Delta e_{ij}-b\Delta e_{ij}^4+b(A-B)e_{ij}^{2,O}+b{S}_{ij}^O\right]\hspace{1em}\hspace{1em}{(S_{ij})}_0\ge S_{\mathrm{S}}, \mathrm{and} t\le t_{\mathrm{P}}$$

(15)

At this moment, the central difference scheme of spherical stress of the model is

$${\sigma }_0^N={\sigma }_0^O+K(\Delta {\epsilon }_{vol}^{}-\Delta {\epsilon }_{vol}^4)$$

(16)

where ${\sigma }_0^N$ and ${\sigma }_0^O$ are the new and old spherical stresses within a certain time increment, respectively. $\Delta {\epsilon }_{vol}^{}$ is the spherical stress increment within the time step $\Delta t$.

(3)

When ${(S_{ij})}_0\ge S_{\mathrm{s}}$ and $t>t_{\mathrm{P}}$, the three-dimensional central difference expression for the new stress deviator is

$$S_{ij}^N=\frac{1}{a}\left[b\Delta e_{ij}-b\Delta e_{ij}^4+b(A-B)e_{ij}^{2,O}+b{S}_{ij}^O\right]{(S_{ij})}_0\ge S_{\mathrm{S}}, \mathrm{and} t>t_{\mathrm{P}}$$

(17)

At this moment, the central difference scheme of spherical stress of the model is

$${\sigma }_0^N={\sigma }_0^O+K\left[\Delta {\epsilon }_{vol}^{}-\frac{<F>}{2{\eta }_3}\Delta t(\exp <\frac{t-t_{\mathrm{p}}}{t_0}{>}^m-1)(\frac{\partial Q}{\partial {\sigma }_{11}}+\frac{\partial Q}{\partial {\sigma }_{22}}+\frac{\partial Q}{\partial {\sigma }_{33}})\right]$$

(18)

In conclusion, the stress–strain relationship of the DNBVP model can be expressed by the three-dimensional central difference expressions for the new stress deviator in Equations (12), (15) and (16). A callable dynamic link library (.dll) is generated through VC++ programming [33], and secondary development of the DNBVP rheological model in FLAC3D can be realized.

4.4. Limitations of the DNBVP Model

In this study, a series of creep tests were performed on argillite after different saturation–dehydration cycles, and the DNBVP model, describing the creep properties after DW cycles, was proposed as well. However, the new model proposed still has some limitations. For example: (i) When the number of saturation–dehydration cycles is greater than 20, the rationality of the model cannot be verified, since no triaxial creep test is performed under such conditions. (ii) The expression of the model is very complex, in which six parameters need to be determined, and some parameters, including the value m, which are independent of the physical property of the sample, do not have definite physical significance. Thus, the variation of these parameters currently cannot be explained by linking the mineralogy, the clay content or other aspects of rock characteristics. (iii) The rock samples used in laboratory tests undergo saturation–dehydration cycles, whereas for rocks that are not in an unsaturated state and whose lithology is not argillite, the rationality of the model needs further verification. Therefore, more samples will be used in our future studies to discuss the above issues to verify and optimize our model. Moreover, the performance of creep tests on samples of other lithology in fully saturated and unsaturated states is also essential, in order to understand the creep behavior of rocks under reservoir water level fluctuations. Finally, it is also a challenge for us in the future to analyze practical engineering cases according to field monitoring data using FLAC3D and the outcomes of secondary development of the DNBVP model, in order to further verify the proposed model.

5. Conclusions

Triaxial creep tests were conducted on argillite after saturation–dehydration cycles and the law in creep tests was analyzed. Furthermore, a nonlinear creep damage model for rocks (the DNBVP model) considering the effect of saturation–dehydration cycles was proposed. The following conclusions can be drawn:

(i)

At the same level of deviatoric stress, the strain after stabilization of the creep curve increases as the number of DW cycles increases. When the number of DW cycle is small, e.g., n = 0, 1 or 5, obvious shear failure zones are observed after sample failure, and conjugate shear fractures occur. Moreover, sample fracture surfaces are approximately X-shaped. However, when n = 10, 15 or 20, the macromorphology of the samples after creep failure is very complex, with many tensile cracks and flaws appearing. It is thus clear that the mechanical properties of the samples are degraded as the number of DW cycles increases, resulting in propagation of many micro-fissures inside the samples.

(ii)

Based on analysis from laboratory triaxial creep tests of the rock samples, a nonlinear damage model considering the effect of saturation–dehydration cycles (the DNBVP model) was established. The key feature of this model is that it introduced a nonlinear viscoplastic body and a damage variable describing saturation–dehydration cycles. The three-dimensional creep equation of the new model was derived and its creep parameters were identified. Comparison between the theoretical curves and the test results shows that the theoretical curves of the DNBVP model can successfully describe the values of rock creep tests after different saturation–dehydration cycles. By comparing classic creep models with the proposed model, it is found that the DNBVP model can accurately reflect the nonlinear characteristics of rocks at the accelerated creep stage.

(iii)

In the DNBVP model, a series of damage variables about the number of cycles n can reflect the effect of saturation–dehydration cycles. The most critical parameter in the DNBVP model is the nonlinear parameter m, which reflects the characteristics of the rock creep curves. Additionally, the creep rate and the creep curve pattern of the DNBVP model strongly depended on the parameter m, and the variation in m fully reflects the nonlinear accelerated creep property of rocks. This suggests the flexibility of the nonlinear parameter m when describing the DNBVP model.

(iv)

To facilitate the implementation of the DNBVP model in FLAC3, three-dimensional central difference expressions necessary for secondary development of the new model were derived in detail. This can provide a reference for the development of other creep models. In addition, the constitutive model developed in this study can improve the geotechnical design of reservoir bank slopes and the control of reservoir bank landslides. It is a challenge for us in the future to analyze practical engineering cases according to field monitoring data using the DNBVP model, in order to further verify the new model we have proposed.