Experimental Investigation and Micromechanics-Based Analytical Modeling of Creep and Relaxation Behaviors of Beishan Granite

Abstract

This paper investigates experimentally and numerically the short- and long-term strength and deformation behaviors of Beishan granite at room temperature. Single-stage creep, relaxation, and conventional triaxial compression tests were performed on cylindrical rock samples. Its typical brittle response is captured and the dependence of peak strength on confining pressure and time-dependent response on deviatoric stress are revealed. For constitutive modeling, a unified micromechanics-based plasticity-damage model is formulated based on the Mori–Tanaka method and the subcritical cracking theory postulate, with the focus on simulating both instantaneous strain and time-dependent deformation process over a broad range of time scales. Its unification is achieved by representing the evolution of damage, which is strongly coupled with plastic deformation induced by frictional sliding along closed cracks, as an internal variable that can be decomposed into instantaneous and time-dependent parts. The performance of the model with analytical predictions is well validated using the experimental results on Beishan granite.

1. Introduction

China has recently developed the program of deep geological disposal of high-level radioactive waste (HLRW), and Beishan granite is considered as the preferred host rock [1,2,3,4]. Similar to such HLRW project, many other deep-ground geotechnical engineering projects such as tunnel building, oil and gas well-boring, and wellbore plugging and abandonment are seeking for more accurate and comprehensive predictions of the long-term behaviors of host-surrounding rocks to ensure the stability and safety of the structures [5,6,7,8]. Thus, the time-dependent mechanical behaviors and deformation characteristics of various potential host rocks have been extensively studied theoretically, experimentally, and numerically [9,10,11,12]. Creep and relaxation tests are two of the main methods of evaluating their time-dependent behaviors, and their function by maintaining either constant strain or constant stress [13,14,15,16]. By parsing out instantaneous deformation and time-dependent deformation, which are two key components of overall deformation in creep and relaxation tests, we expect to describe properly and simulate in a unified way these time-independent and time-dependent mechanical properties.

In the last few decades, extensive experimental research has utilized many different methods to investigate the time-dependent behaviors of the rock [17,18,19,20]. For example, Huang et al. [21] used multi-level loading triaxial-compression creep tests to study the effect of strain history on the deformation of coal samples. Zhang et al. [22] adopted creep and stress–relaxation tests to investigate the effect of asperity degradation on time-dependent behavior, including higher creep and stress–relaxation rates of rock joints. Other researchers [23,24] have investigated the stress–relaxation behaviors of discontinuities under cyclic loading via shear–relaxation tests. Others have proposed advanced techniques, such as digital image correlation, for identifying damage and quantifying full-field strain at various scales. For example, Hou et al. [25] studied the effects of initial damage state on the time-dependent behaviors of sandstone by conducting a creep-test scheme with controlled initial damage.

At the same time, constitutive models of the short- and long-term mechanical properties of multiple types of rock have also been widely studied, both theoretically and numerically [26,27,28,29]. Pietruszczak et al. [30] and Shao et al. [31] presented a mathematical framework for modeling creep in frictional materials, and assumed creep deformation to be associated with progressive evolution of such materials’ micro-structures. Zhou et al. [32] subsequently proposed a unified elastic-plastic and viscoplastic damage model capable of describing instantaneous and time-dependent plastic deformations at various time scales. More importantly, ongoing development of Zhou et al. [32]’s micro-mechanical plastic-damage constitutive model has not only provided more accurate descriptions of the time-dependent behaviors of rock, but has also laid an important foundation for future studies and investigations.

In the present work, driven by the need of a more accurate and comprehensive description of the long term behaviors of the host rock to ensure the integrity and safety of the deep-ground structures, the time-dependent behaviors of Beishan granite are first studied experimentally, through conventional triaxial compression tests and single-stage creep and relaxation tests at various levels of confining pressures and initial deviatoric stress. Then, a unified micromechanics-based plastic damage (UMBPD) analytical model is formulated for simulating instantaneous and time-dependent creep and relaxation behaviors. Specifically, based on Mori–Tanaka homogenization, the representative elementary volume of Beishan granite can be treated as a matrix-inclusion system composed of a solid matrix and randomly oriented penny-shaped microcracks. Second, we will apply a unified plastic-yield criterion to describe time-independent and time-dependent deformation under different levels of high-deviatoric stress. Then, by studying long-term mechanical behaviors, the researchers will ascertain an appropriate time-dependent criterion of Beishan granite for damage due to subcritical cracking, which will be used to simultaneously simulate creep deformation and stress relaxation; and lastly, the simulation results obtained from the constitutive analytical model will be compared against our experimental results.

Throughout the paper, the following notion on tensorial product of any second order tensors $\mathit{A}$ and $\mathit{B}$ will be taken: ${\left(\mathit{A}\otimes \mathit{B}\right)}_{ijkl}=A_{ij}{B}_{kl}$. By using the second order identity tensor $\delta$, the fourth order identity tensor $\mathbb{I}$ and the fourth order hydrostatic projector $\mathbb{J}$ are expressed in the components form as $I_{ijkl}=\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)/2$ and $J_{ijkl}={\delta }_{ij}{\delta }_{kl}/3$, respectively. The fourth order deviatoric projection tensor $\mathbb{K}$ is then introduced.

2. Experimental Study

2.1. Sample Preparations and Testing Procedures

The Beishan granite specimens that are prepared for the triaxial compression test and single-stage creep and relaxation test are quarried from the borehole #BS28 at the Xinchang Site of the underground research laboratory [33]. The dimensions of the cylindrical specimen are 50 mm of diameter and 100 mm of height. The dry bulk density of these specimens are measured by the saturation method and are reported as 2.72 g/cm³. The maximum confining pressure and peak deviatoric stress are taken based on the in situ stress gradient that obtained from the literatures [34,35,36]. Several intermediate stress levels are selected from ambient to the maximum stress to test the creep behaviors’ dependency on the stress levels. By using the rock servo-controlled triaxial testing machine shown in Figure 1, the prepared samples are tested at nine different confining stress levels (0, 1, 2, 5, 10, 15, 25, 30, and 40 MPa) at room temperature.The uniaxial compressive strength (UCS) are obtained from the zero confining stress level and its average value is reported as 128 MPa. The single-stage creep and relaxation tests are carried out at the confining pressure of 15 MPa. The deviatoric stress levels in creep and relaxation tests are different. In the creep test, the three applied deviatoric stress levels are 225.7 MPa, 228 MPa, and 232.5 MPa. In the relaxation test, the three applied deviatoric stress levels are 171 MPa, 199 MPa, and 221 MPa.

The testing procedures start by loading the specimen to a prescribed hydrostatic stress state at a loading rate of 3 MPa/min. The deviatoric stress is then applied in the axial direction in a displacement-controlled mode at a fixed rate of 0.02 mm/min. Once the deviatoric stress arrives the pre-determined value, creep or relaxation test can be started by either maintaining the deviatoric constant and measuring the strain change or maintaining the axial strain constant and measuring the deviatoric change. In both cases, the evolution of strain and stress is measured and recorded by the automatic data collection system equipped by the rock servo-controlled triaxial testing machine.

2.2. Experimental Results

According to the traditional sign convention in geomechanics, the compressive stress is considered here as a positive value. The conventional triaxial compression tests are conducted firstly and the deviatoric stress versus the axial and lateral strains under different confining stresses are shown in Figure 2a. The Young’s modulus can be determined from the linear part of the axial strain and deviatoric stress curve as $E^m=45$ GPa and the Poisson’s ratio can be calculated from the ratio of lateral strain over the axial strain as ${\upsilon }^m=0.15$. Figure 2b presents the curves of the deviatoric stress ${\sigma }_1-{\sigma }_3$ versus volumetric strain ${\epsilon }_v={\epsilon }_1+2{\epsilon }_2$ of Beishan granite under different confining pressures. The transition from compressibility to dilatancy, as well as its dependence on confining pressure are observed on macroscopic scale [37].

The experimental results of single-stage creep tests are shown in Figure 3. A clear trend is observed that creep failure occurs at shorter timeframe when the deviatoric stress is higher. When the deviatoric stress (q) is 225.7 MPa, the creep failure occurs at around 110 h after the deviatoric stress is applied. When the deviatoric stress is increased to 228 MPa, the creep failure occurs at 64 h; as the deviatoric stress increased to 232.5 MPa, the creep failure occurs at around 1.7 h after the deviatoric stress is applied. It is thus concluded that the creep failure time is very sensitive to the levels of the deviatoric stress. When the deviatoric stress is increased by 3%, the creep failure time is reduced by 98%. If we calculate the axial strain rate and use it as the second y axis in the same plot as shown in Figure 3, clear boundaries of three typical creep stages (primary creep stage, steady creep stage, and accelerated creep stage) can be observed in all tested cases with different levels of deviatoric stresses. As summarized in

Table 1. Comparison of single-stage creep test results between Beishan granite.

Deviatoric Stress q	Stage	${\mathit{\epsilon }}_{\mathit{A}}(\%)$	${\mathit{\epsilon }}_{\mathit{L}}(\%)$	${\mathit{\upsilon }}_{\mathit{A}}(\%·{\mathit{h}}^{-1})$	${\mathit{\upsilon }}_{\mathit{L}}(\%·{\mathit{h}}^{-1})$
225.7 MPa	I	1.0$\times {10}^{-3}$	2.1$\times {10}^{-3}$	4.3$\times {10}^{-2}$	2.1$\times {10}^{-2}$
	II	6.9$\times {10}^{-4}$	1.8$\times {10}^{-3}$	7.1$\times {10}^{-4}$	2.0$\times {10}^{-3}$
	III	1.3$\times {10}^{-3}$	4.0$\times {10}^{-3}$	5.9$\times {10}^{-2}$	3.1$\times {10}^{-1}$
228 MPa	I	9.6$\times {10}^{-4}$	1.9$\times {10}^{-3}$	1.8$\times {10}^{-1}$	2.2$\times {10}^{-1}$
	II	6.1$\times {10}^{-4}$	1.5$\times {10}^{-3}$	1.9$\times {10}^{-3}$	6.2$\times {10}^{-3}$
	III	1.2$\times {10}^{-3}$	3.7$\times {10}^{-3}$	4.3$\times {10}^{-1}$	2.1
232.5 MPa	I	8.6$\times {10}^{-4}$	1.1$\times {10}^{-3}$	1.1	1.4
	II	3.5$\times {10}^{-4}$	9.7$\times {10}^{-4}$	1.0$\times {10}^{-1}$	2.8$\times {10}^{-1}$
	III	8.0$\times {10}^{-4}$	2.4$\times {10}^{-3}$	1.7	8.2

, in all three tests, the accelerated creep stage has the highest strain rate in both axial (${\epsilon }_A$) and lateral (${\epsilon }_L$) directions, whereas the steady creep stage has the lowest strain rate. It is noted that, while the lower deviatoric stress case is one where q is equal to 225.7 MPa, the axial strain rates in all three stages are lower than the axial strain rate in the higher deviatoric stress case where q is equal to 232.5 MPa.

The results of the single-stage relaxation tests are shown in Figure 4, where the primary y axis is representing the change of the deviatoric stress and the secondary y axis is representing the change rate of the deviatoric stress. Based on the deviatoric stress changing rate, three relaxation stages, which are rapid relaxation, decelerated relaxation, and steady relaxation, can be observed in all three test cases. It is summarized in

Table 2. Comparison of single-stage relaxation test results between Beishan granite.

q (MPa)	${\mathit{\epsilon }}_0\left({10}^{-3}\right)$	${\mathit{\sigma }}_{\mathit{t}}$ (MPa)	${\mathit{\lambda }}_{12\mathbf{h}}$	$\mathit{\lambda }$	${\mathit{\lambda }}_{12\mathbf{h}}/\mathit{\lambda }$	${\mathit{u}}_{\mathit{A}}$(MPa$·{\mathit{h}}^{-1}$)
171 MPa	3.984	158.1	0.061	0.075	80.4%	0.013
199 MPa	4.584	180	0.067	0.076	88.4%	0.053
221 MPa	5.381	202.6	0.079	0.083	94.6%	0.051

, where q is denoted for initial deviatoric stress, whose value are $q=171$ MPa, $q=199$ MPa and $q=221$ MPa, respectively. ${\epsilon }_0$ represents the initial unchanged axial strain and $u_A$ denotes the relaxation rate of axial stress as shown in in Table 2. Furthermore, ${\sigma }_t$ is the residual deviatoric stress at time t which have been demonstrated in Table 2 Column 3. In this context, the stress relaxation degree $\lambda$ is defined as:

$$\lambda =\frac{q-{\sigma }_t}{q}$$

(1)

According to the residual deviatoric stress at different time, the stress relaxation degree can be divided into different categories ${\lambda }_{12\mathrm{h}}$ and $\lambda ={\lambda }_{36\mathrm{h}}$ in this paper, for reflecting the extent of stress relaxation at 12 h and 36 h [38]. It is concluded from Table 2 that the higher initial deviatoric stress will lead to a higher relaxation degree, and the improvement on the extent of stress relaxation at 12 h.

3. Formulation of the Plastic Damage Model

In this section, a unified micromechanics-based time-dependent damage constitutive model is presented for describing the creep and relaxation mechanical behaviors of Beishan granite. The assumptions of isothermal condition and isotropic material properties are adopted in this paper.

3.1. Free Energy and State Equations

For micro-crack quasi-brittle rocks, the REV (Representative Element Volume) is taken to be a matrix-cracks system, it is assumed composed of an isotropic solid matrix and randomly distributed microcracks. Each micro-crack is geometrically approximated by an oblate spheroid called a penny-shaped crack. When the cracked materials subjected to the a local compressive stress, the friction sliding acting on the surface of closed micro-cracks which induced macroscopic inelastic deformation and micro-crack growth result in the material damage, are two main dissipation processes [39,40].

With the assumption of small deformations, the total macroscopic strain over REV can be divided into two parts: an elastic part ${\epsilon }^e$ attributed to the deformation of isotropic matrix and a plastic part ${\epsilon }^p$ due to the frictional sliding along the surfaces of closed microcracks:

$$\epsilon ={\epsilon }^e+{\epsilon }^p$$

(2)

Apart from the friction-induced plastic strain, another internal damage variable is studied to describe the density parameter of micro-crack which is denoted by d [41]. From the above analysis, one can obtain the free energy for the case of closed cracks [39]:

$$W=\frac{1}{2}\left(\epsilon -{\epsilon }^p\right):{\mathbb{C}}^m:\left(\epsilon -{\epsilon }^p\right)+\frac{1}{2d}{\epsilon }^p:{\mathbb{C}}^b:{\epsilon }^p$$

(3)

Above, the first term on the right-hand side represents the elastic strain energy and the second one represents plastic the strain energy which related to friction sliding along the surface of micro-crack. The fourth tensor ${\mathbb{C}}^b$ can be interpreted as the stiffness modification of damaged matrix–crack system, and it takes the form ${\mathbb{C}}^b=\frac{3{\kappa }^m}{{\alpha }_1}\mathbb{J}+\frac{2{\mu }^m}{{\alpha }_2}\mathbb{K}$, Where ${\kappa }^m$ and ${\mu }^m$ denote the bulk modulus and shear modulus of solid matrix phase. According to the Mori–Tanaka homogenization method [42,43], ${\alpha }_1=\frac{16}{9}\frac{1-{\left({\upsilon }^m\right)}^2}{1-2{\upsilon }^m}$ and ${\alpha }_2=\frac{32}{45}\frac{\left(1-{\upsilon }^m\right)\left(5-{\upsilon }^m\right)}{2-{\upsilon }^m}$ are two constants only dependent on the Poisson’s ratio ${\upsilon }^m$ of the solid matrix.

Under isotropic assumption, ${\mathbb{C}}^m$ is elasticity stiffness tensor of the solid matrix which can be expressed in the following form:

$${\mathbb{C}}^m=3{\kappa }^m\mathbb{J}+2{\mu }^m\mathbb{K}$$

(4)

The macroscopic stress–strain relation can be derived according to the free energy Equation (3):

$$\sigma ={\mathbb{C}}^m:\left(\epsilon -{\epsilon }^p\right)$$

(5)

One then obtains the damage thermodynamic force associated with the internal variables d:

$$Y_d=-\frac{\partial W}{\partial d}=\frac{1}{2d^2}{\epsilon }^p:{\mathbb{C}}^b:{\epsilon }^p$$

(6)

The driving force for plastic strain ${\epsilon }^p$ is also derived:

$${\sigma }^p=-\frac{\partial W}{\partial {\epsilon }^p}=\sigma -\frac{1}{d}{\mathbb{C}}^b:{\epsilon }^p$$

(7)

3.2. Plasticity Criterion and Damage Evolution

A generalized Coulomb-type plastic yield criterion is formulated in terms of the local stress ${\sigma }^p$ acting on crack surfaces.

$$f=\parallel {\mathit{s}}^p\parallel +\eta {\sigma }^p$$

(8)

${\mathit{s}}^p={\sigma }^p:\mathbb{K}$ and ${\sigma }^p=\frac{1}{3}{\sigma }^p:\delta$ which represent the deviatoric and mean parts of the local stress respectively. $\eta$ denotes the friction coefficient of crack’s surfaces.

The evolution of plastic strain can be derived by applying the normality rule:

$${\dot{\epsilon }}^p={\lambda }^p\frac{\partial f}{\partial {\sigma }^p}={\lambda }^p\mathit{D}$$

(9)

with ${\lambda }^p$ is a non-negative plastic multiplier and $\mathit{D}={\displaystyle \frac{{\mathit{s}}^p}{\parallel {\mathit{s}}^p\parallel }}+\frac{1}{3}\eta \delta$ denotes the plastic flow direction.

The damage criterion can be expressed as a function of the damage thermodynamic force:

$$f_d=Y_d-R_d$$

(10)

where $R_d$ represents the damage resistance function which usually depends on the damage state [44]. The following classic form is adopted for $R_d$ in this paper:

$$R_d=r_c\frac{2\xi }{1+{\xi }^2}$$

(11)

with $\xi ={\displaystyle \frac{d}{d_c}}$ which is a dimensionless parameter. Physically, $d_c$ represents the critical value of damage density variable that corresponds to the peak strength. Moreover, $r_c$ is maximal value of damage resistance function $R_d$. The damage increment can be determined by applying the normality rule:

$$\dot{d}={\lambda }^d\frac{\partial f_d}{\partial Y_d}={\lambda }^d$$

(12)

For the case of time-dependent part, the total damage can be divided into two parts: the instantaneous damage $\omega$ induced by external loads and time-dependent damage related to subcritical cracking:

$$d=\omega +\varsigma$$

(13)

where $\omega$ and $\varsigma$ denote the instantaneous damage due to the external mechanical loading, and time-dependent damage which is relative to subcritical cracking due to stress corrosion. The damage criterion (10) is used for instantaneous damage evolution while the time-dependent damage rate can be described in the terms of the distance between the current damage and its equilibrium state $\overline{\varsigma }$. Pietruszczak et al. [30] gave the detailed description of the parameter $\overline{\varsigma }$. The simple linear form is adopted here to describe the time-dependent damage:

$$\dot{\varsigma }=\Gamma \left(\overline{\varsigma }-\varsigma \right)$$

(14)

where $\Gamma$ is the model parameter which controls the rate of accelerated creep and steady relaxation of the presented time-dependent model which will be discussed in Section 4.4.

By setting the initial condition $\varsigma \left(0\right)=0$, the time-dependent damage criterion can be expressed using the Laplacian transformation and convolution operation:

$$\varsigma \left(t\right)={\int }_0^t\overline{\varsigma }\left(\tau \right)\Gamma e^{-\Gamma \left(t-\tau \right)}d\tau$$

(15)

A relevant explicit integration method is taken here [45]:

$${\varsigma }_i={\varsigma }_{i-1}{e}^{-\Gamma \dot{t}}+\frac{1}{2}\Gamma \left({\overline{\varsigma }}_i+{\overline{\varsigma }}_{i-1}\right)e^{-\frac{1}{2}\Gamma \dot{t}}\dot{t}$$

(16)

From above analysis, the parameter $\tilde{\varsigma }$ represents a stationary state of microstructure equilibrium associated with the current damage d. Due to the three-stage feature of creep and relaxation tests on Beishan granite, the following form is adopted:

$$\tilde{\varsigma }=\frac{d^2}{C+d}exp\left(\frac{〈d-d_c〉}{d_c}\right)$$

(17)

where $〈x〉=\frac{x+\left|x\right|}{2}$ is the Macauley bracket and C is a time-dependent parameter whose calibration will de addressed in Section 4.2.

4. Analytical Prediction of Short- and Long-Term Strengths under Conventional Triaxial Compression

The unified micromechanics-based plastic-damage analytical solution under conventional triaxial compression loading is proposed for predicting the time-independent and time-dependent mechanical behaviors including peak strength, creep strain, creep strength, creep rate, relaxation strength, and relaxation rate.

4.1. Short-Term Strength

As a popular loading path, the conventional triaxial compression loading path is conducted here. The plastic strain can be measured from the accumulated multiplier because the flow direction is unchanged in the whole loading process. In this case, the following analytical form is given of the flow direction $\mathit{D}$:

$$\mathit{D}=-\frac{1}{\sqrt{6}}\mathrm{diag}\left(2,-1,-1\right)+\frac{1}{3}\eta \mathrm{diag}\left(1,1,1\right)$$

(18)

Then the accumulated multiplier can be expressed as the following form:

$${\Lambda }^p=\int {\lambda }^p$$

(19)

from which the plastic strain can be rewritten as:

$${\epsilon }^p={\Lambda }^p\mathit{D}$$

(20)

Following the well-established equations, the yield criterion and damage evolution can be reformulated in coincidence with the conventional triaxial compression loading path:

$$\left\{\begin{array}{c}f=\parallel \mathit{s}\parallel +\eta \sigma -\frac{{\Lambda }^p}{d}\chi =0\hfill \\ f_d=\frac{1}{2}{\left(\frac{{\Lambda }^p}{d}\right)}^2\chi -R_d=0\hfill \end{array}\right.$$

(21)

with $\chi =\frac{{\kappa }^m{\eta }^2}{{\alpha }_1}+\frac{2{\mu }^m}{{\alpha }_2}$. By combining these two criteria, we derive the following strength criterion in function of the stresses and damage.

$$f=\parallel \mathit{s}\parallel +\eta \sigma -\sqrt{2R_d\chi }=0$$

(22)

When the loading path of conventional triaxial compression is concerned in particular, the maximum principal stress ${\sigma }_1$ can be expressed in the terms of the confining pressure ${\sigma }_3$ and internal damage variable d.

$${\sigma }_1=\frac{\sqrt{6}+2\eta }{\sqrt{6}-\eta }{\sigma }_3+\frac{3\sqrt{2R_d\chi }}{\sqrt{6}-\eta }$$

(23)

In addition, by combining Equation (5) with Equation (20), the macroscopic strain can be expressed in terms of the macro stress $\sigma$ and the accumulated plastic multiplier ${\Lambda }^p$.

$$\epsilon ={\mathbb{S}}^m:\sigma +{\Lambda }^p\mathit{D}$$

(24)

where ${\mathbb{S}}^m={\left({\mathbb{C}}^m\right)}^{-1}$ is the compliance tensor of the solid matrix phase. With the relation ${\Lambda }^p=d\sqrt{\frac{2R_d}{\chi }}$, the following analytical solution of the short-term strength can be derived [39]:

$$\begin{array}{c}{\epsilon }_1=\frac{1}{E^m}{\sigma }_1-\frac{2{\nu }^m}{E^m}{\sigma }_3-\left(\frac{2}{\sqrt{6}}-\frac{\eta }{3}\right)d\sqrt{\frac{2R_d}{\chi }}\hfill \\ {\epsilon }_2={\epsilon }_3=\frac{1-{\nu }^m}{E^m}{\sigma }_3-\frac{{\nu }^m}{E^m}{\sigma }_1+\left(\frac{1}{\sqrt{6}}+\frac{\eta }{3}\right)d\sqrt{\frac{2R_d}{\chi }}\hfill \end{array}$$

(25)

4.2. Long-Term Strength

In this paper, the crack damage stress ${\sigma }_{cd}$ is used as the long-term strength of Beishan granite which can be defined as the compressibility–dilatancy (C/D) transition stress [18,37,46,47,48]. In our proposed micromechanics-based analytical model, one assumes that there exists a crack critical damage value $d_{cd}$ corresponding to the transition stress ${\sigma }_{cd}$. Then the damage resistance value $R_{cd}$ in terms of the crack critical damage $d_{cd}$ can be derived:

$$R_{cd}=r_c\frac{2{\xi }_{cd}}{1+{\xi }_{cd}^2},\mathrm{with}{\xi }_{cd}=\frac{d_{cd}}{d_c}$$

(26)

Finally, the crack damage stress ${\sigma }_{cd}$ can be expressed in terms of the damage resistance value $R_{cd}$.

$${\sigma }_{cd}=\frac{\sqrt{6}+2\eta }{\sqrt{6}-\eta }{\sigma }_3+\frac{3\sqrt{2R_{cd}\chi }}{\sqrt{6}-\eta }$$

(27)

The short-term peak strength, the crack damage stress and its linear fitting functions under different confining pressures are presented in Figure 5.

It is highlighted that the proposed micromechanics-based analytical time-dependent model contains seven parameters which can be determined from a series of conventional compression tests on Beishan granite. The Young’s modulus $E^m$ and the Poisson’s ratio ${\nu }^m$ are determined from the linear part of the stress–strain curves as shown in Figure 2a. According to Equations (11) and (23), the friction coefficient $\eta$, the critical damage resistance $r_c$ and the crack critical damage $d_{cd}$ can be calibrated by using the linear fitting function of peak stress and transition stress of Beishan granite under different confining pressure, respectively, as shown in Figure 5. The critical damage $d_c$ is related to plastic strain in its peak stress state. In this context, the time-dependent parameter C can be expressed as following form [48]:

$$C=\frac{d_c^2}{d_c-d_{cd}}-d_c$$

(28)

The time-dependent parameter $\Gamma$ can be calibrated by the creep deformation or relaxation stress rate of third stage in triaxial time-dependent test of Beishan granite. The relevant values of model parameters are presented in

Table 3. Parameter’s values of Beishan granite.

Parameters	${\mathit{E}}^{\mathit{m}}$ (MPa)	${\mathit{\nu }}^{\mathit{m}}$	${\mathit{d}}_{\mathit{c}}$	$\mathit{\eta }$	${\mathit{r}}_{\mathit{c}}$	C
Values	45,000	0.15	8	1.71	0.012	2.2

.

A damage-controlled time-dependent model is proposed for the simultaneous account of creep and relaxation mechanical behaviors. The instantaneous damage can be defined by damage evolution $f_d$ while the time-dependent damage can be updated by time-dependent damage criterion $\varsigma \left(t\right)$.

For the case of creep condition, the axial stress ${\sigma }_1$ keeps unchanged. Assuming that during time-dependent stage plastic yield criterion still work, the accumulated multiplier ${\Lambda }^p$ can be expressed by:

$${\Lambda }^p=\frac{\parallel \mathit{s}\parallel +\eta \sigma }{\chi }$$

(29)

According to Equation (25) and plastic yield criterion f, the macroscopic strain-stress relation can be rewritten for both instantaneous and time-dependent conditions:

$$\begin{array}{c}{\epsilon }_1=\frac{1}{E^m}{\sigma }_1-\frac{2{\nu }^m}{E^m}{\sigma }_3+{\mathit{D}}_{11}\frac{{\sigma }_1{\mathit{D}}_{11}+2{\sigma }_3{\mathit{D}}_{33}}{\chi }d\hfill \\ {\epsilon }_2={\epsilon }_3=\frac{1-{\nu }^m}{E^m}{\sigma }_3-\frac{{\nu }^m}{E^m}{\sigma }_1+{\mathit{D}}_{33}\frac{{\sigma }_1{\mathit{D}}_{11}+2{\sigma }_3{\mathit{D}}_{33}}{\chi }d\hfill \end{array}$$

(30)

The second term on the right hand of the equation is plastic strain including instantaneous part and time-dependent part, due to that the damage composed of its instantaneous damage and creep damage. The constants ${\mathit{D}}_{11}$ and ${\mathit{D}}_{33}$ have been presented in Equation (18).

For the case of relaxation condition, the axial strain ${\epsilon }_1$ keeps unchanged. According to Equation (24) and plastic criterion f, the following relations are used for the macroscopic stress ${\sigma }_1$:

$$\left\{\begin{array}{c}f=\parallel \mathit{s}\parallel +\eta \sigma -\frac{{\Lambda }^p}{d}\chi =0\hfill \\ {\epsilon }_1=\frac{{\sigma }_1}{E^m}-\frac{2{\nu }^m}{E^m}{\sigma }_3+{\Lambda }^p{\mathit{D}}_{11}\hfill \end{array}\right.$$

(31)

Accordingly, the accumulated multiplier ${\Lambda }^p$ of relaxation condition can be expressed with the sign convention in geomechanics:

$${\Lambda }^p=\frac{E^m{\epsilon }_1+2{\nu }^m{\sigma }_3-{\sigma }_3\frac{\left(\sqrt{6}+2\eta \right)}{\left(\sqrt{6}-\eta \right)}}{\left[\frac{3}{\left(\sqrt{6}-\eta \right)}\frac{\chi }{d}+E^m\frac{\left(\sqrt{6}-\eta \right)}{3}\right]}$$

(32)

Finally, the macroscopic stress–strain relation of relaxation stage can be expressed as the following form:

$${\sigma }_1=E^m{\epsilon }_1-E^m\frac{\sqrt{6}-\eta }{3}\frac{E^m{\epsilon }_1+2{\nu }^m{\sigma }_3-{\sigma }_3\frac{\left(\sqrt{6}+2\eta \right)}{\left(\sqrt{6}-\eta \right)}}{\left[\frac{3}{\left(\sqrt{6}-\eta \right)}\frac{\chi }{d}+E^m\frac{\left(\sqrt{6}-\eta \right)}{3}\right]}+2{\nu }^m{\sigma }_3$$

(33)

By insertion into Equation (30), one can obtain the macroscopic lateral strain for relaxation stage.

4.3. Computation Algorithm

The strongly nonlinear plastic-damage coupling behaviors can be presented with a strain snap-back mechanical response in brittle rock. A damage-controlled computation algorithm is proved to simulate the short-term mechanical behaviors while a time-controlled method is applied for the time-dependent mechanical responses in this paper. The computation algorithm is detailed as follows.

(1)

Instantaneous stage: The following quantities are assumed to have been completely determined by giving a new damage increment $\Delta d$ for the $(k-1)th$ loading step: ${\sigma }_{k-1}$, ${\epsilon }_{k-1}$, ${\Lambda }_{k-1}^p$.

(2)

Time-dependent stage: Combining the time-dependent damage criterion, the internal variable ${\overline{\varsigma }}_k$, ${\varsigma }_k$, ${\dot{\varsigma }}_k$, $d_k$ can be updated by giving a new time increment $\Delta t$ for the $\left(k\right)th$ loading step.

(3)

Time-dependent stage: We check $f\left({\Lambda }_{k-1}^p,d_k\right)\le 0$ and determine the accumulated plastic multipliers ${\Lambda }_k^p$ if plastic flow occurs.

• The creep stage

  (a).

  ${\sigma }_k={\sigma }_{k-1}$.

  (b).

  Calculate ${\epsilon }_k$ according to Equation (30).

• The relaxation stage

  (a).

  ${\epsilon }_{1,k}={\epsilon }_{1,k-1}$.

  (b).

  Calculate ${\sigma }_{1,k}$, ${\epsilon }_{2,k}$ according to Equations (30) and (33).

4.4. Analytical Prediction

Figure 6 shows the analytical prediction of short-term mechanical response under some classical confining pressures on Beishan granite. Yu et al. [49] had presented the detailed analytical prediction of short-term mechanical behaviors of Beishan granite under all different confining pressures. It is observed that model prediction can correctly capture the typical short-term mechanical characteristics such as the material hardening/softening, confining pressure dependency and volumetric dilatancy, etc. More precisely, the residual strength of model simulation is lager than experimental data due to that the model parameter: friction coefficient $\eta$ keep unchanged in mechanical loading. Zhao et al. [37] gave the residual friction coefficient ${\eta }_r$ to describe the slope of residual strength envelop, from which is considered the degradation of friction coefficient of crack surfaces.

Before carrying out the analytical solutions of time-dependent tests, the influence of the time-dependent parameters $\Gamma$ involved in Equation (14) has been discussed. The model parameter $\Gamma$ has an obvious influence on creep duration time and transitional velocity from the stage I to stage II of relaxation as shown in Figure 7. The tests were carried out in a thermally isolated small room with a constant temperature of $20\pm 0.5$ °C maintained during the experiment because the time-dependent tests are highly influenced by the temperature. Therefore, the time-dependent parameter $\Gamma$ only influenced by microstructure of material, which is related to the loading history [48]. In this paper, different values of $\Gamma$ are adopted for achieving satisfactory simulations due to the effect of loading history.

Analytical prediction of compression time-dependent tests are illustrated and compared with experimental data on Beishan granite under confining pressure $p_c=15$ MPa, as shown in Figure 8 and Figure 9. It is can be seen that the analytical prediction is in good agreement with experimental data on creep strain, relaxation stress, its duration time, and long-term strength of creep and relaxation experiments.

5. Conclusions

By conducting the triaxial compression experiments at room temperature, the important mechanical parameters of Beishan granite were characterized and summarized. More importantly, through the time-dependent tests including single-stage creep and relaxation tests, the the axial/lateral strain–stress curve, time–strain, and time–stress data are also obtained. These experimental results show that the long-term behavior of the Beishan granite is highly dependent on the level of the deviatoric stress. Our experimental data shows that even if a small increment of the deviatoric stress will significantly reduce the creep and relaxation duration time.

Furthermore, the present work is uniquely proposing a unified micromechanics-based plastic-damage (UMBPD) analytical model which can describe both instantaneous strain and time-dependent strain over a broad range of time scales in both creep and relaxation tests. By well benchmarking with the experiment results from long-term tests under different deviatoric stresses on Beishan granite, UMBPD is well demonstrated to be accurate and reliable. The influence of loading history on rock creep deformation and relaxation stress is also captured by the time-dependent parameters $\Gamma$. The future work includes developing a unified micromechanics-based anisotropic time-dependent plastic-damage model. To fulfill this goal, more experiments on various rocks, especially the ones with significant anisotropic properties, should be tested to provide sufficient experimental data on the micro-mechanical constitutive relations under a long-term test time frame.