# A 3D Creep Model considering Disturbance Damage and Creep Damage and Its Application in Tunnel Engineering

Authors: Geng-Feng Wang; Fan Sun; Xiao-Hui Xiong; Lei He; Ke-Hong Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 Figure1), which can simulate the transient creep stage and steady-state creep stage but cannot simulate the tertiary creep stage.Figure 1
Composition of the creep model considering creep damage.In order to improve this defect, a damage element can be added into the model to describe the time-dependent creep damage (Figure1). 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 η3t can be calculated using equation (2):(1)Dc=1−1−ttF1/α,(2)η3t=η3×1−ttF1/α,where tF 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 σs, the model is the same as the Burgers model, and when σ is beyond the long-term strength σs, 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 η3 to 0, and it leads to a sharp increase of strain before failure time tF, which can properly describe the tertiary creep stage:(3)εt=σE1+ση1t+σE21−e−E2/η2tσ<σs,σE1+ση1t+σE21−e−E2/η2t+σ−σsη3×1−t/tF1/αtσ≥σs,where ε is the strain, E1 and η1 are the parameters of the Maxwell model, and E2 and η2 are the parameters of the Kelvin model.Figure 2
Creep curves of the creep model considering creep damage.


## 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 modulusE 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 vp2. So we can measure the longitudinal wave velocity vp to obtain the change rule of E:(4)vp=E1−μρ1+μ1−2μ.By analyzing a large number of measurement results, it is found that the curves ofvp 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):(5)vp=A1−Be−ex−C/D,(6)Ex=E1−Be−ex−C/D2,where A, B, C, and D are the fitting parameters and x is the distance away from the excavation surface.Figure 3
Curve of longitudinal wave velocity in rock.In conclusion, the disturbance damage factorDd can be calculated as follows:(7)Dd=1−1−Be−ex−C/D2.


## 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:(8)e˙1,ij=s˙ij2G1x+sij2η1x,(9)sij=2η2xe˙2,ij+2G2xe2,ij,(10)e˙3,ij=F/F02η3t,x∂F∂σij,(11)e˙ij=e˙1,ij+e˙2,ij+e˙3,ij.The spherical strain can be calculated using the following equation:(12)σ˙m=3Kxε˙m.From all the above equations, the 3D damage creep model can be written as follows:(13)εijt=σmδij3Kx+sij2G1x+sijt2η1x+sij2G2x1−e−G2x/η2xtF<0,σmδij3Kx+sij2G1x+sijt2η1x+sij2G2x1−e−G2x/η2xt+F/F0∂F/∂σij2η3t,xtF≥0,F=J2−σs3,where e˙1,ij is the deviator strain rate of the Maxwell model, e˙2,ij is the deviator strain rate of the Kelvin model, e˙3,ij is the deviator strain rate of the damage element, e˙ij is the total deviator strain rate, ε˙m is the spherical strain rate, sij is the deviator stress, σm is the spherical stress, G1x is the shear modulus of the Maxwell model, η1x is the coefficient of viscosity of the Maxwell model, G2x is the shear modulus of the Kelvin model, η2x is the coefficient of viscosity of the Kelvin model, η3t,x is the coefficient of viscosity of the damage element, Kx is the bulk modulus, F is the yield function, F0 is the initial value of yield function, εijt is the total strain, δij is the Kronecker delta, and J2 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:(14)Δe1,ij=Δsij2G1x+sijnew+sijold4η1xΔt,sijnew+sijold2Δt=2η2xΔe2,ij+2G2xe2,ijnew+e2,ijold2Δt,Δe3,ij=F/F02η3t,x∂F∂σijΔt,Δeij=Δe1,ij+Δe2,ij+Δe3,ij,Δσm=3KxΔεm,J2=sijsij2,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 (Figure4). 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 σs=95MPa according to the isochronous creep curve.Figure 4
Rock samples and experimental apparatus.
(a)
(b)Figure 5
Experimental creep curves.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 timetF. From Figure 6, it can be seen that the failure time tF is approximately inversely proportional to the stress. So in this case, by fitting the curve, we assume that the failure time tF can be calculated by the following equation:(15)tF=356σ−σs.Figure 6
Experimental failure time curve.Thus, the creep damage factor can be written as(16)Dc=1−1−σ−95356t1/α.


## 6. Parameter Inversion Based on PSO

According to the longitudinal wave velocity, we obtained the fitting parametersB = 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:(17)Dd=1−1−Be−ex−C/D2=1−1−0.31e−e−3.05/1.062=0.5.In our triaxial compression test,σ2=σ3=10MPa, so σm=1/3σ1+2σ3, s11=2/3σ1−σ3, and J2=1/3σ1−σ3, and the 3D damage creep model can be converted to the following equations:(18)ε11t=σ1+2σ34.5K+σ1−σ31.5G1+σ1−σ31.5η1t+σ1−σ31.5G21−e−G2/η2tσ1−σ3<95,σ1+2σ34.5K+σ1−σ31.5G1+σ1−σ31.5η1t+σ1−σ31.5G21−e−G2/η2t+σ1−σ3−953η3ttσ1−σ3≥95,η3t=0.5η3×1−σ1−σ3−95356t1/α.By fitting the creep curves, there are seven parameters,K,G1,G2,η1,η2,η3,andα, 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:(19)FK,G1,G2,η1,η2,η3,b=∑i=1nεi˜−εi2,where n is the number of data points, εi˜ is the strain calculated by 3D damage creep model, and εi 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:(20)vijs+1=vijs+c1⋅r1ijs⋅xpbestijs−xijs+c2⋅r2ijs⋅xgbestjs−xijs,(21)xijs+1=xijs+vijs+1,where s is the number of steps, v is the velocity of particle, x 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 (Figure7), 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.Figure 7
Searching process using PSO.Table 1
Fitting parameters of the damage creep model.
K
/
Gpa
G
1
/
Gpa
G
2
/
Gpa
η
1
/
Gpa
⋅
h
η
2
/
Gpa
⋅
h
η
3
/
Gpa
⋅
h
α
68
45
29
9800
0.36
2100
3.8Figure 8
Fitting curves of experimental creep curves.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.Figure 9
Validation of the 3D damage creep model in FLAC3D.


## 7. Application and Discussion

Figure10 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.Figure 10
The profile of the tunnel in Yunnan, China.Table 2
Simulated results of the tunnel.
Model
Burgers
Creep damage
Creep damage and disturbance damage
Plastic zone
Arch settlement
Haunch convergenceFigure 11
Tunnel deformation simulated by different models.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 factorDc=1−1−t/tF1/α 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 asDd=1−1−Be−ex−C/D2.
(3)
The 3D damage creep model is
(22)
ε
i
j
t
=
σ
m
δ
i
j
3
K
x
+
s
i
j
2
G
1
x
+
s
i
j
t
2
η
1
x
+
s
i
j
2
G
2
x
1
−
e
−
G
2
x
/
η
2
x
t
F
<
0
,
σ
m
δ
i
j
3
K
x
+
s
i
j
2
G
1
x
+
s
i
j
t
2
η
1
x
+
s
i
j
2
G
2
x
1
−
e
−
G
2
x
/
η
2
x
t
+
F
/
F
0
∂
F
/
∂
σ
i
j
2
η
3
t
,
x
t
F
≥
0
,
(4)
In this case, the failure timetF 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.



---

*Source: 7073089-2020-04-21.xml*