Abstract
Based on the Kachanov damage theory and elastic wave theory, considering the long-term creep damage from time dimension and instantaneous disturbance damage from space dimension, we presented a 3D damage creep model and converted it to difference expressions in order to write into the finite difference software FLAC3D. Then, according to the results of creep tests, we conducted parameter inversion of our damage creep model with the help of the particle swarm optimization (PSO) method. Finally, the damage creep model was applied in a railway tunnel project in Yunnan to simulate the tunnel deformation. Compared with the Burgers model and the model considering only creep damage, our model which considers both creep damage and disturbance damage yielded more reasonable results.
1. Introduction
Any material under long-term load tends to creep, and rock is no exception, especially if it contains many internal flaws which could exacerbate such process. In geotechnical engineering, the creep property of rock has been widely studied in the past 80 years, including lots of mechanic experiments [1–5] and constitutive models [1, 6–9]. There are many component creep models such as Kelvin, Maxwell, and Burgers which have been applied for simulating the mechanical behavior of rock masses under long-term load, and many improvements were carried out on these traditional creep models in order to improve their accuracy and applicability. One of the improvements is rock damage, which has received great attention since Kachanov proposed the damage theory in 1958 [10] [11–14]. Because it was found that the mechanical property of rock would continuously deteriorate during the creep process, without considering the rock damage, the simulated results showed a lack of security.
Combining the rock creep model and damage theory, many representative results have been obtained. From 1992 to 1997, Chan et al. [15–18] conducted triaxial compressive creep tests on rock salt, considering the contribution of damage to the inelastic strain, and developed a constitutive model for inelastic flow and damage evolution. From 1991 to 2000, Aubertin et al. [19–21] built a viscoplastic-damage model for soft rocks with low porosity on the basis of an internal (or state) variable and made improvement on it. From 2003 to 2010, Shao et al. [22–24] presented a constitutive model for creep deformation considering the progressive degradation of elastic modulus and failure strength and then it coupled with a second-order damage tensor for the description of induced anisotropic damage in brittle rocks. But the damage consists of two parts: long-term creep damage and the instantaneous disturbance damage, which impact the creep property.
So based on the traditional Burgers creep model, the influence of both the creep damage from time dimension and disturbance damage from space dimension will be considered in this paper. Meanwhile, the improved creep model will be written into the finite difference software and a rapid inversion of model parameters will be achieved by using the PSO method. Finally, our damage creep model will be tested in a tunnel project and will be compared with the Burgers model and the model only considering creep damage.
2. Creep Damage
The traditional Burgers model is composed of a Kelvin model and a Maxwell model connected in series (dashed box in Figure 1), which can simulate the transient creep stage and steady-state creep stage but cannot simulate the tertiary creep stage.

In order to improve this defect, a damage element can be added into the model to describe the time-dependent creep damage (Figure 1). Based on the Kachanov damage theory [25], the creep damage factor Dc can be calculated using equation (1); thus, the viscosity coefficient of damage element can be calculated using equation (2):where is the failure time and is the coefficient related to damage process.
The improved Burgers model with a damage element can be shown as the following equations. When the stress is lower than the long-term strength , the model is the same as the Burgers model, and when is beyond the long-term strength , the damage element will start to work. It can be seen form Figure 2 that the viscosity coefficient of damage element decreasing with time forms to 0, and it leads to a sharp increase of strain before failure time , which can properly describe the tertiary creep stage:where is the strain, and are the parameters of the Maxwell model, and and are the parameters of the Kelvin model.

3. Disturbance Damage
The most common disturbance damage is caused by excavation. It can be considered as an instantaneous damage, so it is a function of position instead of time. As we all know, the closer the rock mass is to the excavation surface, the greater the disturbance damage is, so the damage should be related to the distance away from the excavation surface.
Usually, the elasticity modulus E can be used to reflect the damage of rock. Based on the elastic wave theory (equation (4)) [26, 27], if we assume that density of rock and Poisson’s ratio are constants, E should be proportional to . So we can measure the longitudinal wave velocity to obtain the change rule of E:
By analyzing a large number of measurement results, it is found that the curves of are basically consistent, which are shaped like ‘S’ (Figure 3). The S-curve is monotonically increasing which means E is increasing and damage is decreasing, until the distance away from the excavation surface is large enough to ignore the disturbance damage and the curve tends to be stable. And it can also be seen from the curve that the increase rate changes from high to low, so we can use equation (5) to fit the S-curve, and consequently, E can be written as equation (6):where A, B, C, and D are the fitting parameters and x is the distance away from the excavation surface.

In conclusion, the disturbance damage factor Dd can be calculated as follows:
4. 3D Damage Creep Model
We assume the volume deformation is plastic, and the creep property is mainly correlated with the shear deformation. So the deviator strain of Maxwell model, Kelvin model, damage model, and total deviator strain can be calculated by the following equations:
The spherical strain can be calculated using the following equation:
From all the above equations, the 3D damage creep model can be written as follows:where is the deviator strain rate of the Maxwell model, is the deviator strain rate of the Kelvin model, is the deviator strain rate of the damage element, is the total deviator strain rate, is the spherical strain rate, is the deviator stress, is the spherical stress, is the shear modulus of the Maxwell model, is the coefficient of viscosity of the Maxwell model, is the shear modulus of the Kelvin model, is the coefficient of viscosity of the Kelvin model, is the coefficient of viscosity of the damage element, is the bulk modulus, is the yield function, is the initial value of yield function, is the total strain, is the Kronecker delta, and is the second invariant of deviator stress tensor.
In order to be written into FLAC3D, all the equations need to be converted to the difference form:where new means the new value in a time step and old means the old value in a time step.
5. Creep Tests
A group of excavated rock samples were taken from a tunnel excavation project in Yunnan, China, and creep tests were conducted on them (Figure 4). In the test, the confining pressure remained at 10 MPa and the loading process was divided into 5 stages: 50 MPa, 70 MPa, 90 MPa, 110 MPa, and 130 MPa. Each stage was maintained for 24 h and finally the creep curves were obtained (Figure 5). Also, we obtained the long-term strength according to the isochronous creep curve.

(a)

(b)

After that, 5 groups of destructive tests under 100 MPa, 110 MPa, 120 MPa, 150 MPa, and 180 MPa, respectively, were conducted to obtain the failure time . From Figure 6, it can be seen that the failure time is approximately inversely proportional to the stress. So in this case, by fitting the curve, we assume that the failure time can be calculated by the following equation:

Thus, the creep damage factor can be written as
6. Parameter Inversion Based on PSO
According to the longitudinal wave velocity, we obtained the fitting parameters B = 0.31, C = 3.05, and D = 1.06. And our rock samples were taken from the evacuation surface, which means x = 0, so the disturbance damage factor can be obtained as follows:
In our triaxial compression test, , so , , and , and the 3D damage creep model can be converted to the following equations:
By fitting the creep curves, there are seven parameters, , to be determined; it must be a huge quantity of computation. So we constructed an objective function F (equation (19)) and used the PSO method to search for the optimum value:where is the number of data points, is the strain calculated by 3D damage creep model, and is the strain of creep curves.
The PSO method uses particles to simulate the social behavior of birds [28]. In an n-dimensional space (in this case, it is 7-dimensional space), the position of each particle is a possible solution to the problem, and each particle has a fitness value to evaluate its position (in this case, it is F). During the searching process, each particle adjusts its velocity according to its own best position and the group’s best position (equation (20)), and then, the position will be updated (equation (21)). After repeated iteration, particles can finally move to the optimal position:where s is the number of steps, is the velocity of particle, is the position of particle, i means the i-th particle, j means the j-th dimension, pbest means the particle’s best position, gbest means the group’s best position, r1 and r2 are the independent random numbers in the range of (0, 1), and c1 and c2 are the acceleration parameters.
With the help of the PSO method, after just about 1500 iterations (Figure 7), we yielded a satisfactory result. The seven fitted parameters are listed in Table 1, and the fitting curves are shown in Figure 8; compared with the testing curve, it shows high accuracy.


In order to validate the 3D damage creep model in FLAC3D, we rebuilt a cylindrical rock sample which is completely the same as the tested rock sample, including the size, parameters, constraints, and loadings. As shown in Figure 9, the result was consistent with the test data and the difference was always less than 10%, which proved the feasibility of the model.

7. Application and Discussion
Figure 10 shows the profile of the tunnel in Yunnan. On the basis of FLAC3D, we used the Burgers model, the model only considering creep damage, and the model considering both creep damage and disturbance damage to simulate the tunnel, respectively. To consider the disturbance damage, we built a function to assign different values of disturbance factor to each grid cell on the basis of different distances to excavation surface. Finally, the simulated results are listed in Table 2, and also, the arch settlement and haunch convergence of the three models are drawn in Figure 11 in order to compare with each other.


The results show that there was no plastic zone using the Burgers model, and the displacement was far less than in monitoring data, even just a quarter of the monitoring data. When we considered the creep damage, an obvious plastic zone appeared, and the displacement increased significantly, but there was still a big gap to the monitoring data. When both creep damage and disturbance damage were considered, the plastic zone continued to increase, especially the roof of tunnel. The displacement also increased, and the maximum difference with monitoring data was just about 4 mm.
In conclusion, the creep damage and disturbance damage both have a great impact on the simulated results, so they are necessary to be considered to ensure safety.
8. Conclusions
(1)The traditional Burgers model cannot simulate the tertiary creep stage. Adding a damage element based on the creep damage factor can properly describe the tertiary creep stage.(2)The disturbance damage is a function of the distance away from the excavation surface, according to the elastic wave theory and the measurement results of longitudinal wave velocity. The disturbance damage factor can be written as .(3)The 3D damage creep model is(4)In this case, the failure time is approximately inversely proportional to the stress.(5)The PSO method shows high accuracy and high efficiency in the process of searching for optimum fitting parameters.(6)There was no plastic zone using the Burgers model, and the displacement was far less than monitoring data. If we considered the creep damage and disturbance damage, there will be an obvious plastic zone, and the displacement will be more close to the monitoring data.(7)The creep damage and disturbance damage both have a great impact on the simulated results, so they are necessary to be considered to ensure safety.
Data Availability
All the data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was funded by the National Natural Science Foundation of China (No. 51978106).