A New Damage Model of Rock Considering Residual Strength under Seepage Pressure

Seepage pressure changes the rock structure and stress state, inducing rock mass damage and destruction under hydromechanical coupling. To describe the deformation characteristic and the entire stress–strain process, we suppose that rock consists of damaged part, undamaged part, and voids part. The three parts are intermingled to carry external loads and the undamaged part gradually becomes damaged part in compression. A redefinition of the damage variable and evolution equation is proposed incorporating damage mechanics and statistical damage theory. Considering that axial loads of the damage parts equal to residual strength, a new statistical damage softening constitutive model is developed. The model is validated by using the experimental data. The proposed model can reflect the entire stress–strain law very well, especially softening deformation characteristics after the peak. The research results are helpful in the rock mass engineering safety analysis under seepage pressure.


Introduction
The seepage and control of water have applications in hydropower engineering, mining engineering, underground engineering, underground nuclear waste geological disposal, geothermal resource development, and other emerging rock engineering areas. These are related to the coupling effect of rock mass and fluid interaction, mutual connection, and mutual restriction. The analysis of the hydromechanical behavior of rock masses has become an essential topic in rock mechanics.
Extensive hydraulic coupling tests of rock have revealed the mechanical behavior of rock under hydromechanical coupling [1][2][3][4][5][6]. The mechanical behavior of rock is related to confining pressure and has an apparent relationship with seepage pressure. Higher confining pressure can improve rock's peak and residual strength and constrain rock deformation, while seepage pressure can reduce confining pressure so that the strength decreases and the deformation increases. Higher seepage pressure can cause rock damage and create microcracks. Microcracks induced by hydraulic pressure will drastically increase the rock's permeability and cause further damage and weakening [7][8][9]. Therefore, rock mechanics under seepage pressure are complex, and their deformation behavior can only be approximated using statistical approaches.
The two main statistical approaches used in rock mechanics are the micromechanical and phenomenological approaches. The micromechanical approach accounts for the size, shape, interaction mode, and spatial distribution of cracks [10,11]. This method is very suitable for explaining the mechanical characteristics of rock from the perspective of crack structure. However, it is difficult to apply directly to solving practical rock engineering problems because of the complexities of statistical analysis. The phenomenological approach reveals the mechanical behavior of rock using effective stress principles for saturated geomaterial. This approach is feasible for solving practical rock engineering problems, although it fails to consider relevant microstructural features [1,2,12,13].
The statistical damage theory builds effective constitutive models of damaged rock based on the phenomenological approach combining continuum damage mechanics and statistical theory. The statistical damage theory [14] was first proposed to describe randomly distributed internal defects of rock materials. The widespread application of this theory has contributed to its further development. Research [15][16][17][18][19][20] has shown that the SDCM (statistical damage constitutive model) is a continuous damage evolution equation represented by the equivalent stress corresponding to a yield criterion. However, SDCM mainly focused on the stressstrain process under external loads, such as axial stress and circumferential stress. The effect of seepage pressure on rock mechanical behavior is not considered. The new coupling permeability statistical damage model developed by Gao et al. [21] does not adequately reflect the residual strength behavior.
The paper focuses on the entire stress-strain relationship and residual strength behavior for hydromechanical coupling. Assuming that rock consists of damaged area, undamaged area, and voids area, a new SDCM of rock considering residual strength behavior under seepage pressure is proposed. The model parameters were determined from experimental results. The model reasonably describes the deformation behavior of rocks under seepage pressure, and the sensitivity to important parameters was analyzed.

Characteristics of the Entire Stress-Strain Curve for Hydromechanical Coupling
According to triaxial compression test results under hydromechanical coupling in the laboratory, the rock deformation process may be divided into the following five stages: compaction, linear-elastic, plastic, strain-softening, and residual strength. Figure 1 shows a schematic diagram of the entire stress-strain curve of rock in the laboratory. Seepage and confining pressure both have important effects on rock strength and deformation. Rock strength decreases with increasing seepage pressure but increases with increasing confining pressure. The axial deformation of the compression and linear-elastic stages increases significantly [3,4]. Compared with axial deformation, lateral deformation is more sensitive to confining and seepage pressure [12].
Although the rock deformation process theoretically has the above five stages, actual stress-strain curve measurements do not always include the last two stages, strainsoftening and residual strength. On the one hand, in the existing standards and criteria of rock mechanics, peak strength is used as the primary mechanical index because it is not necessary to characterize the postpeak strength and deformation behavior. On the other hand, the strain range of residual strength is about 5-10 times larger than the strain at peak strength, so most conventional compression tests often stop loading before reaching the residual deformation stage [22,23]. However, this does not mean that the last two stages do not exist or that the rock has no residual strength. Therefore, it is necessary to analysis deformation characteristics of all five stages for hydromechanical coupling applications.

SDCM for Hydromechanical Coupling
3.1. The Extended Damage Model for Rock. According to Lemaitre's strain equivalence principle, the damage constitutive relationship of rocks is represented as where σ i denotes the apparent stress acting on the rock, σ u i denotes the effective stress acting on the rock, and D is the damage variable, which varies from the intact state (D = 1) to the fully damaged state (D = 0) of rock. The damage model supposes that the essential cause of damage is the creation and growth of material defects, such as microvoids or microcracks. Once these defects form, the representative volume element will not have any bearing capacity. However, this is inconsistent with the actual situation in rock deformation process. When the rock is fully damaged (D = 1), the apparent axial stress is zero, and the residual strength will also be zero. From Section 2, we know that the rock deformation process has five stages, including the strain-softening and residual strength stages. Some new damage models have been proposed [8,18,21] to describe the complete rock deformation stress-strain curve. These new models assume that the rock is composed of two areas or three areas, as shown in Figures 2(a) and 2(b).
Rock is saturated in indoor triaxial compression tests of hydromechanical coupling, filling the rock's pores and fissures with water. The pressure generated by the water acts on the pores and fissures, causing damage to the rock's internal structure. Therefore, seepage pressure cannot be ignored. Hypothesizing that rock elements consist of three areas-an undamaged area, a damaged area, and voids-the void resistance to loading equals the seepage water pressure (p w ). This provides a feasible model for analyzing deformation under seepage pressure. Therefore, the extended damage model with three areas can represent the stress-strain relationship for hydromechanical coupling.

Basic
Assumptions for the SDCM. The extended damage model is adopted to describe the complete deformation process of rock under seepage pressure. Considering that 2 Geofluids rock is a heterogeneous natural geological body with internal defects, some basic assumptions are necessary to develop the damage model. The major assumptions are as follows: (1) The extended conceptual model for rock consists of three areas: damaged area, undamaged area, and voids, as shown in Figure 2(b). The three areas support the external loads together (2) The damage is isotropic. The undamaged area is elastic, and the generalized Hooke's law can be used (3) When bearing load of mesoscopic elements is larger than their strength, the available mesoscopic elements gradually transform from undamaged into damaged areas, with the total number of elements invariant over increasing external loads. The effective stress on the damaged area can be approximately regarded as the residual strength of the rock in axial direction [8] (4) To represent the heterogeneity, the strength of mesoscopic elements (F) satisfies the Weibull distribution. The probability density function PðFÞ can be written as where F represents the strength of any mesoscopic element, and m is the shape parameter and F 0 is the scale parameter in Weibull distribution's statistical parameters 3.3. Constitutive Stress Equilibrium. According to assumption (1), rock consists of three areas: damaged area, undamaged area, and voids, as shown in Figure 2(b). The quantities A u , A d , and A 0 represent the undamaged, damaged, and voids (pore and fissure) cross-sectional areas, respectively, and A is total cross-section area. The total stress of the cross-section is σ i , consisting of three areas: σ u i , σ d i , and p w are the undamaged, damaged, and void stresses, respectively. As mentioned previously, the voids resist a possible loading equivalent to the seepage water pressure. After employing the static equilibrium condition on the cross-section, the following relations hold: Porosity n is the area fraction of the voids, defined by, and the damage variable (D) is defined as By substituting Equations (4), (5), and (6) into Equation (3), the following relationship is obtained: According to assumption (2), the undamaged area is elastic, and its behavior follows Hooke's law, so σ u i is written as follows: where ði, j, kÞ = ð1, 2, 3Þ, and E u , v u , and ε u i are the effective Young's modulus, effective Poisson's ratio, and effective elastic strain of rock undamaged area, respectively.
Based on the relationship of the material properties between apparent and effective parameters [8]), where E, μ, and ε i are Young's modulus, Poisson's ratio, and macroscale specimen strain, respectively. Therefore, the following equation can be obtained: and the relationships for the conventional triaxial compression test in laboratory can be obtained: Substituting Equation (11) into Equation (10) yields and substituting Equation (12) into Equation (7), the following equation is acquired for the axial direction: According to laboratory rock hydromechanical coupling tests, the seepage pressure is approximately linearly distributed over the sample height, with the maximum seepage pressure at the top of the sample. Because of its contact with the atmosphere, the seepage pressure is zero (relative to atmospheric pressure) at the bottom of the sample. Therefore, the effective confining pressure σ u 3 can be calculated as The measured axial deviatoric stress σ 1t is actually the difference between the effective axial stress and the effective confining pressure in the hydromechanical coupling test. The measured axial deviatoric stress can be expressed as After the axial deviatoric stress is applied, the measured strain ε 1t is equal to the axial strain ε 1 : Substituting Equations (14)- (16) into Equation (13), the measured axial deviatoric stress expression is From assumption (3), once an undamaged area becomes damaged, the axial stress value equals the residual strength. Therefore, Equation (17) can be rewritten as where R is the residual strength of the rock.

Damage Evolution Equation.
According to assumptions (3) and (4), the damage variable D can be defined as the number ratio of failure elements N f to the total elements N. The strength's heterogeneity satisfies the Weibull density distribution function. When the mesoscopic element strength F is zero or greater, the mesoscopic element is failure, and the element becomes damaged state. As the strength changes from F to F + dF, the number of damaged elements can be determined as PðFÞNdF, and the following relationship is obtained: When F is less than zero, the mesoscopic element is undamaged, and D = 0. Therefore, The failure criterion for rock mesoscopic elements can be written in a uniform expression The Mohr-Coulomb criterion is commonly used and can be expressed as Substituting Equations (12), (14), (15), and (16) into Equation (22), this relationship can be expressed as In summary, the SDCM for hydromechanical coupling can be expressed as 3.5. Determination of Damage Model Parameters. The proposed SDCM for hydromechanical coupling, Equation (24), has two types of parameters. One type consists of parameters related to the characteristics of the rock itself, such as porosity n, residual strength R, Young's modulus E, Poisson's ratio μ, friction angle φ, and cohesion c. The characteristics of specific rock types determine these parameters by the triaxial compression tests. The other type consists of statistical parameters, such as m and F 0 , which are discussed below. The porosity n continuously changes during the loading process. It is difficult to directly determine its value, so the relationship between porosity and apparent parameters is expressed indirectly [15]. We adopt the following relationship: where n 0 is the initial porosity. The residual strength R is related to confining pressure σ 3 [24]. It can be calculated based on the Mohr-Coulomb strength criterion at the residual strength stage, as follows: where c r is the residual cohesion and φ r is the residual friction angle. The statistical parameters m and F 0 can be estimated using the extremum method [25]. At the peak point of the stress-strain curve, the slope is zero. That is, the derivative of stress with respect to axial strain is equal to zero, The maximum axial principal stress σ 1p can be expressed in terms of the corresponding axial strain at the peak point, Therefore, m and F 0 can be expressed as where A and B are calculated parameters defined by the following equations:

Validation of the SDCM
Hydromechanical coupling experimental data for fractured limestone from Zhao et al. [24] and Zhao et al. [20] are adopted to validate the new SDCM. The limestone was taken from a construction site in south China. The P-wave velocity of the limestone is in the range of 1.8 to 2.2 km/s. Initial porosity is in the range of 12%-16%, so the initial porosity is taken as 14% for convenience. The rock mechanical parameters of the fractured limestone are shown in Table 1. Friction angle and cohesion, residual friction angle, and residual cohesion can be calculated from the peak strength, residual strength, and confining pressure. These results are shown in Table 2.
The parameters in Equation (24) must be determined first. The porosity n is calculated from Equation (26), and the residual strength R is calculated from Equation (27). The statistical parameters m and F 0 are calculated using Equations (28) and (29) from the extremum method. Therefore, the relationship between axial stress and axial strain can be calculated, and the results are shown in Figure 3. The computed model parameters are shown in Table 3.
In Figure 3, the blue circles represent the experimental curves from the literature [8], and the red squares represent the calculated curves from the proposed model. These results show that the theoretical curve is a good fit for the experimental results. The model developed in this paper reflects the behavior of the complete stress-strain curve of rock under the action of seepage pressure, especially for the postpeak curve. From the results at three seepage pressures, it is also found that the model can represent the linear behavior before the peak but cannot effectively reflect the nonlinear deformation stage at the beginning of loading. This is mainly because there is no damage and only elastic deformation following Hooke's law when F < 0. Because Figure 2: Schematic diagram of the extended damage models for rock: (a) rock is composed of undamaged area and damaged area [18] and (b) rock is composed of undamaged area, damaged area, and voids [21].    6 Geofluids the calculation is carried out using the peak point method, it agrees well with the test curve at the peak point position.
What is more, it can approximately reflect the strainsoftening characteristics in the postpeak stage.
To further validate the proposed model, the calculated results using the effective stress principle, the method had been adopted by Wang et al. [1] and Wang Wei [2] are also plotted in Figure 3. This comparison shows that the new proposed SDCM is superior to the model based on the effective stress principle, and the theoretical curve in this paper is closer to the experimental curve, especially in the postpeak stage.

Parametric Sensitivity Analysis and Discussion
The rock deformation depends on several parameters, as shown in Equation (24), among which the deformation modulus and initial porosity are the most important. A parametric analysis has been carried out to understand the effects of these two parameters better.
The deformation modulus is a comprehensive value related to confining pressure and seepage pressure. According to Zhao et al. [8], deformation modulus continued to rise with increasing confining pressure and reduced gently with increasing differential seepage pressure. The deformation of fractured limestone is typically quantified by calculating the deformation modulus E 50 at the mid of the peak strength. However, the deformation modulus in this paper attempts to reflect the entire stress-strain process, so it is different from that of Zhao et al. [8] when the inversion calculation is carried out.
An analysis of the sensitivity to deformation modulus is conducted by varying its value from 25 GPa to 37 GPa in   3 GPa increments while the other parameters remain constant with a 2 MPa seepage pressure and 7 MPa confining pressure. These results are shown in Figure 4. It is found that different deformation modulus values have a strong influence on the shape of deformation curve. As the deformation modulus increases, the slope of prepeak stage increases. The stress decreases more slowly after the peak, and the strainsoftening gradually transforms to strain-hardening.
The sensitivity to initial porosity n 0 is studied by varying its value from 0.05 to 0.25 in increments of 0.05 while the other parameters remain constant with a 2 MPa seepage pressure and a 7 MPa confining pressure. From Figure 5, the stress-strain curve tends to become convex as the axial strain increases. The rock stiffness is smaller as the initial porosity increases, implying that the voids and pores reduce the strength during continuous compaction. However, the strength decreases quickly after the peak. The larger the initial porosity, the greater the strength attenuation.

Conclusions
The contributions and results of this study can be summarized as follows: (1) A statistical damage model for hydromechanical coupling has been developed using a new damage parameter definition, and the feasibility of the model has been validated against experimental results for fractured limestone. The model can simulate the entire rock deformation process and better reflects the postpeak residual deformation under hydromechanical coupling (2) The deformation modulus E and initial porosity n 0 strongly influence rock deformation. With increasing deformation modulus, the prepeak stress-strain slope becomes larger, the stress decreases more slowly after the peak, and the strain-softening gradu-ally becomes strain-hardening. In addition, it was found that larger initial porosity results in greater strength attenuation

Data Availability
No data were used to support this study.

Conflicts of Interest
The author declares that they have no conflicts of interest.