1 INTRODUCTION
The phenomenon of rock creep refers to the temporal variation in the strain of rock materials under continuous stress circumstances (Sun, 2007). Empirical evidence from engineering studies clearly shows that the creep properties of rocks have a considerable impact on the long-term stability of tunnels. To address the unequal geographical distribution of water resources in Yunnan Province, the Central Yunnan water diversion project is a crucial initiative of immense importance to the province's sustainable economic and social development. The Haidong water conveyance tunnel in this project segment traverses a region with significant topographical variations. The depth of the tunnel mainly falls between 200 and 400 m, with a maximum depth reaching 500 m. After the tunnel excavation is finished, the primary factor responsible for the prolonged distortion of the surrounding rock is the creep phenomena of the rock mass. The diabase in the weak and broken zone of the tunnel has the characteristics of strong weathering, high burial depth, and prominent creep. The tunnel has a long cycle of three-dimensional plastic deformation, and the monitoring data of the elevated arch does not converge after the lining, which adversely affects the safety of tunnel construction and operation. Therefore, it is necessary to study its creep mechanical properties and construct a constitutive model that can reflect the creep deformation and failure law of diabase.
In recent years, many scholars have carried out many creep test studies under different conditions according to actual working conditions (Chen et al., 2024; Xu et al., 2019). In terms of research on the creep characteristics of rocks, Chen, Wang et al. (2023) conducted triaxial unconfined creep tests on deeply buried granite tunnels to investigate the effect of confining pressure on the creep characteristics of rocks. Tao et al. (2020) investigated the relationship between the creep characteristics of thin carbonaceous slate in the Muzailing Tunnel and the angle of slate stratification and water content through uniaxial compression creep tests. Chen, Zhang et al. (2023) investigated the creep characteristics of sandstone in the Yuezhishan Tunnel through a triaxial graded loading‒unloading creep test. Zhou et al. (2022) investigated the crehenep properties of siltstone subjected to triaxial unloading. Currently, there is a lot of research on the creep of road tunnels, but less attention has been paid to the creep characteristics of the surrounding rock of water conveyance tunnels, especially large-diameter tunnels, under deep burial conditions. A fundamental and challenging aspect of the theoretical research of rock creep mechanics is the establishment of a rock creep constitutive model. Scholars have employed numerous theories and methods to characterize the nonlinear curve observed in rock creeps (Liu et al., 2017; Zhao et al., 2017). Zhou et al. (2011) introduced a new creep constitutive model that utilizes the theory of time-based fractional derivative. Yang et al. (2014) constructed a new nonlinear viscoelastic-plastic creep model with a creep threshold and long-term strength. Zhang et al. (2022) employed the fractional Abel dashpot as the viscous element in constructing a temperature-dependent creep constitutive model for salt rock. Given this, a large number of results have been obtained on the creep mechanical properties and constitutive models of rocks. However, the linear element is used in most creep models to describe the attenuation stage of rocks, which cannot be applied to all rocks. The creep constitutive equation has the characteristics of complex form and high parameter inversion time, which brings inconvenience to the identification and calculation of the function.
In view of this, with the diabase in the Haidong water conveyance tunnel of China as the research object, this study performed triaxial graded loading creep tests under different stress levels. The creep mechanical characteristics of diabase were examined by analyzing the test curves obtained at various confining pressures. Based on the theory of fractional derivative, a nonconstant coefficient Abel dashpot considering damage was established, and the simplified equation form was used to facilitate parameter identification and calculation. Combining a nonlinear Kelvin model, a classical element model, and a simplified Abel dashpot that considers damage effects, a constitutive model was established to describe the diabase creep process. The suitability of the model was confirmed by comparing the creep test curves of single-axis and triaxial compression with the theoretical curves obtained by identifying the model parameters.
2 DIABASE CREEP CHARACTERISTIC TEST
2.1 Sample preparation and testing equipment
The sample was taken from the soft and broken zone upstream of Hole 3 of the Haidong water conveyance tunnel. The sample is a deeply weathered diabase of grey-green and grey-black color, composed mainly of pyroxene and basic plagioclase. The structure of the rock mass is shattered and loose. The collected samples were uniformly processed into standard cylinders with a diameter of 50 mm and a height of 100 mm according to the International Society for Rock Mechanics standards.
Samples were subjected to uniaxial and triaxial compression creep testing utilizing the five-channel deep soft rock rheological experimental system at the China University of Mining and Technology-Beijing. The equipment system consists of three essential components: the primary engine, the hydraulic system, and the measurement and control system. This system can monitor and evaluate the complete creep process in real time and with exceptional accuracy.
2.2 Test scheme
According to the actual working conditions of the tunnel, the confining pressure for the conventional triaxial compression test was set at 0, 5, 10, and 15 MPa. The stress gradient for loading in the creep test was chosen using the diabase uniaxial compressive strength value as a guide. The compression creep test was designed to apply constant confining pressure and gradually increase the axial pressure. The initial loading stress for the uniaxial creep test was intended to be 8 MPa, increasing in 2 MPa increments. Each stress level was stabilized for 24 h until the sample was entirely destroyed. The test was set up with a starting loading stress equal to 45% of the triaxial compressive strength value. Each subsequent loading stress was then raised by 10% of the compressive strength value. Each stress level was stabilized for 24 h until the sample was utterly destroyed. Due to the existence of some discrete specimens, the sample grading exceeded the theoretical maximum grades during the actual axial stress loading at the confining pressures of 5 and 15 MPa. The stress loading values for each level were recorded, as shown in Table 1.
Table 1. Loading values of creep tests under step loading. (MPa)
| Deviatoric stress |
Axial stress |
|
|
|
|
|
| First level |
8 |
35.8 |
46.6 |
53.9 |
| Second level |
10 |
44.9 |
59.2 |
68.0 |
| Third level |
12 |
54.0 |
71.2 |
82.2 |
| Fourth level |
14 |
63.0 |
84.3 |
96.4 |
| Fifth level |
16 |
72.1 |
96.0 |
110.5 |
| Sixth level |
18 |
81.2 |
107.7 |
124.7 |
| Seventh level |
20 |
90.2 |
- |
138.9 |
| Eighth level |
22 |
98.0 |
- |
153.0 |
| Ninth level |
24 |
- |
- |
- |
| Tenth level |
26 |
- |
- |
- |
2.3 Analysis of the experimental results
2.3.1 Analysis of the overall creep pattern
The entire diabase creep deformation process was recorded in Figure 1 at 0, 5, 10, and 15 MPa. The sample undergoes conventional compression failure at the last stage under a confining pressure of 15 MPa. Deformations of all samples exhibit positive values for compressive stress and negative values for tensile stress. Based on the creep test curve, it is evident that the sample experiences substantial elastic deformation when subjected to constant axial stress at all levels. The sample then progresses to the attenuation creep stage, where the creep rate steadily declines but remains nonzero. Throughout this phase, the microcracks in the sample gradually undergo expansion. A steady-state creep phase is initiated in the sample, characterized by a generally constant creep rate and a linear rise in strain as time progresses. When the axial stress reaches a specific threshold, the internal microcracks in the sample undergo rapid expansion and penetration. This causes the sample to enter an accelerated creep stage, characterized by a rapidly increasing creep rate, ultimately leading to failure.
Creep strains under different confining conditions. (a)
σ
3 = 0 MPa, (b)
σ
3 = 5 MPa, (c)
σ
3 = 10 MPa, and (d)
σ
3 = 15 MPa.
With confining pressures of 0, 5, and 10 MPa, the sample successively passes through the attenuation, isothermal, and accelerated creep stages at the final stress level. The total time for the entire process was 415, 55, and 25 s, respectively. It can be seen that the whole process of creep damage to diabase lasts a relatively short time. As the confining pressure increases, the total duration gradually becomes shorter with each increase in loading stress. According to the graded loading creep curve, the penultimate stress level was defined as the stress threshold of the sample. The sample was severely creep-damaged at an axial stress level above the stress threshold. According to the test results, the stress thresholds of the sample under four constant confining pressures of 0, 5, 10, and 15 MPa are 24.0, 90.2, 96.0, and 153.0 MPa, respectively, indicating that the ability of the sample to resist deformation can be improved by increasing the confining pressure.
The creep deformation of diabase with confining pressures of 5 and 10 MPa was recorded, as shown in Table 2. The analysis table shows that the shift in axial creep strain is minor at low-stress levels and that the axial creep strain increases with increasing axial stress at high-stress levels. The lateral creep strain always increases with increasing axial stress. Axial creep strain is significantly greater than lateral creep at the initial stress level. As the stress level continues to increase, the lateral creep strain of the sample increases rapidly, and the axial creep strain gradually becomes smaller than the lateral strain. It can be concluded that diabase is mainly deformed by axial creep at low-stress levels and that lateral creep is more evident than axial creep at high-stress levels.
Table 2. Creep deformation of diabase.
| Deviatoric stress |
|
|
| Axial strain (%) |
Lateral strain (%) |
Axial strain (%) |
Lateral strain (%) |
| First level |
0.01260 |
0.0014 |
0.0158 |
0.0010 |
| Second level |
0.01084 |
0.0061 |
0.0100 |
0.0096 |
| Third level |
0.00980 |
0.0080 |
0.0124 |
0.0108 |
| Fourth level |
0.01100 |
0.0099 |
0.0133 |
0.0136 |
| Fifth level |
0.01080 |
0.0122 |
0.0261 |
0.0286 |
| Sixth level |
0.01170 |
0.0157 |
- |
- |
| Seventh level |
0.02200 |
0.0380 |
- |
- |
2.3.2 Long-term strength of the diabase
Long-term strength is the strength value at which a rock remains stable under long-term loading. An isochronous stress–strain curve approach is used to determine the long-term strength of diabase. The stress value corresponding to the starting point of divergence of the curve cluster is the long-term strength of the rock. The isochronous stress–strain curve of diabase under different confining pressures is shown in Figure 2. The long-term strengths of diabase under confining pressures of 0, 5, 10, and 15 MPa were found to be 18.5, 69.6, 77.5, and 111.6 MPa, respectively, and the ratios of these strengths to the creep failure strengths at the corresponding confining pressures were 71.2%, 71.0%, 72.0%, and 72.9%, respectively.
Isochronous deviatoric stress–strain curves under different confining pressures. (a)
σ
3 = 0 MPa, (b)
σ
3 = 5 MPa, (c)
σ
3 = 10 MPa, and (d)
σ
3 = 15 MPa.
3 NONLINEAR CREEP DAMAGE MODEL
3.1 Fractional derivative creep element
Fractional calculus can describe higher-order derivatives and integrals and can be defined in multiple ways. In this paper, the Riemann–Liouville definition of fractional calculus is adopted, and the expression is defined by (Zhou et al.,
2022),
(1)
where
is a function of time,
is the fractional order,
is a variable, and
is the Gamma function, and the expression is defined by
(2)
According to the definition of the fractional derivative, the constitutive relation of the Abel dashpot is given by
(3)
During the creep test, the axial pressure is set to a constant value. By integrating both sides of the equation according to the Riemann–Liouville fractional calculus theory, we obtain
(4)
where
is the viscosity coefficient.
During the accelerated creep stage, damage accumulates rapidly within the rock, and its viscosity coefficient undergoes a continuous shift in line with the progression of the creep process. To describe the deterioration of the viscosity coefficient, it is necessary to adopt a damage variable defined by the time of load application. The viscosity coefficient considering damage can be expressed as
(5)
where
is the damage variable of the rock.
It is assumed that the damage evolution of diabase during creep follows a negative exponential form. The expression can be rewritten as (Zhou et al.,
2012)
(6)
where
a is a coefficient related to the properties of diabase.
By combining Equations (
3), (
5), and (
6), the constitutive relation of the Abel dashpot is given by
(7)
where
is the time-dependent viscosity coefficient.
By performing the Laplace transform on Equation (
7), the constitutive relation can be rewritten as
(8)
As can be seen from Equation (8), it contains an infinite series of terms, which makes it inconvenient to identify and calculate the function. Therefore, by simplifying the function form, the inversion calculation time of the parameters in the formula is greatly reduced.
Taylor's Formula is given by
(9)
Substituting
into the Gamma function, and combining
with Equation (
9), Equation (
8) can be rewritten as
(10)
3.2 Establishment of a one-dimensional nonlinear viscoelastic-plastic damage model
By analyzing the characteristics of the creep curve of diabase under different stress conditions, it can be seen that when the stress level is less than the threshold for sample failure, diabase will produce a certain amount of elastic deformation at the moment of loading under various axial stresses. The constitutive relationship of elasticity can be described using an elastic body as a model. Subsequently, the sample enters the attenuation creep stage, where the creep rate gradually decreases to zero but does not reach zero. The creep strain curve has prominent nonlinear characteristics at this stage. The nonlinear Kelvin model is constructed to describe the attenuation creep stage of the rock. The sample then progresses into a steady-state creep phase characterized by a linear increase in strain with time. In this stage, the stress–strain element model adopts a viscous body. When the stress level exceeds the sample stress threshold, a nonlinear element consisting of a switch and a nonconstant coefficient Abel dashpot in parallel is constructed based on the theory of fractional calculus. This nonlinear element is used to describe the accelerated creep damage stage of the sample. By connecting the above four elements in series, an intrinsic model that can describe the whole process of creep damage of diabase is established, as shown in Figure 3.
Sketch of the nonlinear viscoplasto-plastic creep model.
The expression for total strain is given by
(11)
where
,
,
, and
are the strains of the elastic body, the viscoelastic body, the viscous body, and the viscoplastic body, respectively.
The constitutive relationship of the elastic body is given by
(12)
where
is the elastic modulus.
The creep strain curve has prominent nonlinear characteristics at the attenuation creep stage. Since the parameters in the Kelvin model are constants, this stage cannot be accurately described. A nonlinear Kelvin model that more accurately describes the creep stage of rock attenuation is established by treating the viscosity coefficient in the Kelvin model as a function of the power of time (Zhang & Wang,
2020). The constitutive relation for the nonlinear Kelvin model is given by
(13)
where
and
are the elastic modulus and viscosity coefficient of the viscoelastic body, respectively, and
is a constant.
Solving differential Equation (
11), we can get
(14)
The constitutive relation of the viscous body is given by
(15)
where
is the viscosity coefficient of the Newtonian body.
If
<
, the viscoplastic Abel dashpot fails to initiate the process.
(16)
where
is the stress threshold of the sample.
If
≥
, the constitutive relation of viscoplastic Abel dashpot can be rewritten as
(17)
where
is the viscosity coefficient of the viscoplastic body.
Synthetically, the one-dimensional nonlinear viscoelastic-plastic damage model is given by
(18)
3.3 Establishment of a three-dimensional nonlinear viscoelastic-plastic damage model
Practical engineering usually involves rocks subjected to intricate conditions characterized by numerous stresses. Based on rheological theory, the three-dimensional nonlinear viscoelastic-plastic damage model was obtained by analogy from the one-dimensional creep model.
Under three-dimensional stress conditions, the stress tensor at any point inside a rock can be decomposed into a spherical stress tensor
and a deviatoric stress tensor
, and the expressions are
(19)
where
is the Kronecker function and
is the stress tensor.
The strain tensor can be decomposed into a spherical strain tensor
and a deviatoric strain tensor
, and the expressions are
(20)
The expression for total strain under the 3D stress state is given by
(21)
where
,
,
, and
represent the elastic strain tensor, viscoelastic strain tensor, viscous strain tensor, and viscoplastic strain tensor, respectively.
The generalized Hooke's law is given by
(22)
where
and
, respectively, represent the bulk modulus and shear modulus, and the expressions are
,
.
The constitutive relationship of the elastic body is given by
(23)
The constitutive relation of the viscoelastic body under the 3D stress state is given by
(24)
The constitutive relation of the viscous body under the 3D stress state is given by
(25)
The constitutive relation of the viscoplastic body under the 3D stress state is given by (Jiang et al.,
2013)
(26)
where
is the yield function and
is the initial value and generally taken as 1 (Perzyna,
1966).
is the plastic potential function; according to the correlation flow rule,
can be taken (Nazary Moghadam et al.,
2013).
is the switch function, and the expression is
(27)
where
is the power function,
,
is the material constant and generally taken as 1 (Abu al-Rub et al.,
2013).
The expression for the rock yield function is (Qi et al.,
2012)
(28)
where
is the second invariant of stress deviation.
Under the conditions of a conventional triaxial compression creep test, the lateral confining pressure is the same, that is,
. According to the principle of superposition, the nonlinear viscoelastic-plastic constitutive relation of rocks under the three-dimensional stress states is given by
(29)
4 MODEL PARAMETER IDENTIFICATION AND VERIFICATION
The validity and precision of the developed model are confirmed by doing a fit analysis and parameter identification on the creep test curve and then comparing the resulting curve with the actual test curvature. The creep test curves of the diabase under different test conditions were processed according to the Chen loading method to obtain the axial strain versus time curves under each loading condition. Utilizing the quasi-Newton technique and the broad global optimization approach of the 1stOpt mathematical analysis software, the parameters of the creep model were determined from the results of the creep experiment. Other parameters are determined by the application of a creep constitutive model (Zhang et al., 2011).
If the axial stress level does not exceed the stress threshold and the steady-state creep strain rate is not zero, the second formula in Equations (18) and (29) is utilized to fit and determine the model parameters; if the axial stress level exceeds the stress threshold, the third formulae in Equations (18) and (29) are employed to fit and calculate the model parameters. The model parameters fitted under confining pressures of 0 and 10 MPa are presented in Tables 3 and 4, respectively. The square of the correlation coefficient for the parameter fitting at each stress level with confining pressures of 0 and 10 MPa exceeds 0.98, indicating a high degree of fit.
Table 3. Uniaxial compression creep model parameters.
|
(MPa) |
|
|
|
|
|
|
|
|
|
| 8 |
6.941 |
71.472 |
280.926 |
0.551 |
1.330 × 104 |
- |
- |
- |
0.997 |
| 10 |
7.332 |
63.021 |
184.213 |
0.518 |
4.300 × 103 |
- |
- |
- |
0.999 |
| 12 |
7.767 |
44.283 |
172.774 |
0.542 |
6.990 × 105 |
- |
- |
- |
0.997 |
| 14 |
8.194 |
45.892 |
159.461 |
0.537 |
1.380 × 107 |
- |
- |
- |
0.998 |
| 16 |
8.708 |
51.625 |
126.076 |
0.536 |
1.310 × 106 |
- |
- |
- |
0.999 |
| 18 |
9.110 |
47.849 |
122.196 |
0.515 |
9.340 × 104 |
- |
- |
- |
0.998 |
| 20 |
9.492 |
42.952 |
121.528 |
0.484 |
3.160 × 106 |
- |
- |
- |
0.998 |
| 22 |
9.905 |
41.594 |
100.436 |
0.519 |
1.230 × 107 |
- |
- |
- |
0.998 |
| 24 |
10.296 |
37.229 |
84.278 |
0.509 |
2.540 × 105 |
- |
- |
- |
0.998 |
| 26 |
11.052 |
225.345 |
0.061 |
1.153 |
9.563 |
108.461 |
9.739 |
504.364 |
0.989 |
Table 4. Parameters of the compression creep model with a confining pressure of 10 MPa.
|
(MPa) |
|
|
|
|
|
|
|
|
|
|
| 46.6 |
28.495 |
15.402 |
114.719 |
113.223 |
0.49 |
1.110 × 104 |
- |
- |
- |
0.9980 |
| 59.2 |
29.735 |
16.192 |
95.757 |
146.015 |
0.45 |
1.080 × 104 |
- |
- |
- |
0.9986 |
| 71.2 |
31.109 |
16.535 |
69.252 |
163.363 |
0.40 |
4.500 × 104 |
- |
- |
- |
0.9991 |
| 84.3 |
33.557 |
16.771 |
61.116 |
164.338 |
0.39 |
3.710 × 104 |
- |
- |
- |
0.9995 |
| 96.0 |
36.067 |
16.596 |
47.647 |
125.681 |
0.42 |
2.120 × 104 |
- |
- |
- |
0.9997 |
| 107.7 |
45.604 |
15.484 |
162.837 |
2.947 |
0.57 |
2.035 |
25.510 |
5.44 |
5271.9 |
0.9986 |
The comparison of the test curve and the fitted curve under different confining pressure conditions is shown in Figures 4 and 5. As evident from the diagram, the fitting and test curves exhibit a high degree of agreement in both the one-dimensional and three-dimensional stress conditions. The nonlinear viscoelastic-plastic model established in this paper can accurately describe the characteristics of attenuation and isothermal and accelerated creep failure of diabase under different confining pressures.
Comparison between experimental and fitted curves for stresses below the threshold value. (a)
σ
3 = 0 MPa, (b)
σ
3 = 5 MPa, (c)
σ
3 = 10 MPa, and (d)
σ
3 = 15 MPa.
Comparison of experimental and fitted curves for stresses above the threshold value. (a)
σ
3 = 0 MPa, (b)
σ
3 = 5 MPa, and (c)
σ
3 = 10 MPa.
5 CONCLUSION
This study analyzed the creep characteristics of diabase by carrying out compression creep tests under different confining pressures. A nonlinear viscoelastic-plastic damage creep model is established based on fractional derivatives, which can describe the whole process of creep failure of diabase. The main conclusions are as follows:
1.
When the stress level does not exceed the diabase's stress threshold, the sample exhibits a decay creep stage and a steady-state creep stage with a nonzero creep rate. If the stress level exceeds the threshold, the sample enters an accelerated creep stage with a rapidly increasing creep rate and finally fails. Diabase is dominated by axial creep deformation at low-stress levels, and lateral creep strain is more significant at high-stress levels.
2.
Based on the theory of fractional derivative, a nonconstant coefficient Abel dashpot considering damage was established, and the simplified equation form was used to facilitate parameter identification and calculation. Combining a nonlinear Kelvin model, a classical element model, and a simplified Abel dashpot that considers damage effects, a constitutive model that can describe the diabase creep process was established.
3.
According to the results of the creep experiment, the parameters of the creep model were identified using 1stOpt mathematical analysis software. The triaxial compression creep test curve was compared with the theoretical curve, indicating that this creep damage model can accurately describe the sample's attenuation creep, isothermal creep, and accelerated creep characteristics.
ACKNOWLEDGMENTS
National Natural Science Foundation of China, Grant/Award Number: 42377154.
CONFLICT OF INTEREST STATEMENT
The authors declare no conflict of interest.
Biography
Zhigang Tao is a professor and doctoral supervisor at China University of Mining and Technology-Beijing. His research focuses on large deformation support materials and disaster control in geotechnical engineering. He has presided over or participated in 35 scientific research projects and won 10 scientific research awards, including China Patent Gold Award, Liaoning Provincial Science and Technology Progress Award and China Society for Rock Mechanics and Rock Engineering Science and Technology Progress Award. He has published 59 SCI- or EI-indexed papers as the the first author or corresponding author, with a total of 1378 citations (including 7 Chinese science and technology excellent journals, 17 international high-level journals, 3 high cited papers in ESI and 2 hot papers). And he has published three monographs, authorized 21 Chinese invention patents, 15 utility model patents and 9 software works.
