# Study on Damage Evolution Model of Sandstone under Triaxial Loading and Postpeak Unloading Considering Nonlinear Behaviors

Authors: Tianbai Zhou; Yueping Qin; Jian Cheng; Xiaoyu Zhang; Qiufeng Ma
Journal: Geofluids
 (2021)

---

## Abstract

Conventional triaxial loading and unloading tests were carried out on sandstone samples in the Zigong area, of Sichuan Province, China. The changes in the elastic modulus of the unloading curves under different confining pressures were calculated, and the evolution law of the nonlinear properties of rock was analyzed. The results show that the rock is subjected to nonlinear damage during initial compaction, the elastic phase, destruction, and postpeak unloading. Moreover, the nonlinear behaviors of rock are restrained by the confining pressures. On this basis, a nonlinear stress-strain relationship affected by the average stress is proposed to describe nonlinear behaviors in the initial compaction stage. According to the test data, the evolution laws of various energies inside the rock during loading and unloading cycles are obtained. The results show that the external work is transformed into elastic energy and damage dissipated energy. Based on the energy analysis, the energy balance equation is established according to the law of energy conservation. By deriving the energy balance equation, the damage evolution equation of sandstone under triaxial loading is solved to establish a continuous constitutive model. The calculation results of the model are compared with the test results from two aspects of loading and postpeak unloading. The comparison results show that the proposed model, which reflects the whole stress-strain process and nonlinear properties of rock, could also describe the stress-strain relationship at the postpeak unloading stage to some extent.


---

## Body

## 1. Introduction

In the process of underground mining, the surrounding rock of the working space is always in a complex stress state. The damage and destruction of surrounding rock directly threaten the production efficiency and personnel safety of the working space [1]. Rock is an elastic-plastic material with complex mechanical properties. Research on the damage evolution, stress-strain relationship, rock strength, postpeak bearing capacity, nonlinear deformation characteristics, and rheological properties of rock has always been an important exploration direction in the field of rock mechanics [2–5]. The study of the damage evolution model is mainly based on the theory of continuum mechanics. The damage variable of rock is defined by the weakening phenomenon of some mechanical parameters in macromodeling ideas, and the damage evolution model can be deduced [6–9]. Some studies from the microscopic point of view do not build fracture models to describe constitutive relationships for rock but rather define damage variables based on microscopic parameters such as fractures and defects inside the rock to characterize the destruction process. At present, the establishment of a damage evolution model is a widely used method to reflect the mechanical behavior of rock to provide help for engineering calculations.Many scholars have done much work on damage evolution models [10–15]. Unteregger et al. [16] proposed a damage-plastic coupling model based on continuum theory to describe the damage process of intact rock and established a constitutive model suitable for complex three-dimensional stress states. Yang et al. [17–20] introduced the concept of the yield stress ratio and proposed a damage softening model that combined the M-C yield criterion to describe the whole stress-strain curve. Based on continuous damage mechanics, Poulet et al. [21] established the damage model of a coupled thermal field, fluid field, and stress field and implemented it with numerical software. Ma et al. [22–25] studied the effect of defects on the mechanical behaviors of rock under static or dynamic loads and explored the effect of sample size on the damage evolution of granite by acoustic emission tests and suggested that the height/diameter ratio (H/D) determines the failure form of rock samples. Some studies have analyzed the damage mechanisms from the perspective of rock fracture and proposed optimized rock damage models to provide a basis for engineering calculations and disaster prevention [22, 24, 26, 27]. Cao et al. [28–30] assumed that rock is composed of voids, damaged structures, and undamaged structures. Based on the statistical damage method, a statistical damage constitutive model that can reflect strain softening and hardening behaviors under triaxial compression conditions was established and compared with the experimental observation results. Ma et al. [31] analyzed the mechanical behaviors of porous rock during loading and unloading cycles and established a Poisson’s ratio model for the unloading process according to the evolution law of circumferential strain. The D-P plastic model was introduced to describe the plastic hardening characteristics of porous rock so that a damage model related to equivalent plastic strain was established. Zhang et al. [32] proposed associating the state functions of elastic strain and plastic variables with damage. On this basis, the elastic-plastic damage coupling model for brittle rock was established to simulate specific engineering problems. Many modeling ideas about damage evolution were put forward by these studies, which are also of great significance for engineering calculations [33, 34]. However, from the following two perspectives, there is still room for improvement and optimization in the study. In some models, the gradual evolution of rock damage is ignored, and a two-point theory is adopted; that is, yield or no yield. These models do not consider the gradual damage process of intact rock from the undamaged state to the damaged state. There may be some deficiencies in analyzing the dynamic disasters of coal and rock masses by these models. For example, some statistical damage models describe the damage evolution process by a subsection function, which is equivalent to artificially adding a turning point to define the evolution of rock from the undamaged state to the damaged state. Using a segmented model, there may be some limitations in describing the continuous process of rock damage and failure. On the other hand, energy transformation and dissipation are inevitable in rock damage and failure. Coal and rock dynamic disasters can also be regarded as processes of intense energy release. Therefore, establishing a rock damage evolution model from the perspective of energy is conducive to revealing the damage evolution mechanism in essence, and the modeling idea around energy is more in line with objective physical significance.In recent years, the exploration of rock constitutive models and damage mechanisms from the perspective of energy has been greatly developed, and researchers have done much work [19, 20, 35–40]. Xie et al. [41, 42] explored the energy mechanism in the process of rock deformation and discussed the relationship between energy release and rock failure. On this basis, energy dissipation and damage variables were calculated, and a damage evolution equation was established. Liu et al. [43] introduced reflecting the influence of energy dissipation on the compaction process by modifying the energy release rate. Combined with the Weibull distribution model, a statistical damage model considering energy evolution under triaxial compression was proposed by Liu. Based on thermodynamics theory, Qin et al. [44, 45] interpreted the process of rock damage and failure by an energy evolution mechanism. The process of energy conversion and dissipation was derived under the framework of irreversible thermodynamics, and a rock constitutive model under uniaxial loading was established based on the first law of thermodynamics. Liu et al. [46] carried out a uniaxial cyclic loading test and proposed the concept of the compaction coefficient to modify the damage equation. They established an energy dissipation model of rock, and the model results can reflect the loading and unloading stress-strain curves before the peak stress. At present, there are still some deficiencies in the study of constitutive models based on energy, such as the nonlinear properties and the stress-strain relationship under postpeak unloading, which are also main objectives in this work [23, 25, 47].In this work, the proposed model was established based on the damage evolution constrained by the energy relationship, which could reflect the irreversibility of energy dissipation and damage development. Nonlinear properties of rock that were not considered in some previous models were taken as the main research objective, and a model was proposed to describe nonlinear mechanical behaviors. A continuous constitutive equation that could represent the whole stress-strain process was established. Finally, we explored the simulation of postpeak unloading curves by the proposed model.


## 2. Experimental Process

In this paper, sandstone samples were taken from the Zigong area, Sichuan Province, China. Zigong is located in the eastern part of the Sichuan Basin, which belongs to the Yangzi stratum of the South China stratum according to the lithostratigraphic division in 1997. The accumulation thickness of sandstone in the Zigong area has reached 3000-4000 meters since the Mesozoic period. Zigong sandstone is composed of silicate minerals, the main components of which are feldspar and quartz. After cutting and grinding in the laboratory, the samples were processed into standard cylindrical samples with a size ofφ50×100mm according to the Chinese standard “Standard Test Method for Engineering Rock” (GB/T 50266-2013), as shown in Figure 1. The average density and porosity of the sandstone samples are 2.46 g cm-3 and 14.7%, respectively. In the test, all 36 sandstone samples were taken from the parts with good integrity in the rock block. In the process of sample preparation, intact samples were selected to ensure the quality of sandstone samples.Figure 1
Sandstone samples.The uniaxial compression test, conventional triaxial compression test, and conventional triaxial cyclic loading and unloading test were carried out with a RLJW-2000 electrohydraulic servo test system in the experiment. The confining pressure was set to 5 MPa, 10 MPa, 20 MPa, and 30 MPa, and the loading and unloading rates of the test system were 0.005 mm/s. The test system is shown in Figure2(a), and the failure forms of some samples are shown in Figure 2(b). First, the peak stress under different confining pressures was obtained by a conventional triaxial compression test. Then, a cyclic loading test was carried out before the peak stress was reached with unloading points at 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90% of the peak stress. The postpeak unloading includes two unloading cycles during the softening stage and no less than two unloading cycles in the residual stage. To ensure the accuracy of the test, each test was repeated multiple times.Figure 2
(a) Test system and (b) failure form of sandstone samples.
(a)(b)The mechanical properties of sandstone under different confining pressures obtained in the test are shown in Table1. In the table, the initial elastic modulus E0 can be determined by calculating the tangent modulus of the elastic stage of the outer contour line of the stress-strain curves. The residual strength can be determined by averaging the stress strength under the residual stage of the cyclic loading and unloading curves. In this paper, triaxial cyclic loading test data are mainly used. Because there is too much original data, data smoothing was applied to reflect the evolution law of the stress-strain curve more clearly. The stress-strain curves under different confining pressures are shown in Figure 3.Table 1
Mechanical properties of sandstone.
Confining pressures, MPaInitial elastic modulusE0, GPaPeak strength, MPaResidual strength, MPa514.812630.21018.9161.160.52020.6201.385.83022.2233.2105.8Figure 3
Stress-strain curves under cyclic loading and unloading test.


## 3. Study on Nonlinear Behaviors

### 3.1. Analysis of Nonlinear Properties

To explore the nonlinear properties of rock, the evolution law and elastic modulus of unloading curves under different confining pressures were analyzed and plotted in graphs, as shown in Figures4–7. The tangent modulus of the approximately straight line segment in each unloading curve was calculated to characterize the elastic modulus of rock. In each figure, the left subfigure reflects the evolution of unloading curves. The black curve represents the outer contour line of the stress-strain curve, and the brown solid line and dotted line represent the prepeak and postpeak unloading curves, respectively. The elastic modulus of the unloading curve is shown in the right subfigure, in which the abscissa is the unloading order and the ordinate is the elastic modulus. For example, it can be seen from Figure 4 that under a confining pressure of 5 MPa, a total of 11 loading and unloading cycles were carried out 8 times at the prepeak stage and 3 times at the postpeak stage.Figure 4
Unloading curves and evolution of elastic modulus under 5 MPa confining pressure.Figure 5
Unloading curves and evolution of elastic modulus under 10 MPa confining pressure.Figure 6
Unloading curves and evolution of elastic modulus under 20 MPa confining pressure.Figure 7
Unloading curves and evolution of elastic modulus under 30 MPa confining pressure.The evolution process of the elastic modulus during the whole stress-strain process can be seen in Figures4–7. In the compression and elastic stages, the axial stress-strain curve is concave. In this stage, the primary pores inside the rock are gradually compacted, resulting in an increase in the elastic modulus, which reaches the maximum value before the yield stage. With gradual damage, plastic deformation occurs inside the rock, and the stress-strain curve gradually appears as a convex curve. In the postpeak stage, the rock was destroyed, and the elastic modulus was significantly reduced. In the residual stage, the elastic modulus of unloading curves that slightly change is maintained at a constant value. It can also be observed that each unloading curve shows obvious nonlinear characteristics regardless of the prepeak or postpeak stage. Specifically, the elastic modulus of rock is not a constant. This is because the primary pores inside the rock are compacted during the loading process, which results in an increase in the elastic modulus. During the unloading process, the primary pores gradually become loose from the compaction state, decreasing the elastic modulus. Compared with the unloading curves in Figures 4–7, we see that different confining pressure conditions will affect the degree of nonlinear properties. With increasing confining pressure, the pore state of rock becomes more compact, so the nonlinear properties become nonobvious. In contrast, the pore state of rock becomes relatively loose with decreasing confining pressure, which leads to a more apparent nonlinear property of rock.The definitions of the following parameters are introduced. The maximum value of the elastic modulus of the unloading curve is defined as the peak elastic modulus. The physical meaning of the peak elastic modulus is the maximum value of the elastic modulus after the compaction stage under constant confining pressure, which is represented byEm. The concept of the residual elastic modulus is proposed to reflect the elastic modulus of rock in the residual stage, which is represented by Er. To reflect the relationship between the initial elastic modulus and residual elastic modulus of rock, the elastic modulus loss coefficient is introduced, which can be expressed as(1)α=E0−ErE0,where α is the elastic modulus loss coefficient, which represents the loss rate of elastic modulus after complete damage of rock. E0 is the initial elastic modulus. The value of elastic modulus loss coefficient is between 0 and 1. The loss of elastic modulus increases with α. The maximum value of α is 1, which means the complete loss of elastic modulus. The minimum value of α is 0, which means that elastic modulus of rock is still the initial elastic modulus without any loss.The values of these parameters can be obtained by experiments, as shown in Table2. It can be seen from the table that the initial elastic modulus, peak elastic modulus, and residual elastic modulus are all positively related to confining pressure. The elastic modulus loss coefficient is negatively correlated with confining pressure. The larger the confining pressure is, the smaller the loss of elastic modulus in the residual stage. Therefore, the residual strength increases with the confining pressure.Table 2
Related parameters of elastic modulus.
Confining pressures, MPaInitial elastic modulusE0, GPaPeak elastic modulusEm, GPaResidual elastic modulusEr, GPaElastic modulus loss coefficientα514.819.29.10.381017.723.712.60.292020.125.216.20.193022.229.1200.09


### 3.2. The Establishment of Nonlinear Model

According to the previous analysis, nonlinear behaviors should be an inherent property of rock and are related to the stress state. On this assumption, the concept of average stress is introduced to reflect the stress state of rock. The average stress is defined as the average value of the axial pressure and the confining pressures, which is an important indicator of the stress level of the rock. For example, in the loading process, the average stress increases gradually, and the rock becomes dense, which leads to an increase in the elastic modulus. During the unloading process, the average stress decreases gradually such that the rock becomes relatively loose, so the elastic modulus decreases gradually. For example, at the same stress curve stage, the elastic modulus will increase with the confining pressure condition. In other words, the elastic modulus will increase with the average stress, which is consistent with general physical laws. On this basis, an expression of the elastic modulus is established to describe the nonlinear properties of rock in this paper. It can be assumed that the elastic modulus of rock is positively correlated with the average stress, and the expression of the elastic modulus is mainly based on the functiony=−exp−kx. Drawing on previous studies and our conjecture [31], the expression of the elastic modulus can be expressed as follows:(2)E=Em1−exp−kσm,where E represents elastic modulus of rock. k is the nonlinear coefficient. σm is the average stress, which can be expressed as σm=σ1+σ2+σ3/3. σii=1,2,3 represents the principal stress in three directions.As shown in Equation (2), when σm approaches infinity, the value of elastic modulus E approaches to peak elastic modulus Em. The nonlinear coefficient k is a material parameter which does not change with the loading condition. The value of k is 0.05, which directly affects the degree of nonlinear properties. When the value of k is small, the nonlinear behaviors of rock are significant. On the contrary, the stress-strain curve tends to be linear.Combining Equation (2), the incremental relationship of axial stress and axial strain can be established:(3)dσ1=Edε1=Em1−exp−kσmdε1,where dσ1 is the increment of axial stress and dε1 is the increment of axial strain.According to Equation (3), the axial stress-strain relationship of rock under different confining pressures can be calculated without considering the damage effect. In the compression stage and elastic stage of the stress-strain curve, rock damage is usually negligible. Therefore, the model calculation results are compared with the test results under different confining pressures, as shown in Figure 8. It can be seen from Figure 8 that the model results are roughly consistent with the test results, and the rationality of the elastic modulus expression has been preliminarily verified. The model can reflect the nonlinear properties of rock.Figure 8
Comparison of model results and test results at initial stage of axial stress-strain curves under different confining pressures.


## 3.1. Analysis of Nonlinear Properties

To explore the nonlinear properties of rock, the evolution law and elastic modulus of unloading curves under different confining pressures were analyzed and plotted in graphs, as shown in Figures4–7. The tangent modulus of the approximately straight line segment in each unloading curve was calculated to characterize the elastic modulus of rock. In each figure, the left subfigure reflects the evolution of unloading curves. The black curve represents the outer contour line of the stress-strain curve, and the brown solid line and dotted line represent the prepeak and postpeak unloading curves, respectively. The elastic modulus of the unloading curve is shown in the right subfigure, in which the abscissa is the unloading order and the ordinate is the elastic modulus. For example, it can be seen from Figure 4 that under a confining pressure of 5 MPa, a total of 11 loading and unloading cycles were carried out 8 times at the prepeak stage and 3 times at the postpeak stage.Figure 4
Unloading curves and evolution of elastic modulus under 5 MPa confining pressure.Figure 5
Unloading curves and evolution of elastic modulus under 10 MPa confining pressure.Figure 6
Unloading curves and evolution of elastic modulus under 20 MPa confining pressure.Figure 7
Unloading curves and evolution of elastic modulus under 30 MPa confining pressure.The evolution process of the elastic modulus during the whole stress-strain process can be seen in Figures4–7. In the compression and elastic stages, the axial stress-strain curve is concave. In this stage, the primary pores inside the rock are gradually compacted, resulting in an increase in the elastic modulus, which reaches the maximum value before the yield stage. With gradual damage, plastic deformation occurs inside the rock, and the stress-strain curve gradually appears as a convex curve. In the postpeak stage, the rock was destroyed, and the elastic modulus was significantly reduced. In the residual stage, the elastic modulus of unloading curves that slightly change is maintained at a constant value. It can also be observed that each unloading curve shows obvious nonlinear characteristics regardless of the prepeak or postpeak stage. Specifically, the elastic modulus of rock is not a constant. This is because the primary pores inside the rock are compacted during the loading process, which results in an increase in the elastic modulus. During the unloading process, the primary pores gradually become loose from the compaction state, decreasing the elastic modulus. Compared with the unloading curves in Figures 4–7, we see that different confining pressure conditions will affect the degree of nonlinear properties. With increasing confining pressure, the pore state of rock becomes more compact, so the nonlinear properties become nonobvious. In contrast, the pore state of rock becomes relatively loose with decreasing confining pressure, which leads to a more apparent nonlinear property of rock.The definitions of the following parameters are introduced. The maximum value of the elastic modulus of the unloading curve is defined as the peak elastic modulus. The physical meaning of the peak elastic modulus is the maximum value of the elastic modulus after the compaction stage under constant confining pressure, which is represented byEm. The concept of the residual elastic modulus is proposed to reflect the elastic modulus of rock in the residual stage, which is represented by Er. To reflect the relationship between the initial elastic modulus and residual elastic modulus of rock, the elastic modulus loss coefficient is introduced, which can be expressed as(1)α=E0−ErE0,where α is the elastic modulus loss coefficient, which represents the loss rate of elastic modulus after complete damage of rock. E0 is the initial elastic modulus. The value of elastic modulus loss coefficient is between 0 and 1. The loss of elastic modulus increases with α. The maximum value of α is 1, which means the complete loss of elastic modulus. The minimum value of α is 0, which means that elastic modulus of rock is still the initial elastic modulus without any loss.The values of these parameters can be obtained by experiments, as shown in Table2. It can be seen from the table that the initial elastic modulus, peak elastic modulus, and residual elastic modulus are all positively related to confining pressure. The elastic modulus loss coefficient is negatively correlated with confining pressure. The larger the confining pressure is, the smaller the loss of elastic modulus in the residual stage. Therefore, the residual strength increases with the confining pressure.Table 2
Related parameters of elastic modulus.
Confining pressures, MPaInitial elastic modulusE0, GPaPeak elastic modulusEm, GPaResidual elastic modulusEr, GPaElastic modulus loss coefficientα514.819.29.10.381017.723.712.60.292020.125.216.20.193022.229.1200.09


## 3.2. The Establishment of Nonlinear Model

According to the previous analysis, nonlinear behaviors should be an inherent property of rock and are related to the stress state. On this assumption, the concept of average stress is introduced to reflect the stress state of rock. The average stress is defined as the average value of the axial pressure and the confining pressures, which is an important indicator of the stress level of the rock. For example, in the loading process, the average stress increases gradually, and the rock becomes dense, which leads to an increase in the elastic modulus. During the unloading process, the average stress decreases gradually such that the rock becomes relatively loose, so the elastic modulus decreases gradually. For example, at the same stress curve stage, the elastic modulus will increase with the confining pressure condition. In other words, the elastic modulus will increase with the average stress, which is consistent with general physical laws. On this basis, an expression of the elastic modulus is established to describe the nonlinear properties of rock in this paper. It can be assumed that the elastic modulus of rock is positively correlated with the average stress, and the expression of the elastic modulus is mainly based on the functiony=−exp−kx. Drawing on previous studies and our conjecture [31], the expression of the elastic modulus can be expressed as follows:(2)E=Em1−exp−kσm,where E represents elastic modulus of rock. k is the nonlinear coefficient. σm is the average stress, which can be expressed as σm=σ1+σ2+σ3/3. σii=1,2,3 represents the principal stress in three directions.As shown in Equation (2), when σm approaches infinity, the value of elastic modulus E approaches to peak elastic modulus Em. The nonlinear coefficient k is a material parameter which does not change with the loading condition. The value of k is 0.05, which directly affects the degree of nonlinear properties. When the value of k is small, the nonlinear behaviors of rock are significant. On the contrary, the stress-strain curve tends to be linear.Combining Equation (2), the incremental relationship of axial stress and axial strain can be established:(3)dσ1=Edε1=Em1−exp−kσmdε1,where dσ1 is the increment of axial stress and dε1 is the increment of axial strain.According to Equation (3), the axial stress-strain relationship of rock under different confining pressures can be calculated without considering the damage effect. In the compression stage and elastic stage of the stress-strain curve, rock damage is usually negligible. Therefore, the model calculation results are compared with the test results under different confining pressures, as shown in Figure 8. It can be seen from Figure 8 that the model results are roughly consistent with the test results, and the rationality of the elastic modulus expression has been preliminarily verified. The model can reflect the nonlinear properties of rock.Figure 8
Comparison of model results and test results at initial stage of axial stress-strain curves under different confining pressures.


## 4. Energy Balance Equation

In the energy analysis of the rock damage evolution process, the rock can be regarded as a thermodynamic system, and the rock surface can be regarded as the boundary of the thermodynamic system. According to the first law of thermodynamics, which is the law of energy conservation, the energy entering the rock system minus the dissipated energy of the system is equal to the increase in the stored energy inside the rock system. For the rock system, the energy input to rock is external work. The consumed energy of the rock system is the energy dissipated by the rock damage process. The difference between these two is converted into elastic energy stored inside the rock. Figure9 shows the energy relationship during the loading and unloading cycle of rock. The black curve represents the loading curve, and the red curve represents the unloading curve. It can be seen from Figure 9 that during the loading process, the test system transmits energy to the rock sample by mechanical work, and that the integral area OAC of the loading curve is the energy input by external work. In the unloading process, the elastic energy accumulated in the rock is released, so the integral area ABC of the unloading curve represents elastic energy. The damage dissipated energy is represented by the integral area OAB between the loading and unloading curves.Figure 9
Energy relationship during loading and unloading process.In conclusion, the energy balance equation of rock can be expressed as follows:(4)δW=dE~+δA,where W represents the external work. E~ represents the elastic energy of rock system. A represents the damage dissipated energy. δ and d are differential operators of process quantity and state quantity, respectively.For the unit volume of rock, the energy balance equation can be expressed as(5)δw=de+δa,where δw represents the external work to unit volume rock, de represents the increment of elastic energy, and δa represents the increment of damage dissipated energy.As shown in Figure9, the abscissa is the strain, the ordinate is the stress, and the integral area of curve is the energy density. The integral area OAC of the loading curve is input energy density of external work. The integral area ABC of the unloading curve is the elastic energy density. The integral area OAB between the loading and unloading curves is the damage dissipated energy density. Therefore, the elastic energy density and the damage dissipated energy density at any unloading point can be obtained by the mathematical relationship between the integral area of the loading and unloading curves. The input energy density of external work, elastic energy density, and damage dissipated energy density under each unloading point with different confining pressures both can be obtained by calculating the integral area of stress curve. With loading, the evolution of the energy density is shown in Figure 10. The evolution curves of the ratio of elastic energy and damage dissipated energy to the input energy by external work are shown in Figure 11.Figure 10
Energy evolution curves under different confining pressures ((a) 5 MPa, (b) 10 MPa, (c) 20 MPa, and (d) 30 MPa).
(a)(b)(c)(d)Figure 11
Energy ratio evolution curves of elastic energy and damage dissipated energy under different confining pressures ((a) 5 MPa, (b) 10 MPa, (c) 20 MPa, and (d) 30 MPa).
(a)(b)(c)(d)Figures10 and 11 show that the evolution trends of rock energy density during cyclic loading and unloading cycles are similar under different confining pressure conditions. However, there are differences in the evolution law of the energy density in different stages of the stress-strain curve. During the compaction stage, the primary microcracks inside the rock are compressed and closed, and the energy inputted by external work and the elastic energy increase slowly. In the elastic stage, the growth of these two types of energy accelerates due to the increase in rock rigidity, resulting in the gradual increase in energy storage capacity. During the elastic stage, most of the energy inputted by external work is transformed into elastic energy, which is stored inside the rock. At this time, the damage dissipated energy accounts for a small proportion. In the compaction stage and elastic stage, the proportion of elastic energy is relatively large, and some irregular data may be due to measurement errors, which indicates that most of the inputted energy is converted into elastic energy and stored inside the rock. Before the yield stage, the rock is roughly intact, so elastic energy is dominant. When the rock yields, due to the increase in damage, new cracks are produced inside the rock, and the energy density of the damage dissipated energy increases gradually. In general, the elastic energy density reaches its maximum around the peak stress. During the yield stage, the ratio of damage dissipated energy increases perceptibly, indicating that more energy inputted by external force is dissipated by damage. In the residual stage, the rock is destroyed, and the main fracture is formed. At this stage, the dissipated energy density increases rapidly, and the elastic energy density decreases rapidly and finally tends to be stable.The characteristic values of rock energy under different confining pressures are shown in Table3. In the whole stress-strain process, the maximum elastic energy density is defined as the peak elastic energy density, and the average elastic energy density in the residual stage is defined as the residual elastic energy density. Table 3 shows that the four energy eigenvalues in the table all show noticeable confining pressure effects, which increase with increasing confining pressure. According to the analysis of Figure 10 and Table 2, the input energy required for rock failure increases with increasing confining pressure, and the stored elastic energy and damage dissipated energy also increase.Table 3
Characteristic values of rock energy under different confining pressures.
Confining pressures, MPaPeak elastic energy density, MJ m-3Residual elastic energy density, MJ m-3Dissipated energy density at peak stress, MJ m-3Inputted energy density at peak stress, MJ m-350.320.090.220.55100.440.130.480.93200.560.170.691.15300.880.220.881.75


## 5. Study on Damage Evolution Model

### 5.1. Definition of Damage Variable

In damage mechanics, rock can be regarded as a continuous medium. The evolution process of the mechanical properties of rock under various stress conditions can be studied by introducing damage variables. In short, the damage variable can be understood as a measure of the degree of decay of the mechanical properties. When the material is damaged, there will be a series of changes. For rock, there will be a decrease in the elastic modulus, generation of plastic strain, increase in volume, decrease in density, generation of acoustic emission and electromagnetic radiation, change in temperature, and so on. Therefore, the damage variable can be defined by the changes in these physical quantities. The damage variable is usually defined as the decay of the elastic modulus for rock in domestic and foreign studies, which can be expressed as(6)D=1−EE0,where D is the damage variable of rock.According to the definition of damage, elastic modulus of damaged rock can be expressed as follows:(7)E=1−DE0.However, Equation (7) is more suitable for describing the elastic damage behavior. To describe the rock damage process with both elasticity and plasticity, the expression needs to be modified. For example, the evolution process of the elastic modulus in the whole stress-strain curve has been analyzed in detail. It can be seen from Figure 7 that in the residual stage, there is still a residual elastic modulus. When the value of the damage variable is 1, the value of the elastic modulus calculated according to Equation (7) is 0, which does not conform to the test results. Therefore, the expression of the elastic modulus needs to be modified:(8)E=1−αDE0,where α is the elastic modulus loss coefficient. The physical meaning and value selection of loss coefficient have been introduced in Equation (1).


### 5.2. The Damage Energy Dissipation Rate

The concept of the damage energy dissipation rate is proposed to establish the mathematical relationship between the damage variable and dissipated energy. In unit volume rock, the damage energy dissipation rate is defined as the energy dissipated to produce unit damage. According to the definition, the expression of the damage energy dissipation rate can be obtained by combining Equation (5), as follows:(9)Y=δadD,where Y is the damage energy dissipation rate (expressed in MJ m-3).To obtain the damage energy dissipation rate, the relationship between the damage dissipated energy and damage variable at unloading points in cyclic loading tests is analyzed. According to the test data, the statistical results are drawn into a scatter diagram, as shown in Figure12. In the figure, the abscissa represents the damage variable, and the ordinate represents the damage dissipated energy. It can be seen from the test results that the relationship between the two is approximately linear. The linear model is used for fitting, and the correlation coefficients under different confining pressures are all greater than 0.9. Therefore, it is considered that the relationship between the damage dissipated energy and the damage variable is linear. According to the definition, the slopes of the fitting lines are the values of the damage energy dissipation rate. From this, we can obtain parameter values at 5 MPa–30 MPa, which are 0.86, 1.34, 2.18, and 2.46, respectively.Figure 12
Relationship between damage dissipated energy and damage variable in cycle loading tests.


### 5.3. Establishment of Damage Evolution Model for Sandstone under Triaxial Loading

Because the tests designed in this paper are conventional triaxial loading tests of sandstone, the damage evolution model is suitable for conventional triaxial loading condition. Therefore, suppose first(10)σ2=σ3.The constitutive relationship of sandstone under triaxial loading can be expressed as(11)σ=1−αDM0εe,where σ is the stress tensor, εe is the elastic strain tensor, and M0 is the initial stiffness matrix of rock in the undamaged state.The above formula can be rewritten as(12)σ1σ2σ3=E1−μ1−αD1+μ1−2μ1μ1−μμ1−μ1μ1−μsym1ε1eε2eε3e,where εiei=1,2,3 represents elastic strain in three directions, respectively.Under triaxial loading conditions, the bearing capacity in the direction of maximum principal stress should take into account the effect of residual strength. Therefore, referring to previous studies [5, 29], the expression of residual strength is as follows:(13)R=2crcosφr+σ31+sinφr1−sinφr,where R is the residual strength and cr and φr are the cohesion and the internal friction angle of rock, respectively, which can be obtained by test data.Combining Equations (12) and (13), the constitutive relationship of rock can be expressed as(14)σ1=ηε1e+2βε3e+DRσ2=σ3=βε1e+η+βε3e,where η=E1−μ1−αD/1+μ1−2μ and β=Eμ1−αD/1+μ1−2μ.Based on energy balance equation (Equation (5)), taking unit volume rock as the object of study, three parts of energy calculation are involved in the process of infinitesimal work. The three parts of energy are expressed in incremental form:(1)
Energy inputted by external workUnder triaxial loading, the differential form of external work can be expressed as follows:(15)δw=σdε,where dε is the increment of strain tensor(2)
Elastic energyThe elastic energy stored inside rock can be expressed as(16)E~=12εeMεe,where M is the stiffness matrix after damage, M=1−αDM0.By calculating the differential of Equation (16), we yield(17)de=σdεe−12αεeM0εedD(3)
Damage dissipated energyDamage dissipated energy can be expressed as the product of damage energy dissipation rate and damage variable, as follows:(18)δa=YdD.According to the law of energy conservation, substituting Equations (15), (17), and (18) into Equation (5) yields(19)dD=2σdεp2Y−αεeM0εe,where εp is the plastic strain tensor and εipi=1,2,3 represents the plastic strain in three directions, respectively.


### 5.4. Solution of Damage Evolution Model

According to the above derivation, Equation (19) is the damage evolution equation. Taking the axial strain increment dε1 as the independent variable, it is necessary to solve the axial elastic strain increment dε1e, the axial plastic strain increment dε1p, the circumferential elastic strain increment dε3e, the circumferential plastic strain increment dε3p, and the damage variable dD. Therefore, we solve the following differential equations simultaneously:(1)
The axial strain increment is the sum of the axial elastic strain increment and the axial plastic strain increment, which can be expressed as(20)dε1=dε1e+dε1p(2)
The increment of axial plastic strain can be expressed as(21)dε1p=Ddε1(3)
The increment of circumferential plastic strain can be expressed as(22)dε3p=Ddε3(4)
The increment of circumferential stress is 0,dσ3=0. According to Equation (14) yields(23)η+βdε3e+βdε1e−ηε3e+βε1e+ε3eαdD=0.The above equations are combined with the damage evolution equation Equation (19), which can be expressed by a matrix as follows:(24)110001−D−D0001−D−D0000β0η+β−αηε3e+βε1e+ε3e−σ10−2σ30Y−12αεeM0εedε1pdε1edε3pdε3edD=dε10000.By solving Equation (24), the damage evolution and stress-strain curves of sandstone under triaxial loading can be obtained.


## 5.1. Definition of Damage Variable

In damage mechanics, rock can be regarded as a continuous medium. The evolution process of the mechanical properties of rock under various stress conditions can be studied by introducing damage variables. In short, the damage variable can be understood as a measure of the degree of decay of the mechanical properties. When the material is damaged, there will be a series of changes. For rock, there will be a decrease in the elastic modulus, generation of plastic strain, increase in volume, decrease in density, generation of acoustic emission and electromagnetic radiation, change in temperature, and so on. Therefore, the damage variable can be defined by the changes in these physical quantities. The damage variable is usually defined as the decay of the elastic modulus for rock in domestic and foreign studies, which can be expressed as(6)D=1−EE0,where D is the damage variable of rock.According to the definition of damage, elastic modulus of damaged rock can be expressed as follows:(7)E=1−DE0.However, Equation (7) is more suitable for describing the elastic damage behavior. To describe the rock damage process with both elasticity and plasticity, the expression needs to be modified. For example, the evolution process of the elastic modulus in the whole stress-strain curve has been analyzed in detail. It can be seen from Figure 7 that in the residual stage, there is still a residual elastic modulus. When the value of the damage variable is 1, the value of the elastic modulus calculated according to Equation (7) is 0, which does not conform to the test results. Therefore, the expression of the elastic modulus needs to be modified:(8)E=1−αDE0,where α is the elastic modulus loss coefficient. The physical meaning and value selection of loss coefficient have been introduced in Equation (1).


## 5.2. The Damage Energy Dissipation Rate

The concept of the damage energy dissipation rate is proposed to establish the mathematical relationship between the damage variable and dissipated energy. In unit volume rock, the damage energy dissipation rate is defined as the energy dissipated to produce unit damage. According to the definition, the expression of the damage energy dissipation rate can be obtained by combining Equation (5), as follows:(9)Y=δadD,where Y is the damage energy dissipation rate (expressed in MJ m-3).To obtain the damage energy dissipation rate, the relationship between the damage dissipated energy and damage variable at unloading points in cyclic loading tests is analyzed. According to the test data, the statistical results are drawn into a scatter diagram, as shown in Figure12. In the figure, the abscissa represents the damage variable, and the ordinate represents the damage dissipated energy. It can be seen from the test results that the relationship between the two is approximately linear. The linear model is used for fitting, and the correlation coefficients under different confining pressures are all greater than 0.9. Therefore, it is considered that the relationship between the damage dissipated energy and the damage variable is linear. According to the definition, the slopes of the fitting lines are the values of the damage energy dissipation rate. From this, we can obtain parameter values at 5 MPa–30 MPa, which are 0.86, 1.34, 2.18, and 2.46, respectively.Figure 12
Relationship between damage dissipated energy and damage variable in cycle loading tests.


## 5.3. Establishment of Damage Evolution Model for Sandstone under Triaxial Loading

Because the tests designed in this paper are conventional triaxial loading tests of sandstone, the damage evolution model is suitable for conventional triaxial loading condition. Therefore, suppose first(10)σ2=σ3.The constitutive relationship of sandstone under triaxial loading can be expressed as(11)σ=1−αDM0εe,where σ is the stress tensor, εe is the elastic strain tensor, and M0 is the initial stiffness matrix of rock in the undamaged state.The above formula can be rewritten as(12)σ1σ2σ3=E1−μ1−αD1+μ1−2μ1μ1−μμ1−μ1μ1−μsym1ε1eε2eε3e,where εiei=1,2,3 represents elastic strain in three directions, respectively.Under triaxial loading conditions, the bearing capacity in the direction of maximum principal stress should take into account the effect of residual strength. Therefore, referring to previous studies [5, 29], the expression of residual strength is as follows:(13)R=2crcosφr+σ31+sinφr1−sinφr,where R is the residual strength and cr and φr are the cohesion and the internal friction angle of rock, respectively, which can be obtained by test data.Combining Equations (12) and (13), the constitutive relationship of rock can be expressed as(14)σ1=ηε1e+2βε3e+DRσ2=σ3=βε1e+η+βε3e,where η=E1−μ1−αD/1+μ1−2μ and β=Eμ1−αD/1+μ1−2μ.Based on energy balance equation (Equation (5)), taking unit volume rock as the object of study, three parts of energy calculation are involved in the process of infinitesimal work. The three parts of energy are expressed in incremental form:(1)
Energy inputted by external workUnder triaxial loading, the differential form of external work can be expressed as follows:(15)δw=σdε,where dε is the increment of strain tensor(2)
Elastic energyThe elastic energy stored inside rock can be expressed as(16)E~=12εeMεe,where M is the stiffness matrix after damage, M=1−αDM0.By calculating the differential of Equation (16), we yield(17)de=σdεe−12αεeM0εedD(3)
Damage dissipated energyDamage dissipated energy can be expressed as the product of damage energy dissipation rate and damage variable, as follows:(18)δa=YdD.According to the law of energy conservation, substituting Equations (15), (17), and (18) into Equation (5) yields(19)dD=2σdεp2Y−αεeM0εe,where εp is the plastic strain tensor and εipi=1,2,3 represents the plastic strain in three directions, respectively.


## 5.4. Solution of Damage Evolution Model

According to the above derivation, Equation (19) is the damage evolution equation. Taking the axial strain increment dε1 as the independent variable, it is necessary to solve the axial elastic strain increment dε1e, the axial plastic strain increment dε1p, the circumferential elastic strain increment dε3e, the circumferential plastic strain increment dε3p, and the damage variable dD. Therefore, we solve the following differential equations simultaneously:(1)
The axial strain increment is the sum of the axial elastic strain increment and the axial plastic strain increment, which can be expressed as(20)dε1=dε1e+dε1p(2)
The increment of axial plastic strain can be expressed as(21)dε1p=Ddε1(3)
The increment of circumferential plastic strain can be expressed as(22)dε3p=Ddε3(4)
The increment of circumferential stress is 0,dσ3=0. According to Equation (14) yields(23)η+βdε3e+βdε1e−ηε3e+βε1e+ε3eαdD=0.The above equations are combined with the damage evolution equation Equation (19), which can be expressed by a matrix as follows:(24)110001−D−D0001−D−D0000β0η+β−αηε3e+βε1e+ε3e−σ10−2σ30Y−12αεeM0εedε1pdε1edε3pdε3edD=dε10000.By solving Equation (24), the damage evolution and stress-strain curves of sandstone under triaxial loading can be obtained.


## 6. Model Result and Discussion

### 6.1. Parameters of Model

The parameters involved in the damage evolution model proposed in this paper can basically be obtained through experiments. The initial elastic modulusE0 and Poisson’s ratio μ can be determined experimentally, and the values are shown in Table 1. The elastic modulus loss coefficient α can be obtained by calculating the ratio of the residual elastic modulus to the initial elastic modulus, as shown in Table 2. The nonlinear coefficient k is a material parameter that will not change with loading conditions. The nonlinear coefficient is 0.05 by the least-squares method. The damage energy dissipation rate Y can also be determined through experiments, and the specific method is described in Section 5.2.


### 6.2. Results and Discussion

The damage curves and axial stress-strain curves predicted by the damage evolution model under different confining pressure conditions are shown in Figures13–16 and compared with the test results and the previous models proposed by Cao et al. and Zhu et al. [10, 11, 48].Figure 13
Axial stress-strain curve and damage curve under confining pressure of 5 MPa.Figure 14
Axial stress-strain curve and damage curve under confining pressure of 10 MPa.Figure 15
Axial stress-strain curve and damage curve under confining pressure of 20 MPa.Figure 16
Axial stress-strain curve and damage curve under confining pressure of 30 MPa.The damage evolution model proposed by [48] is as follows:(25)Rd=4rcd/dc1+d/dc2,ε1=1E′σ1−2ν0E′σ3−26−η′3d2Rdχε2=1−ν0E′σ3−ν0E′σ1+16+η′3d2Rdχε3=1−ν0E′σ3−ν0E′σ1+16+η′3d2Rdχ.The model parameters are shown in Table4.Table 4
The parameters of [48].
ν0E′, GPaη′rc, J m-2dc0.25251.50.082The damage evolution model proposed by [10, 11] is as follows:(26)D=1−exp−F/F0mF≥00F<0,εi=γ01−exp−σi−μ1σj+σkE1+1−γ0σi−μ2σj+σkE2exp−F/F0mF≥0γ01−exp−σi−μ1σj+σkE1+1−γ0σi−μ2σj+σkE2F<0.The model parameters are shown in Table5.Table 5
The parameters of [10, 11].
γ0E1, GPaE2, GPaμ1μ20.000414.9250.40.25The nonlinear characteristics of the initial compaction stage are not considered in the previous model, so the mechanical behaviors of rock in the compaction stage cannot be well described by Zhu et al. The previous model established by Cao could represent nonlinear characteristics, but the residual strength is not considered in the model. It can be seen from the comparison results that the damage evolution model proposed in this paper can reflect the whole axial loading process and describe the peak strength and residual strength of rock. The model shows that with increasing confining pressure, the rock strength increases, and the mechanical properties of rock change from brittleness to ductility. It can also be found from the comparison results that the theoretical model can reflect the nonlinear properties of rock, and the nonlinear behaviors gradually change to linear with increasing confining pressure.The relationships between the damage variable and axial strain are shown in Figures13–16. The relationships between the damage variable and axial stress are shown in Figure 17. The damage variable has different evolution laws at different stress stages under constant confining pressure. In the initial stage, the evolution rate of damage is relatively slow. With increasing axial strain, the evolution rate of damage increases. As shown in Figure 17, damage evolution enters the rapid evolution stage from the yield stage. In the residual stage, the damage evolution rate decreases as the damage variable approaches 1. The above theoretical calculation results are consistent with objective physical laws, and the rationality of the model is verified.Figure 17
Damage-axial stress curves under different confining pressures.It can also be seen from the comparison results that the model errors mainly occur in the postpeak softening stage. Under different confining pressures, there will be different degrees of errors. From the comparison results, when the confining pressure condition is larger, the model errors are relatively small. Damage evolution is the main factor affecting the stress-strain relationship in the postpeak softening stage. In the softening stage, when the damage evolution rate increases rapidly, the stress curve decreases rapidly. When the damage evolution rate is relatively gentle, the stress curve decreases slowly.The damage evolution model is used to simulate the postpeak unloading process under different confining pressures. The comparison between the model results and test results is shown in Figure18. In the figure, the red curves represent the model results, and scattered points represent the test results. It can be seen from the comparison results that the model results can reflect the postpeak unloading behaviors of rock to a certain extent and can describe nonlinear properties during the postpeak unloading process. The nonlinear characteristics during the unloading process will also be affected by confining pressures. The smaller the confining pressure is, the more pronounced the nonlinear relationship between stress and strain.Figure 18
Comparison of test results and model results of postpeak unloading under different confining pressures ((a) 5 MPa, (b) 10 MPa, (c) 20 MPa, and (d) 30 MPa).
(a)(b)(c)(d)It can be seen from Figure18 that there are some errors in the model curves, and the more obvious errors in the figure are marked out. Through analysis, we find that the errors of the unloading model mainly exist in two stages. At the beginning of unloading, as the load is gradually unloaded, the strain does not immediately recover but first evolves very slowly and then rebounds gradually. The hysteresis of strain cannot be reflected by the damage evolution model. On the other hand, at the stage after complete unloading, there are also significant errors in the model curves, and the strain values calculated by the model are often smaller than the test results. Specifically, after unloading the load, the theoretically calculated plastic strain is slightly smaller than the test result. During the test, after the load is removed, it is possible that the rock strain does not have enough time to recover before the next loading and unloading. If the specimen can be allowed to stand for a long enough time after each unloading, the strain should rebound further. In addition, the damage evolution model in this paper can only simulate the postpeak unloading process and cannot reflect the stress-strain relationship during prepeak unloading. In the prepeak unloading process, the change in elastic modulus is more complex, i.e., increasing first and then decreasing. We will study this further in future work.


## 6.1. Parameters of Model

The parameters involved in the damage evolution model proposed in this paper can basically be obtained through experiments. The initial elastic modulusE0 and Poisson’s ratio μ can be determined experimentally, and the values are shown in Table 1. The elastic modulus loss coefficient α can be obtained by calculating the ratio of the residual elastic modulus to the initial elastic modulus, as shown in Table 2. The nonlinear coefficient k is a material parameter that will not change with loading conditions. The nonlinear coefficient is 0.05 by the least-squares method. The damage energy dissipation rate Y can also be determined through experiments, and the specific method is described in Section 5.2.


## 6.2. Results and Discussion

The damage curves and axial stress-strain curves predicted by the damage evolution model under different confining pressure conditions are shown in Figures13–16 and compared with the test results and the previous models proposed by Cao et al. and Zhu et al. [10, 11, 48].Figure 13
Axial stress-strain curve and damage curve under confining pressure of 5 MPa.Figure 14
Axial stress-strain curve and damage curve under confining pressure of 10 MPa.Figure 15
Axial stress-strain curve and damage curve under confining pressure of 20 MPa.Figure 16
Axial stress-strain curve and damage curve under confining pressure of 30 MPa.The damage evolution model proposed by [48] is as follows:(25)Rd=4rcd/dc1+d/dc2,ε1=1E′σ1−2ν0E′σ3−26−η′3d2Rdχε2=1−ν0E′σ3−ν0E′σ1+16+η′3d2Rdχε3=1−ν0E′σ3−ν0E′σ1+16+η′3d2Rdχ.The model parameters are shown in Table4.Table 4
The parameters of [48].
ν0E′, GPaη′rc, J m-2dc0.25251.50.082The damage evolution model proposed by [10, 11] is as follows:(26)D=1−exp−F/F0mF≥00F<0,εi=γ01−exp−σi−μ1σj+σkE1+1−γ0σi−μ2σj+σkE2exp−F/F0mF≥0γ01−exp−σi−μ1σj+σkE1+1−γ0σi−μ2σj+σkE2F<0.The model parameters are shown in Table5.Table 5
The parameters of [10, 11].
γ0E1, GPaE2, GPaμ1μ20.000414.9250.40.25The nonlinear characteristics of the initial compaction stage are not considered in the previous model, so the mechanical behaviors of rock in the compaction stage cannot be well described by Zhu et al. The previous model established by Cao could represent nonlinear characteristics, but the residual strength is not considered in the model. It can be seen from the comparison results that the damage evolution model proposed in this paper can reflect the whole axial loading process and describe the peak strength and residual strength of rock. The model shows that with increasing confining pressure, the rock strength increases, and the mechanical properties of rock change from brittleness to ductility. It can also be found from the comparison results that the theoretical model can reflect the nonlinear properties of rock, and the nonlinear behaviors gradually change to linear with increasing confining pressure.The relationships between the damage variable and axial strain are shown in Figures13–16. The relationships between the damage variable and axial stress are shown in Figure 17. The damage variable has different evolution laws at different stress stages under constant confining pressure. In the initial stage, the evolution rate of damage is relatively slow. With increasing axial strain, the evolution rate of damage increases. As shown in Figure 17, damage evolution enters the rapid evolution stage from the yield stage. In the residual stage, the damage evolution rate decreases as the damage variable approaches 1. The above theoretical calculation results are consistent with objective physical laws, and the rationality of the model is verified.Figure 17
Damage-axial stress curves under different confining pressures.It can also be seen from the comparison results that the model errors mainly occur in the postpeak softening stage. Under different confining pressures, there will be different degrees of errors. From the comparison results, when the confining pressure condition is larger, the model errors are relatively small. Damage evolution is the main factor affecting the stress-strain relationship in the postpeak softening stage. In the softening stage, when the damage evolution rate increases rapidly, the stress curve decreases rapidly. When the damage evolution rate is relatively gentle, the stress curve decreases slowly.The damage evolution model is used to simulate the postpeak unloading process under different confining pressures. The comparison between the model results and test results is shown in Figure18. In the figure, the red curves represent the model results, and scattered points represent the test results. It can be seen from the comparison results that the model results can reflect the postpeak unloading behaviors of rock to a certain extent and can describe nonlinear properties during the postpeak unloading process. The nonlinear characteristics during the unloading process will also be affected by confining pressures. The smaller the confining pressure is, the more pronounced the nonlinear relationship between stress and strain.Figure 18
Comparison of test results and model results of postpeak unloading under different confining pressures ((a) 5 MPa, (b) 10 MPa, (c) 20 MPa, and (d) 30 MPa).
(a)(b)(c)(d)It can be seen from Figure18 that there are some errors in the model curves, and the more obvious errors in the figure are marked out. Through analysis, we find that the errors of the unloading model mainly exist in two stages. At the beginning of unloading, as the load is gradually unloaded, the strain does not immediately recover but first evolves very slowly and then rebounds gradually. The hysteresis of strain cannot be reflected by the damage evolution model. On the other hand, at the stage after complete unloading, there are also significant errors in the model curves, and the strain values calculated by the model are often smaller than the test results. Specifically, after unloading the load, the theoretically calculated plastic strain is slightly smaller than the test result. During the test, after the load is removed, it is possible that the rock strain does not have enough time to recover before the next loading and unloading. If the specimen can be allowed to stand for a long enough time after each unloading, the strain should rebound further. In addition, the damage evolution model in this paper can only simulate the postpeak unloading process and cannot reflect the stress-strain relationship during prepeak unloading. In the prepeak unloading process, the change in elastic modulus is more complex, i.e., increasing first and then decreasing. We will study this further in future work.


## 7. Conclusion

Based on cyclic loading and unloading tests of sandstone samples, the evolution process of the elastic modulus and the energy evolution of the rock system are analyzed. A damage evolution model of sandstone under triaxial loading is established. Compared with the test results, the following conclusions can be obtained:(1)
An elastic modulus model with the average stress as the influencing factor is established to describe the nonlinear behaviors during loading(2)
According to the law of energy conservation, an energy balance equation is established, which can be simply expressed as the energy inputted by external work being transformed into elastic energy and damage dissipated energy(3)
The comparison results show that the model can reflect the axial stress-strain relationship and damage the evolution process and nonlinear behaviors of rock(4)
The nonlinear behaviors and stress-strain curve in the postpeak unloading process can be described by the model to a certain extent.



---

*Source: 2395789-2021-11-27.xml*