A Chemical Damage Creep Model of Rock Considering the Influence of Triaxial Stress

In order to accurately describe the characteristics of each stage of rock creep behavior under the combined action of acid environment and true triaxial stress, based on damage mechanics, chemical damage is connected with elastic modulus; thus, the damage relations considering creep stress damage and chemical damage are obtained. The elastic body, nonlinear Kelvin body, linear Kelvin body, and viscoelastic–plastic body (Mogi–Coulomb) are connected in series, and the actual situation under the action of true triaxial stress is considered at the same time. Therefore, a damage creep constitutive model considering the coupling of rock acid corrosion and true triaxial stress is established. The parameters of the deduced model are identified and verified with the existing experimental research results. The yield surface equation of rock under true triaxial stress is obtained by data fitting, and the influence of intermediate principal stress on the creep model is discussed. The derived constitutive model can accurately describe the characteristics of each stage of true triaxial creep behavior of rock under acid environment.


Introduction
Rock rheology is the study of how rocks' mechanical properties change with time, including rheological properties such as creep, stress lax, elastic aftereffect, etc., and is closely related to the long-term stability and safety of rock engineering. A large number of engineering practices and theoretical studies have shown [1][2][3] that the damage and instability found in rock engineering usually occur during the project operation, and the damage to property and loss of life caused by rock rheology are often very serious [4,5]. Therefore, studying the rheological constitutive model of rock has important theoretical significance for the stability and safety of practical engineering. In recent years, many staged results have been achieved in the derivation and establishment of rock creep constitutive models: researchers established various creep models including the empirical model, component combination model, and phenomenological model [6][7][8][9][10]. Hou et al. [11] obtained a four-element combination model considering the damage of sandstone which was well-verified according to the existing experimental data. Jiang et al. [12] obtained a Nishihara model and the correctness of the model was verified by the experimental data of sandstone, which reflected the characteristics of multiple stages of the creep process. Zhou et al. [13] improved the Nishihara model, considering the damage of rock, and accurately verified the model. Wang et al. [14] regarded the macroscopic parameters of rock as a function related to time, and derived a nonlinear model that could reflect the creep Materials 2022, 15, 7590 2 of 13 characteristics of rock. Liu et al. [15] carried out an indoor triaxial creep test on the deep surrounding rock from Fuxin Hengda Coal Mine, and analyzed the creep deformation law of the rock, and a constitutive model that could reflect the nonlinear accelerated creep stage was obtained and verified with experimental data. Zhang et al. [16] pointed out the inaccuracy of parameter identification in the process of establishing the creep model in previous studies, and established a model that could describe the instantaneous elastic strain, decay creep, isokinetic creep, and accelerated creep stages.
The rock creep models established by the above scholars were all obtained under uniaxial or conventional triaxial conditions, and the research on rock creep under the true triaxial stress state is still relatively lacking. In practical engineering, the rock is often in a state of three-phase stress, and the existence of intermediate principal stress will significantly affect the instantaneous and time-dependent properties of the rock [17,18].
At the same time, the environment in which rock formation occurs was often complex and changeable. For example, in the acid rain area of the Sichuan Basin in southwest China, the acidic environment of the rock will have an impact on the mechanical properties of rock engineering. The chemical composition contained in the acidic solution reacts with the rock, changing its microstructure and reducing its bearing capacity. Many studies have shown [19][20][21][22][23] that chemical solutions have a significant impact on the development of damage and the mechanical properties of rock, but these results are mostly aimed at the transient mechanical properties of rocks, including laboratory tests and models; furthermore, research on the creep characteristics of rocks under chemical corrosion is relatively lacking.
Therefore, in order to study the time-dependent characteristics of rock in acidic environment and true triaxial stress, based on damage mechanics, chemical damage is connected with elastic modulus; thus, the damage relations considering creep stress damage and chemical damage are obtained. The elastic body, nonlinear Kelvin body, linear Kelvin body, and viscoelastic-plastic body (Mogi-Coulomb) are connected in series, and the actual situation under the action of true triaxial stress is considered at the same time. Therefore, a damage creep constitutive model considering the coupling of rock acid corrosion and true triaxial stress is established. The parameters of the deduced model are identified and verified with the existing experimental research results. The yield surface equation of rock under true triaxial stress is obtained by data fitting, and the influence of intermediate principal stress on the creep model is discussed. The derived constitutive model can well describe the characteristics of each stage of true triaxial creep behavior of rock under acid environment.

Damage Analysis
Zhang [24] improved the Lemaitre strain equivalent principle by defining the base damage state of rock. When rock is subject to the action of external force, the damage level will expand. Take two damage states: for rock materials in acidic environment, the effective stress of rock under the acid damage state acts on the strain caused by the acid; secondly, there is the stress joint damage state, which is equivalent to the effective stress of rock under the acid. Therefore, the stress joint damage state acts on the strain under the acid damage state: where: σ 1 and E 1 are the effective stress and elastic modulus of rock in acid damage state, respectively, σ 2 and E 2 are the effective stress and elastic modulus of rock under the condition of acid damage and stress damage, respectively. Therefore, the next derivation process can be carried out: Firstly, according to Lemaitre strain equivalence principle, it is obtained that: where: Dc is the chemical damage caused by chemical corrosion, and E 0 is the initial elastic modulus of rock. Then, according to the improved Lemaitre strain equivalence principle, it is obtained that: where: D m is the stress damage caused by stress loading. Substitute Equation (2) into Equation (3) to obtain: According to the Equation (4), we can obtain: So far, we obtain the expression of D s as follows: where: D s is the total damage caused by chemical corrosion and stress loading.

Chemical Damage Variable
When the engineering rock was in an acidic environment, on the one hand, the rock was subjected to a load to form stress damage. On the other hand, the chemical composition in the acidic environment reacted with the rock, which caused the cement inside the rock to dissolve, so that the pores increased, the micro-cracks extended, and chemical damage was formed. The development and penetration of micro-cracks in the rock reflected the damage of the rock. Macroscopically, it showed the decrease in the elastic modulus. Therefore, according to damage mechanics, the damage degree expression of chemical corrosion on rock is obtained: where: E C is the elastic modulus of rock in the environment at pH = C.

Stress Damage Variable
When the stress on the rock is greater than its long-term strength σs, we believe that the rock will suffer from stress damage. According to the Kachanov creep damage law [25], we can obtain: where A and N are the parameters related to the material, and D m is the stress damage variable.
According to Formula (8) and boundary conditions D m = 1, the evolution equation of the damage variable and creep time of rock during accelerated creep process is obtained as:

Comprehensive Damage Variables under Chemical Corrosion and Stress
The damage caused by chemical corrosion exists throughout the creep process, whereas the stress damage occurs when the stress on the rock is greater than its long-term strength σ s . Combining Equations (6), (7), and (9), the comprehensive damage variables of rocks under the action of acid corrosion and stress can be obtained:

Establishment of True Triaxial Creep Model
In the rock creep test, when the initial stress level is less than the long-term strength of the rock, the rock sample will first generate an instantaneous elastic strain during the loading process, since the loading time at this stage accounts for a very small part of the entire creep time. It can be considered that the elastic strain is completed instantaneously. The constitutive relationship at this stage can be described as an elastic body, and considering the true triaxial stress of the rock in the actual situation, the constitutive equation is: where ε e ij is the elastic strain under true triaxial stress, S ij is the deviatoric stress tensor, δ ij σ m is the spherical stress tensor, G 1 is the shear modulus in the elastic strain stage, and K 1 is the bulk modulus in the elastic strain stage.
Continuing to load, the rock enters the decay creep stage, and the creep curve at this stage has obvious nonlinear characteristics; therefore, it is difficult to accurately describe its characteristics by using the classical element model. In this paper, the viscosity coefficient of the sticky pot element in the traditional Kelvin model is modified by introducing the time-function relationship of the damage variable. That is, it is assumed that the viscosity coefficient has an exponential function relationship with the creep time during the rock creep process, so as to construct a nonlinear Kelvin model, and its three-dimensional differential constitutive equation satisfies: where σ ij is the stress, G 2 , η 2 are the shear modulus and viscosity coefficient of the rock in the nonlinear viscoelastic stage, ε nve is the nonlinear viscoelastic strain, . ε nve is the first derivative of the nonlinear viscoelastic strain with time, λ is an undetermined constant, and t is creep time. Separating the variables of the above differential equation and integrating with the initial conditions t = 0 , . ε nve = 0, the creep equation of the nonlinear Kelvin body can be obtained as: As the stress continues to be loaded, viscoelastic strain occurs in the rock. In this paper, the Kelvin model is used to represent this stage in the rheological process of the rock. The constitutive equation is: where: ε ve ij is the viscoelastic strain under true triaxial stress, G 3 is the shear modulus in the viscoelastic strain stage, and η 3 is the viscosity coefficient in the viscoelastic strain stage.
When the applied stress level exceeds the long-term strength of the rock, the rock enters the accelerated creep stage. In the case of three-dimensional stress, if the stress state exceeds the viscoplastic yield surface, viscoplastic strain will begin to occur. The rate of change of viscoplastic strain is represented by the limit stress flow law of Perzyna [26,27]: where: where F is the yield function of the rock, Q is the plastic potential function, and η 4 is the time-dependent viscosity coefficient. Based on Equations (9), (15), and (16), the three-dimensional damage constitutive equation in the accelerated creep state can be obtained [28]: The constitutive models of the four stages are combined, as shown in Figure 1. where: where F is the yield function of the rock, Q is the plastic potential function, and η4 is the time-dependent viscosity coefficient.
Based on Equations (9), (15), and (16), the three-dimensional damage constitutive equation in the accelerated creep state can be obtained [28]: The constitutive models of the four stages are combined, as shown in Figure 1.
In a true triaxial stress environment, the rock is stressed in three directions. For true triaxial rock creep tests, there are: In the process of triaxial compression, the failure mode of rock is mainly the compression shear failure sliding along the failure surface. Due to the formation and further development of cracks in the rock and the expansion of internal pores, the viscosity coefficient of the rock will be affected. Therefore, considering the influence of damage on shear modulus and viscosity coefficient, ignoring its influence on bulk modulus [29], there are: Combining Equations (18)- (20), the creep equation of the viscoelastic-plastic chemical damage creep model of rock under true triaxial stress state can be obtained as:  [31] analyzed multiple sets of rock true triaxial compression test data and found that the linear criterion proposed by Mogi was suitable for the true triaxial compression test of rock; this criterion is called the Mogi-Coulomb criterion, and its expression is as follows: where: where τ oct is the octahedral stress, σ m,2 is the effective normal stress, and δ 1 , δ 2 are experimental constants, which can be obtained by data fitting.
When the stress state of the rock reaches its stress yield surface, the rock will have obvious time-dependent behavior. Based on the Mogi-Coulomb criterion under threedimensional stress, the damage stress yield surface of the rock in three-dimensional space is obtained as: where: where: I 1 is the first invariant of stress, J 2 is the second invariant of deviatoric stress, and θ is the Lode angle of stress.

Parameter Identification
Due to the lack of laboratory creep test data considering the simultaneous action of acid corrosion and true triaxial stress, the conventional triaxial rock creep test after acid solution corrosion and the rock creep test under true triaxial stress were used for parameter identification and model verification.
Wang [21] conducted a conventional triaxial creep test of sandstone after acid corrosion. From this experiment, we observed that the elastic modulus of rock decreased obviously after the action of different pH values; this is the embodiment of chemical damage. We take the experimental data when the confining pressure is 6 MPa and pH = 3 to verify, the stress and strain are shown in the Figure 2  According to Equation (7), we can obtain the value of the Dc when pH = 3 is 0.312. Wang conducted a conventional triaxial creep test of rock, which is a special case of σ2 = σ3. Therefore, the creep constitutive model deduced in this paper can be greatly simplified.
Taking the long-term strength σs and strain rate ε  during the rock creep process as the judgment conditions for each stage of the constitutive equation, the models that need to be determined are G1, K1, λ, G2, η2, G3, η3, η4, and N. The creep failure time tF of the rock can be determined according to the test data. G1, K1 are determined according to the elastic modulus E and Poisson's ratio v obtained from the conventional triaxial compression test of the rock under the same confining pressure, namely: When in Formula (21) t → ∞ , we obtain: Subtract the first Formula of (21) from Formula (30), and take the logarithm of both sides to obtain: Then order: The Formula (31) becomes: Finally, the nonlinear fitting analysis is performed on the experimental data, and the fitting parameters a, b, and c can be obtained, then According to Equation (7), we can obtain the value of the Dc when pH = 3 is 0.312. Wang conducted a conventional triaxial creep test of rock, which is a special case of σ 2 = σ 3 . Therefore, the creep constitutive model deduced in this paper can be greatly simplified.
Taking the long-term strength σ s and strain rate · ε during the rock creep process as the judgment conditions for each stage of the constitutive equation, the models that need to be determined are G 1 , K 1 , λ, G 2 , η 2 , G 3 , η 3 , η 4 , and N. The creep failure time t F of the rock can be determined according to the test data. G 1 , K 1 are determined according to the elastic modulus E and Poisson's ratio v obtained from the conventional triaxial compression test of the rock under the same confining pressure, namely: When in Formula (21) t → ∞ , we obtain: Subtract the first Formula of (21) from Formula (30), and take the logarithm of both sides to obtain: Then order: The Formula (31) becomes: Finally, the nonlinear fitting analysis is performed on the experimental data, and the fitting parameters a, b, and c can be obtained, then λ, G 2 , η 2 are: The creep failure time t F of the rock can be determined according to the rock creep test, and the remaining model parameters G 3 , η 3 , η 4 , and N are determined by inversion based on the principle of least squares. The Levenberg-Marquardt algorithm that comes with Origin software is used for nonlinear regression analysis to identify the parameters of the rock creep model. The algorithm improves the least squares method and introduces a damping factor d to avoid iterative non-convergence [32]. Table 1 shows the rock creep model parameters considering chemical damage. Figure 3 shows the comparison between the rock creep test curve and the constitutive curve considering chemical damage.
The creep failure time tF of the rock can be determined according to the rock creep test, and the remaining model parameters G3, η3, η4, and N are determined by inversion based on the principle of least squares. The Levenberg-Marquardt algorithm that comes with Origin software is used for nonlinear regression analysis to identify the parameters of the rock creep model. The algorithm improves the least squares method and introduces a damping factor d to avoid iterative non-convergence [32]. Table 1 shows the rock creep model parameters considering chemical damage. Figure 3 shows the comparison between the rock creep test curve and the constitutive curve considering chemical damage.    Zhao [17] carried out the creep test of rock under true triaxial stress with Jinping marble, and obtained the mechanical behavior of rock creep under different intermediate principal stresses, which can be regarded as a special case where the chemical damage factor is 0. Therefore, it is used to verify the applicability of the constitutive model derived in this paper under the true triaxial stress state. The creep data under the action of true triaxial stress are shown in Table 2 below. According to the data, we can obtain the linear relationship between the octahedral shear stress τoct and the average normal stress σ m,2 , and thus obtain the size of the constant in formula (22), as shown in Figure 4, C 1 = 0.608, C 2 = 17.2, and substituting it into formula (25) to obtain the expression of the yield surface.   0  60  253  5  5  225  5  30  265  5  50  281  5  100  306  10  50  325  15  15  285  15  30  289  20  50  356  30  30  348  30  40  373  30  50  384  30  80  409  30  105  413  30  120  425  30  150  441  40  50  415  40  100  461  40 200 511 The determination of the parameters of the constitutive model under the act true triaxial stress is consistent with the above process. In the case of three-dimen stress, this paper selects rock creep test data under different intermediate pri stresses and rock creep data under different maximum principal stresses, respective verification. The creep parameters of marble under the true triaxial stress state we tained as shown in Table 3. The experiments conducted by Zhao did not consid chemical damage of the rock. Therefore, for the rock creep equation derived in this considering chemical damage, the chemical damage factor Dc is 0. Combined wi yield surface constants C1, C2 obtained above and the determined parameters  The determination of the parameters of the constitutive model under the action of true triaxial stress is consistent with the above process. In the case of three-dimensional stress, this paper selects rock creep test data under different intermediate principal stresses and rock creep data under different maximum principal stresses, respectively, for verification. The creep parameters of marble under the true triaxial stress state were obtained as shown in Table 3. The experiments conducted by Zhao did not consider the chemical damage of the rock. Therefore, for the rock creep equation derived in this paper considering chemical damage, the chemical damage factor D c is 0. Combined with the yield surface constants C 1 , C 2 obtained above and the determined parameters of the constitutive model, the comparison between the test curve and the model curve of the rock under the action of true triaxial stress is shown in Figure 5.

Constitutive Model Verification
Comparing the experimental and theoretical curves in Figures 3 and 5, we can conclude that the true triaxial creep model curve under the action of rock acid corrosion deduced in this paper is in agreement with the experimental curve, which not only reflects the effect of chemical damage on rock creep characteristics, but also reflects the creep behavior of rock under the action of true triaxial stress, including elastic stage, decay stage, constant velocity stage, and acceleration stage of rock creep. Different intermediate principal stresses will have an important influence on the creep behavior of rocks. It can be seen that when the intermediate principal stress is too large or too small, the conformity between the constitutive curve and the test curve will decrease to a certain extent. At the same time, it can be seen that when the rock creep is in the acceleration stage, the data of the constitutive curve and the test curve often have a certain deviation, and the deviation is about 10%. This is because the stress on the rock at this stage is greater than its long-term strength, and the cracks generated on the surface and inside of the rock begin to extend and penetrate, resulting in uneven distribution of stress inside the rock at this stage. On the whole, the true triaxial creep model considering chemical corrosion established in this paper can accurately reflect the creep behavior of rock, which verifies the rationality of the model establishment and the correctness of the parameter verification. constitutive model, the comparison between the test curve and the model curve of the rock under the action of true triaxial stress is shown in Figure 5.

Constitutive Model Verification
Comparing the experimental and theoretical curves in Figures 3 and 5, we can conclude that the true triaxial creep model curve under the action of rock acid corrosion deduced in this paper is in agreement with the experimental curve, which not only reflects the effect of chemical damage on rock creep characteristics, but also reflects the creep behavior of rock under the action of true triaxial stress, including elastic stage, decay stage, constant velocity stage, and acceleration stage of rock creep. Different intermediate principal stresses will have an important influence on the creep behavior of rocks. It can be seen that when the intermediate principal stress is too large or too small, the conformity between the constitutive curve and the test curve will decrease to a certain extent. At the same time, it can be seen that when the rock creep is in the acceleration stage, the data of the constitutive curve and the test curve often have a certain deviation, and the deviation is about 10%. This is because the stress on the rock at this stage is greater than its longterm strength, and the cracks generated on the surface and inside of the rock begin to extend and penetrate, resulting in uneven distribution of stress inside the rock at this stage. On the whole, the true triaxial creep model considering chemical corrosion established in this paper can accurately reflect the creep behavior of rock, which verifies the rationality of the model establishment and the correctness of the parameter verification.

Influence of Model Parameters on Intermediate Principal Stress
The rock creep constitutive model established in this paper considers the effect of true triaxial stress. Therefore, it is necessary to explore the influence of different intermediate principal stresses on the model parameters. The creep model parameters of rock under different intermediate principal stresses are shown in Figure 6

Influence of Model Parameters on Intermediate Principal Stress
The rock creep constitutive model established in this paper considers the effect of true triaxial stress. Therefore, it is necessary to explore the influence of different intermediate As can be seen from Figure 6, the value of the shear modulus G1, as an elastic parameter, increases gradually with the increase in the intermediate principal stress. This shows that with the increase in the intermediate principal stress, the original cracks, holes, and other micro-defects in the rock are closed, so that the rock sample is hardened to a certain extent. The viscosity coefficients η2, η3 are used to describe the ability of rock samples to resist deformation, and their values show an upward trend with the increase in the intermediate principal stress, which means that the increase in the intermediate principal stress inhibits the time-dependent deformation of the rock, so at the macro level, it is manifested as a decrease in the creep rate; meanwhile, the material parameter shows a decreasing As can be seen from Figure 6, the value of the shear modulus G 1, as an elastic parameter, increases gradually with the increase in the intermediate principal stress. This shows that with the increase in the intermediate principal stress, the original cracks, holes, and other micro-defects in the rock are closed, so that the rock sample is hardened to a certain extent. The viscosity coefficients η 2 , η 3 are used to describe the ability of rock samples to resist deformation, and their values show an upward trend with the increase in the intermediate principal stress, which means that the increase in the intermediate principal stress inhibits the time-dependent deformation of the rock, so at the macro level, it is manifested as a decrease in the creep rate; meanwhile, the material parameter λ shows a decreasing trend with the increase in the intermediate principal stress.

Conclusions
Based on damage mechanics, chemical damage is connected with elastic modulus; thus, the damage relations considering creep stress damage and chemical damage are obtained. The elastic body, nonlinear Kelvin body, linear Kelvin body, and viscoelasticplastic body are connected in series, and the actual situation under the action of true triaxial stress is considered at the same time. Therefore, a damage creep constitutive model considering the coupling of rock acid corrosion and true triaxial stress is established, which is verified with the experimental data, and the following conclusions are obtained: 1.
The constitutive model established in this paper can accurately reflect the creep characteristics of rock under acid corrosion and true triaxial stress state, and the model's fitting degree is more than 90%, which verifies its rationality and the correctness of parameter determination. Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.