**1. Introduction**

During the excavation of coal mine shafts and roadways, the stress of the surrounding rock is redistributed, and the surrounding rock is repeatedly disturbed by construction. Therefore, there are a large number of micro defects such as cracks and cavities in the rock mass, resulting in the nonlinear mechanical characteristics of the rock mass under various stress states [1,2]. The mechanical properties of rocks are particularly important for the stability analysis of surrounding rocks [3,4], so it is necessary to analyze the mechanical properties of rock under different stress states. Establishing the constitutive model of rock, using statistical damage theory, has become one of the important methods to study the nonlinear mechanical properties of rock [5–9].

Experts at home and abroad have produced much research on the statistical damage constitutive model of rock materials. The key to statistical damage theory is the reasonable measurement of rock micro element strength [10]. Tang Chunan [11] proposed to measure the rock micro-element strength through axial strain, which has achieved satisfactory results, but has ignored the influence of the stress state of the rock micro-element on its strength, therefore, this method has some shortcomings. Cao Wengui et al. [12] first proposed a new rock micro element strength measurement method based on Mohr-Coulomb

**Citation:** Zhou, R.; Guo, L.; Hong, R. Study on Energy Evolution and Damage Constitutive Model of Siltstone. *Crystals* **2021**, *11*, 1271. https://doi.org/10.3390/ cryst11111271

Academic Editors: Peter Taylor, Per-Lennart Larsson, Yifeng Ling, Chuanqing Fu and Peng Zhang

Received: 14 September 2021 Accepted: 14 October 2021 Published: 20 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

criterion. Mohr-Coulomb criterion can accurately reflect the bearing capacity of rock, but it ignores the influence of intermediate principal stress, so the constitutive model established also has some limitations. Xu Weiya et al. [13] assumed that rock micro-element failure conforms to the Drucker–Plager criterion and established a statistical damage constitutive model of rock based on Weibull distribution. However, the Drucker–Plager criterion is relatively conservative, which limits the rationality of the rock micro-element strength that is determined based on this criterion. In view of the shortcomings of the above yield criterion, the Hoek–Brown criterion proposed by Hoek et al. [14,15] can better reflect the nonlinear failure characteristics of rock and the confining pressure effect of rock strength; Cao Ruilang et al. [16] established a three-dimensional statistical damage constitutive model under the condition that the micro-element failure of rock conforms to the Hoek–Brown criterion, which can reflect the post peak softening characteristics of rock by using the post-peak residual strength of rock to modify the damage variable. Most of the above studies use the mechanical parameters of rock as the basis for establishing the statistical damage constitutive equation, whereas only a few studies introduce the energy evolution characteristics of rock into the constitutive relationship [17]. Gao Wei et al. [18] established the damage constitutive model of granite under uniaxial compression by using the principle of minimum energy consumption. Sun Mengcheng et al. [19] established a new damage constitutive model by introducing the rock unified energy yield criterion as the energy consumption constraint condition, which is based on the principle of minimum energy consumption and continuous damage theory, but did not consider the influence of the post-peak softening stage on the damage constitutive model. It can be seen that, when studying the nonlinear constitutive relationship of rock under complex stress conditions, the minimum energy consumption principle is rarely applied to the establishment of rock damage evolution equation, and the influence of post-peak strain softening stage on the constitutive model is rarely considered in the existing damage constitutive model based on the energy principle.

In summary, much progress has been made in the research on rock damage mechanics, but there are still some problems, specifically the following: (1) there are weaknesses in the definition of the damage variable. Some scholars establish damage variables according to the statistical distribution theory, but the statistical distribution function cannot accurately describe the crack development and deformation in rock materials. Therefore, the rock damage constitutive model based on the statistical distribution theory is not in accordance with the actual situation of rock damage and deformation. (2) There are weaknesses in the establishment of the damage model. Hooke's law can only describe the elastic stress–strain relationship of rock in the pre-peak stage, but cannot describe the stress–strain relationship in the post-peak strain softening stage. Some studies do not consider the strain softening stage, and therefore there are weaknesses in the damage constitutive model.

In order to avoid the above problems, a constitutive model that can accurately describe the characteristics of rock damage deformation and residual strength is constructed. Based on the principle of minimum energy consumption, this paper introduces the improved three-shear energy yield criterion as the energy consumption constraint, considers the influence of the post-peak strain softening stage and deduces the damage evolution equation of rock. According to the principle of effective stress, the damage constitutive model of rock under a conventional triaxial compression is established. In addition, the conventional triaxial compression tests of siltstone under different confining pressures are carried out, and the model parameters are identified according to the test data. The rationality and feasibility of the model in this paper are verified by comparing the model results with the test results.

### **2. Conventional Triaxial Compression Test of Siltstone**

### *2.1. Test Process*

The test material used is Permian upper Shihezi Formation siltstone, which is taken from matoumen—at a 425 m depth of an east air shaft in Yuandian No. 2 mine, Huaibei

City, Anhui Province. According to the standards for test methods of engineering rock mass (GBT50266-2013), the sample size is a cylinder with a diameter of 50 mm, a length of 100 mm and a length diameter ratio of 2. 2S-200 using a vertical coring machine, DQ-4 rock cutter and SHM-200 double face grinder to drill, cut and grind the siltstone rock sample. The disparity between the sample is less than 0.05 mm and the diameter error of the sample is no more than 0.3 mm. City, Anhui Province. According to the standards for test methods of engineering rock mass (GBT50266-2013), the sample size is a cylinder with a diameter of 50 mm, a length of 100 mm and a length diameter ratio of 2. 2S-200 using a vertical coring machine, DQ-4 rock cutter and SHM-200 double face grinder to drill, cut and grind the siltstone rock sample. The disparity between the sample is less than 0.05 mm and the diameter error of the sample is no more than 0.3 mm.

The test material used is Permian upper Shihezi Formation siltstone, which is taken from matoumen—at a 425 m depth of an east air shaft in Yuandian No. 2 mine, Huaibei

*Crystals* **2021**, *11*, 1271 3 of 18

**2. Conventional Triaxial Compression Test of Siltstone**

*2.1. Test Process*

The conventional triaxial compression test of siltstone samples is performed using a ZTCR-2000 low temperature rock triaxial system (Figure 1). During the test, the siltstone samples are preloaded to 0.5 MPa, and then the confining pressure is applied to the predetermined value at a loading speed of 50 N/s according to the load control mode. The test confining pressure is set to 0, 5, 10, 15 and 20 MPa. After the confining pressure reaches the predetermined value, it is stabilized for 30 s, and finally the press applies an axial pressure at the loading speed of 0.06 mm/min in the way of displacement control until the siltstone sample is damaged. The conventional triaxial compression test of siltstone samples is performed using a ZTCR-2000 low temperature rock triaxial system (Figure 1). During the test, the siltstone samples are preloaded to 0.5 MPa, and then the confining pressure is applied to the predetermined value at a loading speed of 50 N/s according to the load control mode. The test confining pressure is set to 0, 5, 10, 15 and 20 MPa. After the confining pressure reaches the predetermined value, it is stabilized for 30 s, and finally the press applies an axial pressure at the loading speed of 0.06 mm/min in the way of displacement control until the siltstone sample is damaged.

**Figure 1.** Rock mechanics test system. **Figure 1.** Rock mechanics test system.

### *2.2. Analysis of Test Results 2.2. Analysis of Test Results*

The siltstone samples were divided into three groups with five samples in each group. The tests were repeated three times under the same loading type, and the representative test data were selected for analysis. Figure 2a shows the stress–strain curve of siltstone samples under different confining pressures. It can be seen from Figure 2a that when the confining pressure is at 0 MPa in the pre-peak stage, due to the lack of confining pressure constraints when the loading exceeds the elastic limit of the sample the siltstone sample quickly reaches the failure state, and the pre-peak curve of the sample has no obvious yield stage. After reaching the peak strength, the stress decreases rapidly, and the post-peak curve is steep, displaying an obvious strain softening. When the confining pressure is at 20 MPa, because the confining pressure effectively limits the propagation speed of microcracks in the sample and slows down the damage degree of the sample, the yield stage is obvious in the pre-peak curve. The post-peak curve and stress drop trend tend to be gentle, and the strain softening phenomenon decreases in the post-peak stage. With the increase of the confining pressure, the peak stress and peak strain increase gradually, indicating that the bearing capacity of the siltstone sample increases gradually, and it becomes more difficult for the sample to enter the failure state. The failure characteristics of siltstone samples under different confining pressures are shown in The siltstone samples were divided into three groups with five samples in each group. The tests were repeated three times under the same loading type, and the representative test data were selected for analysis. Figure 2a shows the stress–strain curve of siltstone samples under different confining pressures. It can be seen from Figure 2a that when the confining pressure is at 0 MPa in the pre-peak stage, due to the lack of confining pressure constraints when the loading exceeds the elastic limit of the sample the siltstone sample quickly reaches the failure state, and the pre-peak curve of the sample has no obvious yield stage. After reaching the peak strength, the stress decreases rapidly, and the postpeak curve is steep, displaying an obvious strain softening. When the confining pressure is at 20 MPa, because the confining pressure effectively limits the propagation speed of microcracks in the sample and slows down the damage degree of the sample, the yield stage is obvious in the pre-peak curve. The post-peak curve and stress drop trend tend to be gentle, and the strain softening phenomenon decreases in the post-peak stage. With the increase of the confining pressure, the peak stress and peak strain increase gradually, indicating that the bearing capacity of the siltstone sample increases gradually, and it becomes more difficult for the sample to enter the failure state. The failure characteristics of siltstone samples under different confining pressures are shown in Figure 2b.

Figure 2b. The basic mechanical parameters of siltstone are obtained according to the total stress– strain curve. See Table 1.

3 

Notes:

3

*W*<sup>E</sup>

where

 1 3 ,

**Figure 2.** Conventional triaxial test results of siltstone samples: (**a**) Full stress–strain curves of siltstone under different confining pressures; (**b**) Failure characteristics of siltstone samples. **Figure 2.** Conventional triaxial test results of siltstone samples: (**a**) Full stress–strain curves of siltstone under different confining pressures; (**b**) Failure characteristics of siltstone samples.


The basic mechanical parameters of siltstone are obtained according to the total **Table 1.** Mechanical parameters of siltstone under conventional triaxial test.

5 75.57 4.32 34.86 5.93 17.6 0.26 10 95.45 4.72 61.29 6.37 20.0 0.26 Notes: *σ*<sup>3</sup> is the confining pressure, *σsc* is the peak stress, *εsc* is the peak strain, *σ r* 1 is the residual stress, *ε<sup>r</sup>* is the strain corresponding to the residual stress, *E* is the elastic modulus and *µ* is the Poisson's ratio.

### 15 115.09 5.26 64.54 8.17 21.0 0.25 20 134.38 5.81 86.40 9.80 21.5 0.25 **3. Energy Analysis of Triaxial Compression Process of Siltstone**

### is the confining pressure, is the peak stress, is the peak strain, *3.1. Theoretical Analysis of Energy Evolution*

 *sc sc* 1 *r* strain corresponding to the residual stress, *E* is the elastic modulus and is the Poisson's ratio. **3. Energy Analysis of Triaxial Compression Process of Siltstone** Assuming that heat exchange does not occur between the rock system and the external environment, according to the first law of thermodynamics, during the process of rock deformation and failure, the input energy *W*<sup>F</sup> is equal to the elastic strain energy *W*<sup>E</sup> plus the dissipative energy *W*<sup>D</sup> [20,21].

*3.1. Theoretical Analysis of Energy Evolution* The above-mentioned energy relationship is shown in Formula (1):

$$\mathcal{W}\_{\rm F} = \mathcal{W}\_{\rm E} + \mathcal{W}\_{\rm D} \tag{1}$$

is the residual stress,

is the

of rock deformation and failure, the input energy *W*<sup>F</sup> is equal to the elastic strain energy plus the dissipative energy *W*<sup>D</sup> [20,21]. The work performed by the axial force and confining pressure during the test is [22]:

Assuming that heat exchange does not occur between the rock system and the

$$W\_{\rm F} = \frac{\pi}{4} D^2 H \left( \int\_0^{\varepsilon\_1} \sigma\_1 \mathbf{d} \varepsilon\_1 + 2 \int\_0^{\varepsilon\_3} \sigma\_3 \mathbf{d} \varepsilon\_3 \right) = V \mathcal{U}\_{\rm F} \tag{2}$$

*r*

The work performed by the axial force and confining pressure during the test is [22]: 1 3 <sup>2</sup> F 1 1 3 3 F 0 0 d 2 d 4 *π <sup>ε</sup> <sup>ε</sup> W D H <sup>σ</sup> <sup>ε</sup> <sup>σ</sup> <sup>ε</sup> VU* (2) where *σ*1, *σ*<sup>3</sup> are the axial pressure and confining pressure, respectively, *ε*1, *ε*<sup>3</sup> are the axial and radial strain, respectively, *D*, *H* are the diameter and height of the siltstone specimen, respectively, *V* is the volume of the siltstone sample and *U*<sup>F</sup> is the input energy density. In the same way, the elastic strain energy and dissipative energy are as follows:

$$\begin{cases} \begin{array}{l} \text{W}\_{\text{E}} = \frac{\pi}{4} D^2 H \text{U}\_{\text{E}} = V \text{U}\_{\text{E}}\\ \text{W}\_{\text{D}} = \frac{\pi}{4} D^2 H \text{U}\_{\text{D}} = V \text{U}\_{\text{D}} \end{array} \tag{3}$$

siltstone specimen, respectively, is the volume of the siltstone sample and *U*<sup>F</sup> is the input energy density. In the same way, the elastic strain energy and dissipative energy are as follows: where *U*E, *U*<sup>D</sup> are the elastic strain energy density and dissipative energy density, respectively. According to the elastic theory [22], the elastic strain energy density is:

$$\mathcal{U}\_{\rm E} = \frac{1}{2} (\sigma\_1 \mathfrak{e}\_1^{\rm e} + \sigma\_2 \mathfrak{e}\_2^{\rm e} + \sigma\_3 \mathfrak{e}\_3^{\rm e}) \tag{4}$$

*V*

The three-dimensional constitutive relationship of siltstone is:

$$
\varepsilon\_{ij}^{\varepsilon} = \frac{1+\mu}{\mathbb{E}\_{ij}} \sigma\_{ij} - \frac{\mu}{\mathbb{E}\_{ij}} \sigma\_{kk} \delta\_{ij} \tag{5}
$$

where *ε* e *ij*(*i*, *j* = 1, 2, 3) is the elastic strain in the direction of the main stress, *σij*(*i*, *j* = 1, 2, 3) is the main stress, *σkk* = *σ*<sup>1</sup> + *σ*<sup>2</sup> + *σ*3, *δij* is the Kronecker tensor, *Eij*(*i*, *j* = 1, 2, 3) is replaced by the initial elastic modulus and *E* [22], *µ* is the Poisson's ratio.

Equations (4) and (5) can be substituted into Equation (3) to obtain the elastic strain energy *W*E:

$$\mathcal{W}\_{\rm E} = \frac{1}{2E} \left( \sigma\_1^2 + \sigma\_2^2 + \sigma\_3^2 - 2\mu\sigma\_1\sigma\_2 - 2\mu\sigma\_1\sigma\_3 - 2\mu\sigma\_2\sigma\_3 \right) V \tag{6}$$

In the conventional triaxial compression test where *σ*<sup>2</sup> = *σ*3, Formula (6) can be simplified to:

$$\mathcal{W}\_{\rm E} = \frac{1}{2E} \left[ \sigma\_1^2 + 2(1 - \mu)\sigma\_3^2 - 4\mu\sigma\_1\sigma\_3 \right] V \tag{7}$$

Substituting Formula (7) into Formula (1) and combining this with Formula (2) to obtain the dissipative energy results in the following:

$$\mathcal{W}\_{\rm D} = \left\{ \int\_0^{\varepsilon\_1} \sigma\_1 \mathbf{d} \varepsilon\_1 + 2 \int\_0^{\varepsilon\_3} \sigma\_3 \mathbf{d} \varepsilon\_3 - \frac{1}{2E} \left[ \sigma\_1^2 + 2(1 - \mu)\sigma\_3^2 - 4\mu \sigma\_1 \sigma\_3 \right] \right\} V \tag{8}$$

### *3.2. Principle of Minimum Energy Consumption*

For dissipative materials such as geotechnical materials, the stress comes from internal variables such as strain and temperature. Therefore, the internal variables reflecting the internal changes of materials need to be considered to truly establish the constitutive relationship of dissipative materials [23].

According to the internal variable theory, for any infinitesimal dissipative micro element, the volume dissipation rate is:

$$\mathfrak{q}\rho\_0 = \sigma\_{\dot{\imath}} : \boldsymbol{\varepsilon}^N + \boldsymbol{Y} : \boldsymbol{D} + \sum\_{i=1}^k \mathsf{R}\_i \overset{\bullet}{\gamma}\_i - q\frac{\mathcal{g}}{T} \tag{9}$$

where, *ρ*<sup>0</sup> is the unit weight of rock, *ϕ* is the energy consumption rate of rock, *D* is the damage variable of material, *Y* is the dual variable corresponding to the damage variable, *Ri* is the corresponding dual variable of *γ<sup>i</sup>* .

Without considering heat dissipation, the above formula is simplified as:

$$
\rho \rho\_0 = \sigma\_{\hat{\imath}} : \stackrel{\bullet}{\varepsilon^N} + \mathcal{Y} : D + \sum\_{i=1}^k R\_i \stackrel{\bullet}{\gamma\_i} \tag{10}
$$

It is generally assumed that dissipative materials meet the following energy dissipation constraints in the process of energy dissipation:

$$\begin{cases} \, \, \_F\_1(\sigma, \, \_Y\mathbb{R}\_1, \dots \mathbb{R}\_k) = 0\\ \, \_F\_m(\sigma, \, \_Y\mathbb{R}\_1, \dots \mathbb{R}\_k) = 0 \end{cases} \tag{11}$$

According to the principle of minimum energy consumption, Formula (10) takes the stationary value under the condition of Formula (11) and introduces the Lagrange multiplier *λ* to obtain:

$$\frac{\partial(\varrho + \lambda\_i F\_i)}{\partial \xi\_i^\*} = 0 \ (i = 1, \ldots, m) \tag{12}$$

where, *ξ<sup>i</sup>* takes *σ*,*Y*, *R*1, . . . *R<sup>k</sup>* respectively. The Lagrange multiplier *λi*(*i* = 1, . . . , *m*) is introduced into Equation (10) to obtain:

$$\begin{cases} \begin{aligned} \stackrel{\bullet}{\varepsilon}\_{N} + \sigma &: \stackrel{\bullet}{\frac{\partial \bullet}{\partial \sigma}}\_{\delta \sigma} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial \sigma} = 0 \\ \stackrel{\bullet}{D} + Y &: \stackrel{\bullet}{\frac{\partial \bullet}{\partial Y}}\_{\delta Y} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial Y} = 0 \\ \stackrel{\bullet}{\gamma\_{1} + R\_{1}} &: \stackrel{\bullet}{\frac{\partial \bullet}{\partial R\_{1}}}\_{\delta R\_{1}} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial R\_{1}} = 0 \\ & \cdot \\ & \cdot \\ & \cdot \\ \stackrel{\bullet}{\gamma\_{k} + R\_{k}} &: \stackrel{\bullet}{\frac{\partial \bullet}{\partial R\_{k}}}\_{k} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial R\_{k}} = 0 \end{aligned} \tag{13}$$

Assuming *ρ*0*ϕ* is a potential function, we get:

$$\begin{cases} \stackrel{\bullet}{e^{N}} + \stackrel{m}{\sum}\_{i=1} \lambda\_{i} \frac{\partial F\_{i}}{\partial \sigma} = 0\\ \stackrel{\bullet}{D} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial Y} = 0\\ \stackrel{\bullet}{\gamma\_{1}} + \sum\_{i=1}^{m} \lambda\_{i} \frac{\partial F\_{i}}{\partial R\_{1}} = 0\\ \vdots\\ \stackrel{\bullet}{\gamma\_{k}} + \stackrel{\bullet}{\sum}\_{i=1} \lambda\_{i} \frac{\partial F\_{i}}{\partial R\_{k}} = 0 \end{cases} (14)$$

The above formula is the internal variable evolution equation of energy dissipation of dissipative materials derived based on the principle of minimum energy dissipation.

### *3.3. Relationship between Energy Evolution and Axial Strain*

The triaxial compression test results are substituted into Equations (2), (7) and (8) to obtain the evolution curves of the input energy, elastic energy and dissipative energy of siltstone with an axial strain under different confining pressures, as shown in Figure 3.

It can be seen from Figure 3 that under different confining pressures, the input energy and dissipative energy increase with the increase of axial strain, and the elastic strain energy increases at first and then decreases. At the initial stage of sample loading, the initial pores in the rock are gradually closed as a result of the action of the external load, most of the work performed by the external load is transformed into elastic energy and stored in the sample and the elastic strain energy gradually increases with axial strain, the dissipative energy at this stage is very small, and the evolution curves of input energy, elastic energy and dissipative energy concave upward. Alongside the increase of the external load, the siltstone sample enters the linear elastic deformation stage, and the external work is essentially transformed into elastic strain energy, and the slope of the three energy evolution curves reaches the maximum. This stage is the main stage of energy storage in the overall process of rock failure. When the external load reaches the yield limit of rock, new cracks will appear in the sample, and part of the energy is required to be dissipated for its initiation and diffusion. Therefore, the slope of the elastic strain energy curve decreases gradually at this stage. When the external load reaches the peak strength, the elastic strain energy reaches the maximum value, and as a result the siltstone sample is damaged, and the elastic strain energy stored in the pre peak stage is released rapidly. Therefore, the elastic strain energy after the peak decreases gradually with the axial strain, and most of the input energy is dissipated in the process of the mutual penetration of cracks to form a macro-fracture surface. When the sample reaches its peak strength, the

*Crystals* **2021**, *11*, 1271 7 of 18

*3.3. Relationship between Energy Evolution and Axial Strain*

elastic strain energy decreases gradually until the sample failure reaches the minimum value, and the dissipative energy increases gradually until the sample failure reaches the maximum value. The triaxial compression test results are substituted into Equations (2), (7) and (8) to obtain the evolution curves of the input energy, elastic energy and dissipative energy of siltstone with an axial strain under different confining pressures, as shown in Figure 3.

**Figure 3.** Energy evolution curve of siltstone under different confining pressures: (**a**) 0 MPa; (**b**) 5 MPa; (**c**) 10 MPa; (**d**) 15 MPa; (**e**) 20 MPa. **Figure 3.** Energy evolution curve of siltstone under different confining pressures: (**a**) 0 MPa; (**b**) 5 MPa; (**c**) 10 MPa; (**d**) 15 MPa; (**e**) 20 MPa.

### It can be seen from Figure 3 that under different confining pressures, the input energy *3.4. Relationship between Energy Evolution and Confining Pressure*

and dissipative energy increase with the increase of axial strain, and the elastic strain energy increases at first and then decreases. At the initial stage of sample loading, the initial pores in the rock are gradually closed as a result of the action of the external load, most of the work performed by the external load is transformed into elastic energy and stored in the sample and the elastic strain energy gradually increases with axial strain, the dissipative energy at this stage is very small, and the evolution curves of input energy, elastic energy and dissipative energy concave upward. Alongside the increase of the external load, the siltstone sample enters the linear elastic deformation stage, and the external work is essentially transformed into elastic strain energy, and the slope of the three energy evolution curves reaches the maximum. This stage is the main stage of energy storage in the overall process of rock failure. When the external load reaches the yield limit of rock, new cracks will appear in the sample, and part of the energy is required to be dissipated for its initiation and diffusion. Therefore, the slope of the elastic strain energy curve decreases gradually at this stage. When the external load reaches the peak strength, the elastic strain energy reaches the maximum value, and as a result the siltstone sample is damaged, and the elastic strain energy stored in the pre peak stage is released rapidly. Therefore, the elastic strain energy after the peak decreases gradually with the axial strain, and most of the input energy is dissipated in the process of the mutual penetration of cracks to form a macro-fracture surface. When the sample reaches its peak Based on the triaxial compression test data of siltstone under different confining pressures, the characteristic energy results of siltstone are calculated by substituting Equations (2), (7) and (8), as shown in Table 2. It can be seen from Table 2 that when the confining pressure is at 0 MPa, the input energy *W*F(pre), elastic energy *W*E(pre) and dissipative energy *W*D(pre), corresponding to the peak stress of the siltstone sample, are 13.97 J, 11.97 J and 2.00 J, respectively. Most of the input energy before the peak of the sample is transformed into storable elastic energy, and only a small amount of the energy is dissipated in the process of damage deformation and of crack propagation of the sample, indicating that the energy storage capacity of the sample is strong, Therefore, the elastic property is the source power of the specimen failure. When the confining pressure is at 20 MPa, the input energy *W*F(pre), elastic energy *W*E(pre) and dissipative energy *W*D(pre), corresponding to the peak stress of the siltstone sample, are 91.46 J, 52.01 J and 39.45 J respectively. The elastic energy that was stored before the peak accounts for 56.8% of the input energy. The energy dissipated via the damage and deformation of the sample accounts for 43.2% of the input energy, indicating that the energy storage capacity of the sample is weakened at this time. The energy required for specimen failure is involved in the work provided by the external testing machine, and the self-sustaining fracture ability of the specimen is weak.


**Table 2.** Calculation results of characteristic energy of siltstone sample at peak stress.

The evolution law of input energy, elastic strain energy and dissipative energy of the siltstone sample with confining pressure is shown in Figure 4. It can be seen from Figure 4, that with the increase of confining pressure, the three energies increase at different rates, the input energy and dissipative energy are in an exponential function proportional relationship with the confining pressure, the elastic energy is in a linear function proportional relationship with the confining pressure, and the storage rate of elastic strain energy gradually increases, indicating that the confining pressure has an obvious restrictive effect on crack propagation, thus limiting the siltstone sample to the failure state. *Crystals* **2021**, *11*, 1271 9 of 18

**Figure 4.** Characteristic energy of siltstone under different confining pressures. **Figure 4.** Characteristic energy of siltstone under different confining pressures.

### **4. Criterion for Determining Rock Damage Threshold 4. Criterion for Determining Rock Damage Threshold**

The yield criterion of rock materials generally included the following three categories: the yield criterion established from the angle of stress, strain and energy [24]. There are abundant yield criterion that are established from the angle of stress, mainly the Mises criterion, the Drucker-Plager criterion, the Mohr–Coulomb criterion, the Hoek– Brown criterion and the double shear strength criterion. Among them, the Mohr– Coulomb criterion is the most widely used, is applicable to both plastic and brittle rocks, and can reflect the characteristic that the tensile strength of rock is less than the compressive strength. However, the criterion does not consider the influence of the intermediate principal stress on rock yield, and the envelope line of the criterion on the meridional plane is a straight line, indicating that the internal friction angle does not change with hydrostatic pressure; this is inconsistent with the test results of rock mechanics [25,26]. The Drucker–Plager criterion considers the influence of intermediate principal stress on rock yield, but it cannot distinguish the difference between the tensile meridian and compressive meridian of rock, which is inconsistent with the triaxial test results of rock [27]. Double shear strength criterion is only applicable to materials for which the shear , tensile and compressive strength meet a certain relationship [28]. The accuracy of establishing yield criterion from the point of view of stress depends on the description of yield stress, and the expression does not include material parameters, so its application range is limited. The establishment of a yield criterion from the perspective of energy can accurately express the yield state of rock with the assistance of the description The yield criterion of rock materials generally included the following three categories: the yield criterion established from the angle of stress, strain and energy [24]. There are abundant yield criterion that are established from the angle of stress, mainly the Mises criterion, the Drucker-Plager criterion, the Mohr–Coulomb criterion, the Hoek–Brown criterion and the double shear strength criterion. Among them, the Mohr–Coulomb criterion is the most widely used, is applicable to both plastic and brittle rocks, and can reflect the characteristic that the tensile strength of rock is less than the compressive strength. However, the criterion does not consider the influence of the intermediate principal stress on rock yield, and the envelope line of the criterion on the meridional plane is a straight line, indicating that the internal friction angle does not change with hydrostatic pressure; this is inconsistent with the test results of rock mechanics [25,26]. The Drucker–Plager criterion considers the influence of intermediate principal stress on rock yield, but it cannot distinguish the difference between the tensile meridian and compressive meridian of rock, which is inconsistent with the triaxial test results of rock [27]. Double shear strength criterion is only applicable to materials for which the shear, tensile and compressive strength meet a certain relationship [28]. The accuracy of establishing yield criterion from the point of view of stress depends on the description of yield stress, and the expression does not include material parameters, so its application range is limited. The establishment of a yield criterion from the perspective of energy can accurately express the yield state of rock with the assistance of the description of energy evolution. Therefore, the establishment

of energy evolution. Therefore, the establishment of an energy yield criterion across a wide range of applications is a fundamental method to analyze the yield deformation of

The yield criterion is established from the perspective of energy, mainly by determining the functional relationship between the strain energy during material yield [29]. The mises energy criterion defines the maximum shape change of specific energy as a constant, but this criterion cannot reflect the characteristics of the different tensile and compressive strengths of rock materials. Gao Hong et al. [30] proposed the three-shear energy yield criterion of rock materials based on the maximum strain specific energy, and introduced the expression of shear strain energy by considering the internal friction angle. The criterion is simple and is a generalization of many common yield criteria. However, the three-shear energy yield criterion assumes that the sum of the shear strain energy of a composite sliding surface at the time of the material yield is a constant, and does not consider the influence of hydrostatic pressure on the material yield, which leads to a

In conclusion, it is highly necessary to establish a yield criterion that can reflect the internal friction characteristics and hydrostatic pressure effects of rock materials. Based on the test results, this paper explores the functional relationship between the sum of the

considerable difference between and the test results of rock materials.

rock.

of an energy yield criterion across a wide range of applications is a fundamental method to analyze the yield deformation of rock.

The yield criterion is established from the perspective of energy, mainly by determining the functional relationship between the strain energy during material yield [29]. The mises energy criterion defines the maximum shape change of specific energy as a constant, but this criterion cannot reflect the characteristics of the different tensile and compressive strengths of rock materials. Gao Hong et al. [30] proposed the three-shear energy yield criterion of rock materials based on the maximum strain specific energy, and introduced the expression of shear strain energy by considering the internal friction angle. The criterion is simple and is a generalization of many common yield criteria. However, the three-shear energy yield criterion assumes that the sum of the shear strain energy of a composite sliding surface at the time of the material yield is a constant, and does not consider the influence of hydrostatic pressure on the material yield, which leads to a considerable difference between and the test results of rock materials.

In conclusion, it is highly necessary to establish a yield criterion that can reflect the internal friction characteristics and hydrostatic pressure effects of rock materials. Based on the test results, this paper explores the functional relationship between the sum of the shear strain energy of the three composite sliding surfaces and the hydrostatic pressure during rock yield, and establishes an improved three-shear energy yield criterion.

### *4.1. Three-Shear Energy Yield Criterion*

The three-shear energy yield criterion assumes that the sum of the shear strain energy *w<sup>s</sup>* of the three composite sliding surfaces is constant when the rock material yields. For the rock material with internal friction characteristics, *w<sup>s</sup>* plays a vital role in the rock yield. The following is a simple derivation of *w<sup>s</sup>* . It agrees that the compressive stress is positive in the derivation process.

The limit value of the shear strain energy of the three yield sliding surfaces of rock mass is:

$$w\_{12} = \frac{1}{2G} \left[ \frac{\sigma\_1 - \sigma\_2}{2 \cos \varphi\_{12}} - \frac{\sigma\_1 + \sigma\_2}{2} \tan \varphi\_{12} \right]^2 \tag{15}$$

$$w\_{23} = \frac{1}{2G} \left[ \frac{\sigma\_2 - \sigma\_3}{2 \cos \varphi\_{23}} - \frac{\sigma\_2 + \sigma\_3}{2} \tan \varphi\_{23} \right]^2 \tag{16}$$

$$w\_{13} = \frac{1}{2G} \left[ \frac{\sigma\_1 - \sigma\_3}{2 \cos \varphi\_{13}} - \frac{\sigma\_1 + \sigma\_3}{2} \tan \varphi\_{13} \right]^2 \tag{17}$$

where, *G* is the shear modulus.

The three-shear energy yield criterion considers that the rock material begins to yield when the sum of the shear strain energy of the three maximum friction angle action surfaces reaches a certain value, as follows:

$$w\_s = w\_{12} + w\_{23} + w\_{13} = k\_0 \tag{18}$$

In the case of uniaxial compression, *σ*<sup>2</sup> = *σ*<sup>3</sup> = 0, the value of *k*<sup>0</sup> is obtained:

$$k\_0 = \frac{c\_0^2}{G} \tag{19}$$

where, *c*<sup>0</sup> is the cohesion of rock at yield.

### *4.2. Improved Three-Shear Energy Yield Criterion*

The research results of Hao Tiesheng et al. [31] reveals that the sum of shear strain energy of three composite sliding surfaces during rock yield is not constant, and rock yield has an obvious hydrostatic pressure effect. Therefore, the functional relationship between the sum of the shear strain energy during rock yield and the hydrostatic pressure is established as follows:

$$F(w\_{\rm s}, p) = 0 \tag{20}$$

where, *p* is the hydrostatic pressure and *W<sup>s</sup>* is the shear strain energy of rock.

In order to establish the functional relationship between the sum of shear strain energy and hydrostatic pressure, by considering the sum of shear strain energy *W<sup>s</sup>* as the ordinate and hydrostatic pressure *p* as the abscissa, the variation law of the sum of shear strain energy of the three composite sliding surfaces with hydrostatic pressure during rock yield can be obtained. The two approximate accordance with the first-order functional relationship, so the energy yield criterion of siltstone is as follows [31]:

$$F(\sigma\_1, \sigma\_3) = \frac{1}{4G} \left[ \frac{\sigma\_1 - \sigma\_3}{\cos \varphi\_s} - (\sigma\_1 + \sigma\_3) \tan \varphi\_s \right]^2 + a(\sigma\_1 + 2\sigma\_3) + b \tag{21}$$

where, *a* and *b* are material parameters, representing the numerical relationship between shear strain energy and hydrostatic pressure and the rock yield, *ϕ<sup>s</sup>* is the internal friction angle when rock yield and *G* is the shear modulus of rock.

### **5. Establishment of Damage Constitutive Model of Siltstone**

## *5.1. Damage Constitutive Relationship*

Assuming that the rock is an isotropic material, according to J. Lemaitre's [32] strain equivalence hypothesis, the constitutive relationship of the rock is established as follows:

$$
\sigma\_i = \sigma\_i'(1 - D) + \sigma\_i'' D \ (i = 1, 2, 3) \tag{22}
$$

where, *σ<sup>i</sup>* is the nominal stress of the rock micro-element, *σ* 0 *i* is the effective stress of the rock micro-element, *σ* 00 *i* is the stress on the damaged part of the rock micro-element, and *D* is the damage variable.

Since the damaged part of the rock is closely associated with the undamaged part of the rock, according to the deformation coordination principle, *ε<sup>i</sup>* = *ε* 0 *<sup>i</sup>* = *ε* 00 *i* , the undamaged part of the rock obeys Hooke's law, and its stress is:

$$
\sigma\_i' = E\varepsilon\_i + \mu(\sigma\_j' + \sigma\_k') \tag{23}
$$

Equation (23) can be written as:

$$\begin{cases} \sigma\_1' = E\varepsilon\_1 + \mu(\sigma\_2' + \sigma\_3') \\ \sigma\_2' = E\varepsilon\_2 + \mu(\sigma\_1' + \sigma\_3') \\ \sigma\_3' = E\varepsilon\_3 + \mu(\sigma\_1' + \sigma\_2') \end{cases} \tag{24}$$

The damage of the rock micro-element can be defined as the result of the reduction of stiffness caused by the change of physical properties of undamaged micro-element. The stress of damaged micro-element and undamaged the micro-element under the external load is related to the stiffness, and *γ* is defined as the damage correction coefficient [33]. Based on the rock damage mechanism, the damage models of rock micro-elements in different states can be established.

(1) When the rock is not damaged:

$$
\sigma\_{\mathbf{i}}^{\prime} = \sigma \quad , \quad \sigma\_{\mathbf{i}}^{\prime\prime} = 0 \tag{25}
$$

(2) When random damage occurs to rock:

$$
\sigma = \sigma\_i'(1 - D) + \sigma\_i''D \tag{26}
$$

From the stress distribution relationship between damaged micro-elements and undamaged micro-elements, it can be concluded that:

$$
\sigma\_i'' = \gamma \cdot \sigma\_i' \tag{27}
$$

Substituting to formula (26) to obtain:

$$
\sigma\_i = \sigma\_i'(1 - D + \gamma D) \tag{28}
$$

(3) When the rock is completely damaged:

$$
\sigma\_{\rm i}^{\prime\prime} = \sigma \,\,\sigma\_{\rm i}^{\prime} = 0 \tag{29}
$$

Substituting Equation (23) into Equation (28) yields:

$$
\sigma\_1 = E\varepsilon\_1(1 - D + \gamma D) + \mu(\sigma\_2 + \sigma\_3) \tag{30}
$$

Equation (30) is the damage constitutive equation of rock. To establish the damage constitutive model of rock, the damage variable must first be determined.

## *5.2. Damage Evolution Equation of Siltstone*

In each stage of rock damage and failure, the thermodynamic change of the system is regarded as a process in which the equilibrium is broken and then reconstructed. When the rock is damaged and deformed, the equilibrium state of the system is broken. When the energy dissipation tends to be stable, the system reaches a new equilibrium state. The energy dissipation process of rock damage and deformation conforms to the minimum energy dissipation principle. According to the minimum energy dissipation principle, all energy dissipation processes occur along the minimum energy dissipation path under corresponding constraints. It illustrates that the instantaneous energy dissipation rate is at the minimum at any time in the energy dissipation process. In the elastic damage model, it is assumed that the irreversible strain caused by damage is the only energy dissipation mechanism in the rock failure process [23]. The energy consumption rate of rock is defined as:

$$
\varphi = \sigma\_i \mathfrak{e}\_i \tag{31}
$$

where, *ϕ* is the rock energy consumption rate, *σ<sup>i</sup>* is the nominal stress of rock micro element, and • *εi* is the irreversible strain rate caused by damage.

The three-dimensional constitutive relationship of rock can be obtained from Equation (30) as follows:

$$\begin{cases} \begin{array}{l} \varepsilon\_{1} = \frac{\sigma\_{1} - \mu(\sigma\_{2} + \sigma\_{3})}{[1 - D(t) + \gamma D(t)]E} \\ \varepsilon\_{2} = \frac{\sigma\_{2} - \mu(\sigma\_{1} + \sigma\_{3})}{[1 - D(t) + \gamma D(t)]E} \\ \varepsilon\_{3} = \frac{\sigma\_{3} - \mu(\sigma\_{1} + \sigma\_{2})}{[1 - D(t) + \gamma D(t)]E} \end{array} \end{cases} \tag{32}$$

where, *D*(*t*) is the damage variable at *t* time.

The irreversible strain rate caused by the damage variable is:

$$\begin{cases} \begin{aligned} \mathcal{E}\_{1}^{\bullet} &= \frac{(1-\gamma)D'(t)[\sigma\_{1}-\mu(\sigma\_{2}+\sigma\_{3})]}{\left[1-D(t)+\gamma D(t)\right]^{2}E} \\ \mathcal{E}\_{2}^{\bullet} &= \frac{(1-\gamma)D'(t)[\sigma\_{2}-\mu(\sigma\_{1}+\sigma\_{3})]}{\left[1-D(t)+\gamma D(t)\right]^{2}E} \\ \mathcal{E}\_{3}^{\bullet} &= \frac{(1-\gamma)D'(t)[\sigma\_{3}-\mu(\sigma\_{1}+\sigma\_{2})]}{\left[1-D(t)+\gamma D(t)\right]^{2}E} \end{aligned} \tag{33}$$

By substituting Equation (33) into Equation (31), the energy consumption rate of rock is found as follows:

$$\varphi = \frac{(1-\gamma)D'(t)}{\left[1 - D(t) + \gamma D(t)\right]^2 E} \left[\sigma\_1^2 + \sigma\_2^2 + \sigma\_3^2 - 2\mu(\sigma\_1\sigma\_2 + \sigma\_2\sigma\_3 + \sigma\_1\sigma\_3)\right] \tag{34}$$

For a conventional triaxial compression test, *σ*<sup>1</sup> > *σ*<sup>2</sup> = *σ*3, Equation (34) can be written as:

$$\varphi = \frac{(1-\gamma)D'(t)}{\left[1 - D(t) + \gamma D(t)\right]^2 E} \left[\sigma\_1^2 + 2(1-\mu)\sigma\_3^2 - 4\mu\sigma\_1\sigma\_3\right] \tag{35}$$

*5.3. Establishment of Damage Constitutive Model*

Substitute Equation (21) into Equation (12) to obtain [23]:

$$\frac{\partial(\varrho + \lambda F(\sigma\_1, \sigma\_3))}{\partial \sigma\_i} = 0 \ (i = 1, 3) \tag{36}$$

$$\begin{cases} \frac{\partial (\varrho + \lambda F(\sigma\_1 \rho\_3))}{\partial \sigma\_1} = 0\\ \frac{\partial (\varrho + \lambda F(\sigma\_1 \rho\_3))}{\partial \sigma\_3} = 0 \end{cases} \tag{37}$$

where, *∂ϕ ∂σ*1 = (1−*γ*)*D*<sup>0</sup> 1−*D*+*γD* · 2*ε*1, *∂ϕ ∂σ*3 = (1−*γ*)*D*<sup>0</sup> 1−*D*+*γD* · 4*ε*3, *∂F ∂σ*1 = (1−sin *ϕ*) 2 2*G* cos<sup>2</sup> *ϕ σ*<sup>1</sup> + <sup>−</sup> <sup>1</sup> 2*G σ*<sup>3</sup> + *a*, *∂F ∂σ*3 = <sup>−</sup> <sup>1</sup> 2*G σ*<sup>1</sup> + (1+sin *ϕ*) 2 2*G* cos<sup>2</sup> *ϕ σ*<sup>3</sup> + 2*a*.

The damage evolution equation is derived as follows:

$$D(t) = \frac{1}{1-\gamma} - \frac{1}{1-\gamma} \exp\left(\frac{\lambda \cdot H}{R} t + \mathcal{C}\_0\right) \tag{38}$$

where, *A* = (1−sin *ϕ*) 2 2*G* cos<sup>2</sup> *ϕ* , *<sup>B</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup> 2*G* , *C* = (1+sin *ϕ*) 2 2*G* cos<sup>2</sup> *ϕ* , *H* = *B* <sup>2</sup> <sup>−</sup> *AC σ*<sup>3</sup> + (*B* − 2*A*)*a*, *R* = 2*Bε*<sup>1</sup> − 4*Aε*3, *λ* and *C*<sup>0</sup> are parameters related to material properties.

Under a constant loading rate and loading path, the axial strain of the rock sample is directly proportional to time. Assuming that *λ* · *t* = *λ* ∗ · *ε*1, the damage evolution equation is simplified as:

$$D(t) = \frac{1}{1-\gamma} - \frac{1}{1-\gamma} \exp\left(\frac{\lambda^\* \cdot H}{R} \varepsilon\_1 + \mathcal{C}\_0\right) \tag{39}$$

### *5.4. Parameter Identification of Constitutive Model*

According to the triaxial compressive stress–strain curve of rock, when the rock is completely damaged (*D =* 1), Formula (30) is written as:

$$
\sigma\_1^r = E \varepsilon\_1^r \gamma + \mu (\sigma\_2 + \sigma\_3) \tag{40}
$$

$$\gamma = \frac{\sigma\_1^r - \mu(\sigma\_2 + \sigma\_3)}{E\varepsilon\_1^r} \tag{41}$$

Two boundary conditions can be determined from the extreme value characteristics of the rock total stress–strain curve. The peak stress of the rock is *σsc* and the corresponding strain is *εsc*. The two boundary conditions are as follows:

where *ε*<sup>1</sup> = *εsc*, *σ*<sup>1</sup> = *σsc*, d*σ*1/d*ε*<sup>1</sup> = 0.

Substitute conditions into Formula (30) to obtain:

$$
\sigma\_{\rm sc} = E \varepsilon\_{\rm sc} \exp\left(\frac{\lambda^\* \cdot H}{R} \varepsilon\_{\rm sc} + \mathbb{C}\_0\right) + 2\mu \sigma\_3 \tag{42}
$$

In the triaxial compression test, the peak stress *σsc* and peak strain *εsc* of rock are not the peak strength *σ* 0 *sc* and corresponding strain *ε* 0 *sc* of the rock deviatoric stress–strain curve, because the starting point of the test curve provides that the rock is subjected to deviatoric stress, ignoring the initial strain of the rock under hydrostatic pressure, therefore, the relationship between the theoretical peak value and the test peak value of rock stress–strain curve is as follows:

$$
\sigma\_{\text{sc}} = \sigma\_{\text{sc}}' + \sigma\_{\text{3}} \tag{43}
$$

$$
\varepsilon\_{\rm sc} = \varepsilon\_{\rm sc}^{\prime} + \varepsilon\_{\rm c} \tag{44}
$$

The initial strain can be obtained from the first equation in formula (24):

$$
\varepsilon\_{\mathcal{L}} = \frac{\sigma\_{\mathcal{3}}(1 - 2\mu)}{E} \tag{45}
$$

The values of *λ* ∗ and *C*<sup>0</sup> are as follows:

$$\lambda^\* = \frac{\left(2B\varepsilon\_{1(\text{sc})} - 4A\varepsilon\_{3(\text{sc})}\right)^2}{4AH\varepsilon\_{1(\text{sc})}\varepsilon\_{3(\text{sc})}} = \frac{R^2}{4AH\varepsilon\_{1(\text{sc})}\varepsilon\_{3(\text{sc})}}\tag{46}$$

$$\mathcal{C}\_{0} = \ln\left[\frac{\sigma\_{1(\text{sc})} - 2\mu\sigma\_{3}}{E\varepsilon\_{1(\text{sc})}}\right] - \frac{B\varepsilon\_{1(\text{sc})} - 2A\varepsilon\_{3(\text{sc})}}{2A\varepsilon\_{3(\text{sc})}} = \ln\left[\frac{\sigma\_{1(\text{sc})} - 2\mu\sigma\_{3}}{E\varepsilon\_{1(\text{sc})}}\right] - \frac{R}{4A\varepsilon\_{3(\text{sc})}}\tag{47}$$
