Study on Creep Characteristics and Constitutive Relation of Fractured Rock Mass

Abstract

In order to effectively describe the whole creep process of fractured rock mass, triaxial unloading creep tests were carried out on prefractured coal samples using constant axial pressure and graded unloading confining pressure, and the axial and lateral creep laws of fractured coal samples with different dip angles were studied. Combined with the characteristics of creep curve and based on Kachanov’s creep damage theory, the damage variable is introduced into the constitutive relation and creep equation, and the evolution equation of damage variable with time in the whole creep process is derived. At the same time, a new method to calculate the initial damage is proposed. The elastoplastic body with damage variable is connected with the Burgers model in series. Meanwhile, lade criterion and switch element are introduced into the creep model to establish a new fracture damage creep model. The one-dimensional and three-dimensional damage creep equations are derived. The damage creep equation is obtained according to the superposition principle. A simple and feasible method for parameter identification of the model is given based on the characteristics of creep curve. The applicability of the model is verified by comparing the creep test curve of fractured coal sample with the theoretical curve. The results show that the two models are in good agreement. The model can not only accurately reflect the nonlinear characteristics of creep curves in the attenuation and isokinetic stages but also describe the accelerated creep characteristics of fractured rocks.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Rock burst refers to the local high stress concentration in coal (rock) caused by the mining of working face and the tunneling of roadway, which makes the local elastic energy in coal (rock) accumulate stably in time and transfer to weak surface through easy propagation path, resulting in continuous generation of weak structure in coal and rock, causing local fracture, damage, instability, and sliding of coal (rock). The phenomenon that the elastic potential energy stored in coal (rock) is released instantaneously after some external disturbance [1–5]. This phenomenon can only last for a few seconds to tens of seconds, which will produce severe vibration, accompanied by a large number of coal (rock) thrown out, seriously threatening the safety of coal production [6–8]. With the increase of coal mining depth, when the rock stratum is relatively soft and the coal seam with strong risk is impacted, the strength, failure range, and destructive force of shock wave of rock burst have obvious change characteristics. Rock burst will also occur under the condition of no disturbance. The sudden release of this type of energy is because of creep-type rock burst [9–11]. Therefore, it is of great significance for the prevention and control of coal body impact disaster that the creep-type impact risk can be predicted by quickly judging the creep deformation of coal body on-site. The existing prevention and control methods of rock burst mainly include the following: (1) energy releasing drilling method; (2) high-pressure water injection in coal seam; (3) deep hole controlled blasting; and (4) hydraulic slotting pressure relief method [12–15]. Li Chao et al. used the method of cutting coal body with high-pressure water jet to prevent rock burst and monitored the test results with KBD5 electromagnetic radiation monitor. It shows that the effect of cutting on stress unloading of shallow coal is the best, and the stress is redistributed around the seam, so that the high stress is transferred from shallow to deep. The existing prevention methods of rock burst are to release the accumulated energy in coal and rock mass by releasing confining pressure, so as to reduce the stress concentration in coal and rock mass and prevent rock burst [16–20].

The rheological property of rock is an important factor affecting the long-term stability of rock mass engineering. There are many factors affecting the creep process of rock, and one of the important factors is the influence of rock damage on creep. In the process of rock creep, the internal damage cracks of rock are produced continuously, which makes the bearing capacity of rock decrease continuously, making the rock to enter the accelerated creep stage. In addition, under the coupling effect of engineering disturbance and rock damage fracture penetration and fusion, roof falling, rock burst, and other disasters are likely to occur. Therefore, it is necessary and significant to introduce damage body to simulate the whole creep process curve of rock. At present, the existing rock creep model simulation has made great achievements, but there are still many problems in the simulation of accelerated creep process [21–25]. Therefore, this paper establishes a creep model to simulate accelerated creep process by introducing a damage body.

In order to simulate the accelerated creep process of rock, scholars at home and abroad have done a lot of research on nonlinear rheological theory and achieved good results: one is that three different material element models (such as Hooke body, Saint Venant body, and Newton body) can be combined in series and parallel to establish a new creep constitutive model, and a lot of rock creep can be established by this method of constitutive models, such as Maxwell model, Kelvin model, and Xiyuan model; Zhao Yanlin constructed a new rheological model, which was composed of components in series, and could well describe the creep law under the state of graded stress; Xia Caichu et al. proposed a unified rheological mechanical model including four basic rheological mechanical forms: viscoelastic, viscoplastic, viscous, and viscoelastic-plastic, and gave the identification method of the model.

Although there are many methods to construct traditional components, they still have some limitations. They can only reflect the characteristics of viscoelastic-plasticity and cannot simulate the process of accelerated creep well. Therefore, some nonlinear models are proposed. The others are based on the damage theory because the basic parameters of creep can be changed, so it can be realized by using the damage factor. The new component is combined with the traditional element in series, and the creep constitutive equation which can simulate the whole process of creep is established. Murakammi and Ohno proposed a damage model for jointed rock mass with multiple sets of planar fractures based on the description of the geometric size of the fracture system; Kawamoto et al. further hypothesized and improved the abovementioned model for jointed rock mass; Zhu Weishen and Qi Yinping introduced a second-order tensor to describe the damage of rock mass based on the geometric damage principle and the definition of the abovementioned damage variables. Based on Kachanov’s damage theory, Wu Zhulin et al. simulated the initial damage factor in series with Burgers element model and achieved good results. However, most of the existing creep models are based on the combination of original or improved components. Either there are many components, or the formula is complex, or cannot be solved. It can be seen that the simulation method of accelerated creep stage still has some shortcomings, which needs to be improved. Therefore, the damage factor is introduced into the creep model to simulate the whole creep process. There are two methods of damage factor: the first is to define the damage factor from the angle of energy damage; the second is to define the damage factor based on the effective bearing area of the structure in terms of geometric damage. However, the nonlinear characteristics of rock creep are not only in the accelerated creep stage but also in the whole creep process. However, most of the creep damage models focus on the accelerated creep process after the stress level exceeds the long-term strength of the rock, and there is little research on the nonlinear characteristics of the early stable creep. In this paper, the initial damage is defined on the basis of existing damage mechanics, and the evolution equation of damage variable is derived according to the Kachanov damage rate. At the same time, based on the lade failure criterion, the elastoplastic damage element is introduced to simulate the plastic damage body and then is combined with the Burgers model in series, a constitutive model with few parameters, easy to determine, uncomplicated components, and able to simulate accelerated creep is constructed. The research in this paper can be used for reference for the stability of underground engineering in the future, especially for the stability of support in deep mining and the prevention of rock burst [26–29].

2. Materials and Methods

2.1. Specimen Preparation

The coal sample was taken from Wang JiaLing coal mine in Shanxi Province. The coal is hard and anthracite. In order to avoid the influence of anisotropy on the test results, rock samples were processed along the same direction when preparing coal samples. At the same time, in order to ensure the authenticity and comparability of the test results, the appearance of the processed coal samples was carefully observed before the test to determine that there were no obvious weak planes such as joints and cracks so as to ensure that there is no obvious difference between the test rock samples in macroscopic. According to the requirements of the International Society for Rock Mechanics (ISRM), the rock samples are processed into standard cylindrical rock samples with a diameter of 50 mm × 100 mm.

The precast fractured rock sample is shown in Figure 1, in which the fracture length is 2a and the fracture dip angle is α. The high-speed electric cutting machine is used to process three-dimensional cracks. The cutting wheel is 0.3 mm thick ultrathin diamond saw blade. The crack thickness is 0.3–0.5 mm, and the crack is filled with gypsum. Change the fracture dip angle α (30°, 45°, and 60°), fracture length 2a = 20 mm, and fracture spacing 2b = 30 mm.

[figures omitted; refer to PDF]

2.2. Test System and Loading Process

The triaxial unloading creep test was performed using the ROCK 600-50 triaxial test system manufactured by Topouan, France (Figure 2). The function and application scope of the equipment are as follows: it can be used for mechanical experiment of rock, rock-soil, and concrete under normal condition and temperature bar with a sample size of φ 50 mm × 100 mm, rheological test of rock mechanics, permeability test of rock mechanics, and mechanical coupling test of rock under conventional condition and temperature condition. The main technical indexes of the equipment are as follows: the maximum axis is 380 MPa, the maximum confining pressure is 60 MPa, the maximum pore pressure is 40 MPa; the axial range of deformation measurement is 0–10 mm, and the circumferential range is 0–4 mm; and the measurement resolution is 0.00001 mm, which fully meets the working requirements of long-time high load.

[figure omitted; refer to PDF]

In the test, constant axial pressure is used to release confining pressure by stages, and load control mode is adopted. Firstly, according to the hydrostatic pressure condition σ₁ = σ₂ = σ₃ = 10 MPa, confining pressure is applied to 10 MPa at the loading rate of 0.5 MPa/s, and then the axial load is continuously applied to 75% of the predicted ultimate bearing capacity under the confining pressure, which exceeds the uniaxial compressive strength of the coal sample and then maintains the principal stress difference σ₁–σ₃. Under the same condition, the confining pressure is gradually reduced by 0.5 MPa/s, and the unloading is carried out step by step from 10 MPa to 8 MPa, 6 MPa, 4 MPa, and 2 MPa; each stage is maintained for 24 h until the coal sample is completely destroyed. The specific control parameters are shown in Table 1.

Table 1

Testing schemes of grading unloading pressure in the triaxial creep test.

Number	Fracture characteristics	σ3 (MPa)	σ3 (MPa)
			Level 1	Level 2	Level 3	Level 4	Level 5
A00	Complete coal sample	45	10	8	6	4	2
A30	30° angle	40	10	8	6	4	2
A45	45° angle	35	10	8	6	4	2
A60	60° angle	33	10	8	6	4	2

2.3. Test Results

In this section, the unloading creep test results of prefractured coal samples obtained by using the triaxial rheometer are analyzed, and the rheological deformation and failure law of coal samples under different fracture dip angles are studied. In order to ensure the reliability of the test data, the creep test is divided into three groups, four specimens are tested in each group, and the creep time of each stage is 24 h. The test results are shown in Table 2.

Table 2

Triaxial unloading creep test data of prefabricated cracked coal samples.

Number	σ1 (MPa)	σ2 (MPa)	εi1 (%)	εi3 (%)	${\epsilon }_{iv}$ (%)	εc1 (%)	εc3 (%)	${\epsilon }_{cv}$ (%)
A00	45	10	1.682	0.615	0.452	0.080	0.011	0.058
		8	3.444	1.248	0.948	0.085	0.016	0.054
		6	5.211	1.900	1.411	0.089	0.027	0.035
		4	6.982	2.605	1.773	0.095	0.032	0.031
		2	8.759	3.363	2.032	—	—	—
A30	40	10	1.793	0.644	0.506	0.102	0.020	0.062
		8	3.690	1.361	0.968	0.111	0.022	0.067
		6	5.595	2.151	1.293	0.121	0.029	0.064
		4	7.510	3.043	1.424	0.137	0.034	0.069
		2	9.469	4.034	1.401	—	—	—
A45	35	10	1.970	0.699	0.572	0.119	0.019	0.081
		8	4.129	1.466	1.198	0.131	0.038	0.055
		6	6.299	2.292	1.716	0.142	0.083	−0.023
		4	8.487	3.204	2.079	—	—	—
A60	33	10	2.084	0.757	0.570	0.188	0.012	0.164
		8	4.286	1.534	1.217	0.196	0.061	0.073
		6	6.501	2.380	1.741	0.209	0.096	0.016
		4	8.751	3.309	2.134	—	—	—

σ₁ and σ₃ are the axial pressure and confining pressure; ε_i1, ε_i3, and ${\epsilon }_{iv}$ are the axial, lateral, and volumetric instantaneous strain increments, respectively; ε_c1, ε_c3, and ${\epsilon }_{cv}$ are the axial, lateral, and volumetric creep strain increments, respectively.

The creep law of the prefractured coal sample is obtained through the creep test of unloading confining pressure step by step under constant axial pressure, as shown in Figure 3.

[figures omitted; refer to PDF]

The creep curve obtained under the condition of step unloading creep is ladder shape, which cannot be directly used and needs to be transformed. According to Chen’s loading treatment method, the whole curve is transformed into a creep curve under various loads, as shown in Figures 4–7.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

It can be seen from Figures 3–7 that different types of coal samples show significant creep deformation under different stress levels and axial deformation and lateral deformation are composed of instantaneous deformation and creep deformation caused by unloading.

2.3.1. Analysis on Creep Deformation Law of Complete Coal Sample

Taking A00 sample as an example, the instantaneous deformation and creep deformation of A00 sample increased with the increase of deviator stress (Figures 3 and 4). The larger the deviator stress is, the greater the increase range of creep deformation is. When the confining pressure is 10 MPa, obvious instantaneous strain is produced. With the decrease of confining pressure, the axial and lateral instantaneous strain values increase gradually. When the confining pressure is reduced from 10 MPa to 2 MPa, the increment of axial instantaneous strain is about 1.77 × 10⁻², and the lateral increment is about 3.5 × 10⁻⁴. Since the second stage unloading, the axial and lateral creep deformation of each stress level is significantly increased, and the axial strain increase is about 0.5 × 10⁻². When the confining pressure is 2 MPa, the axial and lateral creep deformation increases significantly, and the accelerated creep phenomenon appears. Within 1 hour, the axial strain increased by 4.05 × 10⁻⁴ in the direction almost parallel to the axial strain axis, and the lateral deformation increased by 0.2 × 10⁻⁴. The results show that the internal bearing capacity of the sample is weakening with time, and the through fracture surface under the level of fracture confining pressure leads to the failure of the sample.

2.3.2. Analysis on Creep Deformation Law of Fractured Coal Sample

According to the triaxial unloading creep test results of coal samples with different fracture dip angles, the influence of fractures on the creep deformation of coal samples with different inclination angles is analyzed, and the creep mechanical properties of coal samples with different fracture dip angles are discussed. In Figures 3(b)–3(d), the creep test curves of coal samples with fracture dip angles of 30°, 45°, and 60° under a confining pressure of 10 MPa and stress level of 75% triaxial compressive strength are given, respectively; unloading creep curve clusters of fractured coal samples are given in Figures 5–7, and the number above the curve indicates the value of confining pressure at all levels.

According to Figures 3(b)–3(d), coal samples with different fracture dip angles still show three typical stages of rock mass rheology in the process of creep deformation, which are primary attenuation creep stage, steady-state uniform flow stage, and accelerated creep stage, and are destroyed in the final nonlinear accelerated creep stage. Compared with the intact coal, the fractured coal samples have significant creep deformation under high confining pressure. Taking A45 sample as an example, when the axial pressure is 33 MPa and the confining pressure is 10 MPa, the instantaneous strain of the specimen is larger. With the decrease of confining pressure, the increment of each stage of axial instantaneous strain is 2.16 × 10⁻², 2.17 × 10⁻², and 2.19 × 10⁻², and the increment of lateral instantaneous strain is 2.4 × 10⁻⁴, 2.6 × 10⁻⁴, and 2.9 × 10⁻⁴. Compared with the intact coal sample, the deformation increases significantly. When the confining pressure is 10 MPa, the axial strain increased by 1188.8 × 10⁻⁶, and the lateral strain increased by 5.88 × 10⁻⁶ after 24 hours at a constant stress level; when the confining pressure drops to 8 MPa, the axial strain increased by 1306 × 10⁻⁶, and the lateral strain increased by 11.84 × 10⁻⁶ after 24 hours at the constant stress level. When the confining pressure is reduced to 6 MPa, the axial strain and lateral strain increased by 1422 × 10⁻⁶, and 25.94 × 10⁻⁶, respectively, after 24 hours at constant stress level. When the confining pressure is 4 MPa, the creep deformation characteristics of the samples are very obvious. There are three typical creep stages in the whole stage, and the creep deformation characteristics of nonlinear acceleration are presented. In a short time, the axial and lateral creep deformation increased rapidly by 739.6 × 10⁻⁶ and 184.2 × 10⁻⁶.

3. Establishment of Rock Damage Creep Model

3.1. Damage Analysis

In the process of rock formation, there are a lot of microdefects (joints, holes, etc.) in the process of rock formation, so the initial damage is not zero before the creep test. Damage variable is used to measure the damage state of materials, and it is an important link to establish the damage mechanics model of materials. It can be established from microscopic (number of voids, length, area, volume, etc.) and macroscopic (elastic constant, yield stress, tensile strength, etc.). Considering that the angle and length of some joints and fissures in rock mass can be observed in practical engineering, this paper describes the initial damage of rock mass from the geometric damage principle.

The fracture geometry distribution of the rock sample is shown in Figure 8. The damage calculation method of the fractured rock sample is as follows: firstly, each fracture is projected on the plane perpendicular to x₁, and the projection area S₃₀ is calculated; secondly, the distance L₃₀ between the center normal vectors of the two projection planes is calculated; thirdly, considering the influence of the fracture spacing on the mechanical properties of the rock sample, the spacing L₀ with the minimum influence of the fracture on the mechanical properties of the coal sample is calculated; and fourthly, the damage coefficient λ is introduced. The specific calculation formula is$$\begin{array}{c}\text{(1)}& D_{30}=\frac{V_{30}}{V}\times \lambda =\frac{S_{30}\times a}{V}\times \frac{L_{30}}{L_0}\\ =\frac{564.882\times 0.5}{\pi \times {25}^2\times 100}\times \frac{41.51}{38.08}\\ =1.568\times {10}^{-3},\end{array}$$where D₃₀ is the damage value of 30° fractured coal sample; V₃₀ is the projected volume of crack in the x₁ direction; λ is the damage coefficient, λ = L₃₀/L₀; and V is the volume of the rock sample.

[figures omitted; refer to PDF]

According to the nonlinear characteristics of rock creep stage, many researches on the rock creep model have been carried out in recent years. There are four methods to establish the nonlinear creep model:

 

① According to the indoor or field creep test and measured data, the empirical model can be obtained by using the creep formula or empirical formula fitting regression, but its universality is poor

 

② The parameters of the linear element are treated unsteady to make it a nonlinear element, and then the modified nonlinear creep model is obtained

 

③ Based on the damage theory of creep and damage quantity, the damage model is established

 

④ The creep model is obtained considering the external environment, such as the nonlinear creep model coupled with temperature and chemical conditions

Under the action of stress, the accumulation and evolution of internal crack damage is the essential reason for the final deformation and failure of rock. However, the fractured rock mass itself has some damage before the test, and then the creep test is carried out to further expand the damage. Therefore, the creep damage of fractured rock mass does not only exist in the accelerated creep stage but also exist throughout the whole creep process. The research idea of this section is based on the time-dependent evolution mechanism of unloading creep damage of fractured rock mass, introducing the damage factor and combining with the existing creep element model, to establish a creep model which can describe the nonlinear deformation characteristics of each stage of unloading creep process of fractured rock mass, including attenuation creep, steady-state creep, and accelerated creep.

According to the Kachanov damage rate [30], the evolution equation of damage variable is as follows:$$\begin{array}{}\text{(2)}& \frac{\text{d}D}{\text{d}t}=C{\frac{\sigma }{1-D}}^n,\end{array}$$where C and n are the material parameters. By integrating the above formula, we can get$$\begin{array}{}\text{(3)}& C{\sigma }^{n}t=K-\frac{{1-D}^{n+1}}{n+1}.\end{array}$$

The fractured rock mass itself contains different degrees of initial cracks, and the rock mass contains initial damage during loading. Therefore, the initial damage is not 0 at the moment of loading.$$\begin{array}{}\text{(4)}& D_{t=0}=D_0\ne 0.\end{array}$$

By introducing equation (4) into (3), we can get$$\begin{array}{}\text{(5)}& K=\frac{{1-D_0}^{n+1}}{n+1}.\end{array}$$

By introducing equation (5) into (3), we can get$$\begin{array}{}\text{(6)}& C{\sigma }^{n}t=\frac{{1-D_0}^{n+1}}{n+1}-\frac{{1-D}^{n+1}}{n+1}.\end{array}$$

The evolution law of damage variable of fractured rock mass in the whole creep process is as follows:$$\begin{array}{}\text{(7)}& D=1-{{1-D_0}^{n+1}-n+1C{\sigma }^{n}t}^{1/n+1}.\end{array}$$

The Burgers model is composed of the Maxwell model and the Kelvin model in series. The Kelvin body describes the viscoelastic properties of materials and reflects the attenuation creep of rocks; the Maxwell volume describes the elastic and viscous characteristics of the material, reflecting the instantaneous deformation of rock in the loading stage and the constant rate creep deformation in the steady-state creep stage. The Burgers model has the characteristics of elastic viscoelastic viscosity, which can better describe the attenuation and steady-state creep process of general rock mass under low stress level. In this paper, the Lade criterion considering plastic deformation is introduced into the Burgers model, and the damage variable D and switching element are introduced to construct the damage creep model which can describe the instantaneous elastic strain, viscoelastic strain, viscous strain, and nonlinear elastic-plastic strain of rock. The mechanical model is shown in Figure 9.

[figure omitted; refer to PDF]

3.2. Damage Creep Model of Fractured Rock Mass

According to the series parallel relationship of the model, we can know that$$\begin{array}{}\text{(8)}& \begin{array}{c}\sigma ={\sigma }_e={\sigma }_{ve}={\sigma }_v={\sigma }_{vp},\\ \epsilon ={\epsilon }_e+{\epsilon }_{ve}+{\epsilon }_v+{\epsilon }_{vp},\\ {\sigma }_e=E_1{\epsilon }_e,\\ {\sigma }_{ve}=E_2{\epsilon }_{ve}+{\eta }_2\stackrel{.}{{\epsilon }_{ve}},\\ {\sigma }_v={\eta }_1\stackrel{.}{{\epsilon }_v},\\ {\sigma }_{vp}=E_{3}1-D{\epsilon }_{ve}+{\sigma }_s,\end{array}\end{array}$$where ${\sigma }_e,{\sigma }_{ve},{\sigma }_v,\hspace{0.17em}\text{and}\hspace{0.17em}{\sigma }_{vp}$ are the stress of the elastic body, Kelvin body, viscous body, and damaged elastoplastic body, respectively, and ${\epsilon }_e,{\epsilon }_{ve},{\epsilon }_v,\hspace{0.17em}\text{and}\hspace{0.17em}{\epsilon }_{vp}$ are the strain of the elastic body, Kelvin body, viscous body, and damaged elastic-plastic body, respectively.

In order to compare with the test results, the following reasonable assumptions are made firstly:

 

① The rock is isotropic, and the damage is assumed to be consistent in all directions

 

② The damage time and creep time are consistent in the accelerated creep stage

The evolution law of the one-dimensional Burgers model based on damage variable in the whole creep process is as follows:$$\begin{array}{}\text{(9)}& \epsilon =\frac{\sigma }{E_1}+\frac{\sigma }{{\eta }_1}t+\frac{\sigma }{E_2}1-\exp -\frac{E_2}{{\eta }_2}t+\frac{\sigma }{E_D}{{1-D_0}^{n+1}-n+1C{\sigma }^{n}t}^{-1/n+1}.\end{array}$$

3.3. Three-Dimensional Creep Equation of Damage Creep Model of Fractured Rock Mass

In the three-dimensional stress state, the stress tensor inside the material can be decomposed into the spherical stress tensor σ_m and the deviatoric stress tensor S_ij:$$\begin{array}{c}\text{(10)}& \begin{array}{c}{\sigma }_m=\frac{1}{3}{\sigma }_1+{\sigma }_2+{\sigma }_3=\frac{1}{3}{\sigma }_{kk},\\ S_{ij}={\sigma }_{ij}-{\delta }_{ij}{\sigma }_m={\sigma }_{ij}-\frac{1}{3}{\delta }_{ij}{\sigma }_{kk},\end{array}{\sigma }_{ij}=S_{ij}+{\delta }_{ij}{\sigma }_m,\end{array}$$where δ_ij is the Kronecker function. It is generally believed that the spherical stress tensor can only change the volume of the material, but not its shape, while the partial stress tensor can only cause the shape change. Similarly, the strain tensor can be decomposed into spherical strain tensor and deviatoric strain tensor, as shown in the following formula:$$\begin{array}{c}\text{(11)}& \begin{array}{c}{\epsilon }_m=\frac{1}{3}{\epsilon }_1+{\epsilon }_2+{\epsilon }_3=\frac{1}{3}{\epsilon }_{kk},\\ e_{ij}={\epsilon }_{ij}-{\delta }_{ij}{\epsilon }_m={\epsilon }_{ij}-\frac{1}{3}{\delta }_{ij}{\epsilon }_{kk},\end{array}{\epsilon }_{ij}=e_{ij}+{\delta }_{ij}{\epsilon }_m.\end{array}$$

For the Hooke body under three-dimensional stress state,$$\begin{array}{}\text{(12)}& \begin{array}{c}{\sigma }_m=3K{\epsilon }_m,\\ S_{ij}=2G{e}_{ij},\end{array}\end{array}$$where K is the bulk modulus and G is the shear modulus. Under the condition of isotropic material, it is assumed that the elastic strain is caused by the stress sphere tensor and the creep is caused by the stress deviator. The creep equation of the material under three-dimensional stress state is$$\begin{array}{}\text{(13)}& \begin{array}{c}{\epsilon }_m=\frac{{\sigma }_m}{3K},\\ e_{ij}=\frac{S_{ij}}{2G_1}+\frac{S_{ij}}{2{\eta }_1}t+\frac{S_{ij}}{2G_2}1-\exp -\frac{G_2}{{\eta }_2}t,\end{array}\end{array}$$where K is the bulk modulus and G₁ and G₂ are the shear modulus. K, G₁, G₂, η₁, and η₂ are all the model parameters, which are determined by fitting the test results:$$\begin{array}{}\text{(14)}& \begin{array}{c}K=\frac{E_1}{31-2\mu },\\ G_1=\frac{E_1}{21+\mu },\\ G_2=\frac{E_2}{21+\mu },\end{array}\end{array}$$where μ is Poisson’s ratio. Therefore, the parameters K, G₁, and G₂ that need to be fitted can be transformed into E₁ and E₂:

The three-dimensional creep equation of the damaged plastic body is as follows:$$\begin{array}{}\text{(15)}& {\epsilon }_{ij}^{vp}=\frac{1}{G_3}{{1-D_0}^{n+1}-n+1C{\sigma }^{n}t}^{-1/n+1}\Phi \frac{F}{F_0}\frac{\partial Q}{\partial {\sigma }_{ij}},\end{array}$$where $\mathrm{\Phi }F/F_0=\begin{array}{cc}\mathrm{\Phi }F/F_0& F\ge 0\\ 0& F<0\end{array}$, $•$ is a switch function, F is the yield function of rock; F₀ is the initial value of the rock yield function, which is generally taken as 1 for the convenience of calculation; $\mathrm{\Phi }F/F_0$ is the power function, $\mathrm{\Phi }F/F_0={F/F_0}^m$; M is the experimental constant, generally take 1; and Q is the plastic potential function. For the convenience of calculation, the relevant flow rule is adopted, F = Q.

The Mohr–Coulomb yield criterion and the von Mises yield criterion are widely used. However, the Mohr–Coulomb criterion cannot consider the effect of intermediate principal stress on rock yield, and von Mises yield criterion ignores the influence of spherical stress on rock creep characteristics. In this paper, the yield criterion is adopted, which has the following advantages:

(1)

The results show that the three curves are smooth on the meridian plane, which is conducive to the numerical calculation and the determination of the direction of plastic strain increment

(2)

The S-D effect of different tensile and compressive strengths of geotechnical materials is considered

(3)

The material parameters are few and easy to be tested

(4)

It can be widely used in sandy soil, rock-soil, concrete, and over consolidated clay

(5)

It reflects the nonlinear relationship between the shear expansion yield curve and the hydrostatic pressure

(6)

The effect of hydrostatic pressure P on the yield of geotechnical materials and the nonlinear characteristics of yield and failure are considered

The yield criterion is expressed as$$\begin{array}{}\text{(16)}& F=\frac{I_1^3}{I_3}-27{\frac{I_1}{p_a}}^m-k,\end{array}$$where I₁ and I₃ are the first and third invariants of the stress tensor; p_a is the atmospheric pressure; and m and k are the stress level parameters.

For conventional creep tests,$$\begin{array}{}\text{(17)}& \begin{array}{c}{\sigma }_1>{\sigma }_2={\sigma }_3,\\ {\sigma }_m=\frac{{\sigma }_1+{\sigma }_2+{\sigma }_3}{3}=\frac{{\sigma }_1+2{\sigma }_3}{3},\\ S_{11}={\sigma }_1-\frac{{\sigma }_1+2{\sigma }_3}{3}=\frac{2{\sigma }_1-{\sigma }_3}{3},\\ I_1={\sigma }_1+{\sigma }_2+{\sigma }_3={\sigma }_1+2{\sigma }_3,\\ I_3={\sigma }_1{\sigma }_2{\sigma }_3={\sigma }_1{\sigma }_3^2,\\ \frac{\partial F}{\partial {\sigma }_{11}}=\frac{1}{p_a}\frac{2{\sigma }_1}{{\sigma }_3^2}+\frac{6}{{\sigma }_3}-\frac{8{\sigma }_3}{{\sigma }_1^2}-27.\end{array}\end{array}$$

The creep equation of the nonlinear creep model of fractured rock mass under three-dimensional stress state is obtained by using simultaneous equations (15)–(17):$$\begin{array}{}\text{(18)}& \epsilon =\begin{array}{c}\frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_1}+\frac{{\sigma }_1-{\sigma }_3}{3{\eta }_1}t+\frac{{\sigma }_1-{\sigma }_3}{3G_2}1-\exp -\frac{G_2}{{\eta }_2}t\\ +\frac{{\sigma }_1-{\sigma }_3}{3G_3}{{1-D_0}^{n+1}-n+1C{\frac{2{\sigma }_1-{\sigma }_3}{3}}^{n}t}^{-1/n+1},\hspace{1em}\sigma \le {\sigma }_s,\\ \frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_1}+\frac{{\sigma }_1-{\sigma }_3}{3{\eta }_1}t+\frac{{\sigma }_1-{\sigma }_3}{3G_2}1-\exp -\frac{G_2}{{\eta }_2}t\\ +\frac{{\sigma }_1-{\sigma }_3-{\sigma }_s}{G_3}\frac{1}{p_a}\frac{2{\sigma }_1}{{\sigma }_3^2}+\frac{6}{{\sigma }_3}-\frac{8{\sigma }_3}{{\sigma }_1^2}-27{{1-D_0}^{n+1}-n+1C{\frac{2{\sigma }_1-{\sigma }_3}{3}}^{n}t}^{-1/n+1},\hspace{1em}\sigma >{\sigma }_s.\end{array}\end{array}$$

4. Results and Discussion

4.1. Parameter Identification

In the process of model parameter identification, the selection of the mathematical model identification method is very important. At present, the main methods to determine rock creep model parameters are as follows: least squares method, regression inversion method, and displacement back analysis method. The least squares fitting curve is widely used to determine the relevant parameters, but the disadvantage is that the unreasonable selection of the initial value will lead to divergence of results or slow convergence speed. In this paper, the parameters of the rock creep model are identified by the mathematical optimization analysis software 1stOpt. The software does not need the user to provide appropriate initial parameters, which overcomes the shortcomings of ordinary least squares method. Levenberg–Marquardt + general global optimization algorithm is adopted to find the optimal model parameters. Before fitting, equation (19) is complex, which is not conducive to the preparation of the user-defined function, so it is necessary to simplify the coefficient. The results are as follows:$$\begin{array}{}\text{(19)}& \epsilon =\begin{array}{c}A+Bt+E1-\exp -Dt+F{{1-D_0}^{n+1}-n+1C{\frac{2{\sigma }_1-{\sigma }_3}{3}}^{n}t}^{-1/n+1},\hspace{1em}\sigma \le {\sigma }_s,\\ A+Bt+E1-\exp -Dt+G{{1-D_0}^{n+1}-n+1C{\frac{2{\sigma }_1-{\sigma }_3}{3}}^{n}t}^{-1/n+1},\hspace{1em}\sigma >{\sigma }_s,\end{array}\end{array}$$where$$\begin{array}{c}\text{(20)}& A=\frac{{\sigma }_1+2{\sigma }_3}{9K}+\frac{{\sigma }_1-{\sigma }_3}{3G_1},B=\frac{{\sigma }_1-{\sigma }_3}{3{\eta }_1},\\ E=\frac{{\sigma }_1-{\sigma }_3}{3G_2},\\ D=\frac{G_2}{{\eta }_2},\\ F=\frac{{\sigma }_1-{\sigma }_3}{3G_3},\\ G=\frac{{\sigma }_1-{\sigma }_3-{\sigma }_s}{G_3}\frac{1}{p_a}\frac{2{\sigma }_1}{{\sigma }_3^2}+\frac{6}{{\sigma }_3}-\frac{8{\sigma }_3}{{\sigma }_1^2}-27.\end{array}$$

According to the fitting results, the coefficients A∼G are obtained, and then the model parameters are converted.

4.2. Model Validation

In order to verify the effect of the creep damage model established in this paper on the whole creep stage, the traditional Burgers model and the F-D model considering the initial damage are used to fit the test curve at the same time. The fitting results are shown in Table 3, and the comparison of fitting curves is shown in Figure 10.

Table 3

Identification results of creep model parameters.

Model	K	G1 (MPa)	η1 (MPa·h)	G2 (MPa)	η2 (MPa·h)	G3 (MPa)	C	n	R2
Burgers	1405.173	1282.984	970486.977	24829.455	69348.095	—	—	—	0.997
	690.772	630.705	943793.638	25684.545	67545.498	—	—	—	0.994
	460.658	420.601	928555.343	26136.347	63703.257	—	—	—	0.991
	347.036	316.859	842878.145	25978.658	55542.856	—	—	—	0.990
	265.569	249.237	798621.687	24036.263	49851.892	—	—	—	0.969
F-D	982.525	897.283	635384.554	16248.256	45454.447	191574.563	7.713 × 10−9	0.055	0.998
	485.155	443.064	630126.484	17148.671	45070.881	104343.923	7.929 × 10−9	0.015	0.996
	325.072	296.869	631538.602	17776.032	43292.905	68274.135	8.471 × 10−9	0.003	0.994
	245.983	224.642	582988.970	17966.628	38380.009	51212.810	2.469 × 10−9	0.292	0.993
	198.108	180.921	704302.959	15649.894	27820.801	55882.069	1.040 × 10−3	1.093	0.994

[figures omitted; refer to PDF]

Combined with Table 3 and Figure 10, it can be seen that both Burgers model and F-D model can accurately describe viscoelastic creep of fractured rock mass. In terms of fitting effect, adding damage elements to the Burgers model in series and introducing the Lade criterion can make the simulation effect of the improved model (F-D) more accurate and better reflect the creep deformation characteristics of fractured rock mass. Moreover, the correlation coefficient R² of model identification is more than 0.99, and the fitting accuracy of the model meets the requirements. The experimental curve and the theoretical curve are in good agreement; the theoretical curve of the model is basically consistent with the test curve. The F-D model can not only accurately describe the creep characteristics of attenuation and constant velocity stage but also describe the accelerated creep characteristics of fractured rock mass. The deformation characteristics of three typical stages of rock rheology can be well described, which verifies the applicability of the model.

5. Conclusions

Rock burst is one of the major disasters in coal mines because of its great destructiveness. Because of the mining depth, creep-type rock burst occurs frequently, and the number of coal mines and coal seams is increasing. This kind of rock burst disaster has seriously restricted the safety production of coal mines in China. The prediction and prevention of this kind of rock burst is different from the conventional type. If the mining boundary is complex, the self-weight stress field is large, and the residual tectonic stress field exists, and the superposition of these multiple factors will make the coal (rock) body produce “strong” creep failure. If disturbed by the dynamic stress field of the working face, it is easy to produce large-scale and serious rock burst. The existing prevention methods of rock burst are to release the energy accumulated in coal and rock mass by releasing confining pressure, so as to reduce the stress concentration in coal and rock mass and prevent rock burst phenomenon. Based on this, triaxial unloading confining pressure creep tests are carried out on fractured coal samples with different angles, and the fracture damage creep model is established. The conclusions are as follows:

(1)

Under the same stress state, the different fracture dip angle significantly affects the creep deformation of coal samples. With the decrease of confining pressure, the axial instantaneous strain and creep strain increase with the increase of fracture dip angle, each curve is basically straight line, and the deformation increase amount of the same type coal sample under each stress level is basically the same. The order of axial instantaneous strain and creep strain from large to small are 60° fracture dip angle, 45° fracture dip angle, 30° fracture dip angle, and complete coal sample. Under the same stress state, with the decrease of confining pressure, the order of lateral instantaneous strain from large to small is 60° fracture angle, 45° fracture angle, 30° fracture angle, and complete coal sample; compared with complete coal sample, due to the existence of prefabricated cracks, the bearing capacity of fractured coal sample in creep process decreases, and the lateral expansion deformation increases. The larger the fracture angle, the greater the lateral strain.

(2)

Based on the geometric damage principle, a second-order tensor is introduced to describe the initial damage of rock mass. According to Kachanov’s creep damage theory, the damage variable with initial damage is introduced into the constitutive relation and creep equation, and the evolution equation of damage variable with time in the whole creep process is derived.

(3)

The elastoplastic body with damage variable is connected with the Burgers model in series. Meanwhile, the Lade criterion and switch element are introduced into the creep model. A new crack damage creep model is established. The one-dimensional and three-dimensional damage creep equations are derived, and the model parameter identification method is given.

(4)

The applicability of the fracture damage creep model is verified by using the creep test data of fractured coal samples. The results show that the creep model can not only accurately reflect the creep characteristics of the decay and constant velocity creep curves of fractured coal samples but also describe the accelerated creep characteristics and effectively simulate the whole creep process of fractured coal samples.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (grant number 51574156) and the Key Development Program for Research of Shandong Province (grant number 2018GNC110023). This research was also partially funded by the Shandong Province Higher Educational Science and Technology Program (grant number J18KA195).

References

[1] P. Y. Shan, Z. M. Tao, X. B. Ye, "The analysis of rock burst in coalmine," Journal of Coal Science and Engineering, vol. 2 no. 1, pp. 3535-3544, 1996.

[2] R. S. Faradonbeh, A. Taheri, L. R. E. Sousa, M. Karakus, "Rock burst assessment in deep geotechnical conditions using true-triaxial tests and data-driven approaches," International Journal of Rock Mechanics and Mining Sciences, vol. 128, 2020.

[3] P. Wang, H. Jia, P. Zheng, "Sensitivity analysis of bursting liability for different coal-rock combinations based on their inhomogeneous characteristics," Geomatics, Natural Hazards and Risk, vol. 11 no. 1, pp. 149-159, DOI: 10.1080/19475705.2020.1714754, 2020.

[4] M. He, F. Ren, D. Liu, "Rockburst mechanism research and its control," International Journal of Mining Science and Technology, vol. 28 no. 5, pp. 829-837, DOI: 10.1016/j.ijmst.2018.09.002, 2018.

[5] J. Zhou, X. Li, X. Shi, "Long-term prediction model of rockburst in underground openings using heuristic algorithms and support vector machines," Safety Science, vol. 50 no. 4, pp. 629-644, DOI: 10.1016/j.ssci.2011.08.065, 2012.

[6] C. A. Tang, J. Wang, J. Zhang, "Preliminary engineering application of microseismic monitoring technique to rockburst prediction in tunneling of Jinping II project," Journal of Rock Mechanics and Geotechnical Engineering, vol. 2 no. 3, pp. 193-208, DOI: 10.3724/sp.j.1235.2010.00193, 2010.

[7] J. Zhou, X. Li, H. S. Mitri, "Evaluation method of rockburst: state-of-the-art literature review," Tunnelling and Underground Space Technology, vol. 81, pp. 632-659, DOI: 10.1016/j.tust.2018.08.029, 2018.

[8] R. S. Faradonbeh, A. Taheri, "Long-term prediction of rock burst hazard in deep underground openings using three robust data mining techniques," Engineering with Computers, vol. 35 no. 2, pp. 659-675, DOI: 10.1007/s00366-018-0624-4, 2019.

[9] X. Z. Shi, J. Zhou, L. Dong, H. Y. Hu, H. Y. Wang, S. R. Chen, "Application of unascertained measurement model to prediction of classification of rock burst intensity," Chinese Journal of Rock Mechanics and Engineering, vol. 29 no. 1, pp. 2720-2726, 2009.

[10] ASTM D5312, "Standard test method for evaluation of durability of rock for erosion control under freezing and thawing conditions," ASTM Standards, vol. 12 no. 6, 2013.

[11] K. Shivakumar, M. V. M. S. Rao, C. Srinivasan, K. Kusunose, "Multifractal analysis of the spatial distribution of area rockbursts at Kolar gold mines," International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 33 no. 2, pp. 167-172, DOI: 10.1016/0148-9062(95)00066-6, 1996.

[12] J.-A. Wang, H. D. Park, "Comprehensive prediction of rockburst based on analysis of strain energy in rocks," Tunnelling and Underground Space Technology, vol. 16 no. 1, pp. 49-57, DOI: 10.1016/s0886-7798(01)00030-x, 2001.

[13] M. Aubertin, D. E. Gill, R. Simon, "On the use of the brittleness index modified (BIM) to estimate the post-peak behavior or rocks," Aqua Fennica, vol. 23, pp. 24-25, 1994.

[14] C. Q. Zhang, J. J. Lu, J. Chen, H. Zhou, F. J. Yang, "Discussion on rock burst proneness indexes and their relation," Rock and Soil Mechanics, vol. 38 no. 5, pp. 1397-1404, 2017.

[15] J. P. Liu, X. T. Feng, "Case histories of rock burst at coal mines," pp. 93-113, .

[16] M. P. J. Schopfer, G. Zulauf, "Strain-dependent rheology and the memory of plasticine," Tectonophysics, vol. 354 no. 2, pp. 85-99, 2002.

[17] K. Miura, Y. Okui, "Micromechanics-based prediction of creep failure of hard rock for long-term safety of high-level radioactive waste disposal system," Mechanics of Materials, vol. 35 no. 6, pp. 587-601, DOI: 10.1016/s0167-6636(02)00286-7, 2003.

[18] D. Griggs, "Creep of rocks," The Journal of Geology, vol. 47 no. 3, pp. 225-251, 1939.

[19] K. Terzaghi, "Stability of steep slopes on hard unweathered rock," Géotechnique, vol. 12 no. 4, pp. 251-270, DOI: 10.1680/geot.1962.12.4.251, 1962.

[20] T.-K. Tan, W.-F. Kang, "Locked in stresses, creep and dilatancy of rocks, and constitutive equations," Rock Mechanics Felsmechanik Mecanique des Roches, vol. 13 no. 1,DOI: 10.1007/bf01257895, 1980.

[21] Y. Z. Jiang, M. M. Zhang, L. Q. Li, "Study on nonlinear viscoelasto-plastic creep model of rock and its parameters identification," Chinese Journal of Rock Mechanics and Engineering, vol. 27 no. 4, pp. 833-839, 2008.

[22] H. T. Xiao, T. C. Qiang, W. Y. Zhou, "Rheological damage study of three gorges ship lock high slope and analysis on in-situ measured deformation," Chinese Journal of Rock Mechanics and Engineering, vol. 18 no. 5, pp. 497-502, 1999.

[23] Y. L. Zhao, W. Wan, W. J. Wang, "Compressive-shear rheological fracture of rock-like cracks and subcritical crack propagation test and fracture mechanism," Chinese Journal of Geotechnical Engineering, vol. 34 no. 6, pp. 1050-1059, 2012.

[24] K. Matsuki, E. Q. Wang, K. Sakaguchi, K. Okumura, "Time-dependent closure of a fracture with rough surfaces under constant normal stress," International Journal of Rock Mechanics and Mining Sciences, vol. 38 no. 5, pp. 607-619, DOI: 10.1016/s1365-1609(01)00022-3, 2001.

[25] C. Yang, D. Huang, R. Q. Huang, "Experimental research on the triaxial unloading creep characteristics of stone with two pre-existing cracks and its damage creep model," Journal of China Coal Society, vol. 41 no. 9, pp. 2203-2211, 2016.

[26] S. Murakami, N. Ohno, "A continuum theory of creep and creep damage," Proceedings of the 3rd IUTAM Symposium on Creep in Structure, pp. 422-443, .

[27] T. Kawamoto, Y. Ichikawa, T. Kyoya, "Deformation and fracturing behaviour of discontinuous rock mass and damage mechanics theory," International Journal for Numerical and Analytical Methods in Geomechanics, vol. 12 no. 1,DOI: 10.1002/nag.1610120102, 1988.

[28] W. S. Zhu, Y. P. Qi, Y. H. Guo, "Three dimensional damage rheological analysis of surrounding rock deformation and fracture of Jinping I hydropower station," Journal of Rock Mechanics and Engineering, vol. 31 no. 5, pp. 865-872, 2012.

[29] J. Lemaitre, "How to use damage mechanics," Nuclear Engineering and Design, vol. 80 no. 3, pp. 233-245, 1984.

[30] M. Kachanov, "Effective elastic properties of cracked solids: critical review of some basic concepts," Applied Mechanics Reviews, vol. 45 no. 8, pp. 304-335, DOI: 10.1115/1.3119761, 1992.