# Shear Creep Simulation of Structural Plane of Rock Mass Based on Discontinuous Deformation Analysis

Authors: Guoxin Zhang; Zhengqi Lei; Heng Cheng
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Numerical simulations of the creep characteristics of the structural plane of rock mass are very useful. However, most existing simulation methods are based on continuum mechanics and hence are unsuitable in the case of large displacements and deformations. The discontinuous deformation analysis method proposed by Genhua is a discrete one and has a significant advantage when simulating the contacting problem of blocks. In this study, we combined the viscoelastic rheological model of Burgers with the discontinuous deformation analysis (DDA) method. We also derived the recurrence formula for the creep deformation increment with the time step during numerical simulations. Based on the minimum potential energy principle, the general equilibrium equation was derived, and the shear creep deformation in the structural plane was considered. A numerical program was also developed and its effectiveness was confirmed based on the curves obtained by the creep test of the structural plane of a rock mass under different stress levels. Finally, the program was used to analyze the mechanism responsible for the creep features of the structural plane in the case of the toppling deformation of the rock slope. The results showed that the extended DDA method is an effective one.


---

## Body

## 1. Introduction

The rheological properties of a rock mass are among its inherent mechanical properties. In the case of weak rock and rock masses with a clay interlayer and a fracture zone in the fault, the rheological properties are even more significant; for hard rock, if there are a large number of joints and fissures, the shear creep deformation can be quite high [1]. The rheological properties determine the deformation and instability of rock mass engineering underground as well as that of rock slope engineering. As an example, consider a nonferrous metal mine in the northern part of the Qinghai-Tibet Plateau. After the excavation of the tunnel, the surrounding rock deformed continuously, and the cumulative deformation was as high as 400 mm after 140 days. The maximum deformation and convergence value was 1000 mm. This had an adverse effect on mining and lane development [2]. In another case, a long and deep-buried tunnel in Wushaoling crosses Qilian Mountain in China. Owing to the combined effects of the high ground stress and excavation disturbances, the rock surrounding the tunnel has undergone continuous rheological deformation. The maximum settlement of the vault is 1053 mm, and the deformation is not yet convergent [3]. The left abutment of Longtan Hydropower Station exhibits typical toppling rock creep. The total volume of the deformable body is approximately 11.90 million m3. Further, the long-term creep deformation has resulted in multiple instances of collapse and cracking and is a direct threat to the safe operation of the power station [4]. Similarly, after a huge hydropower station in the upper reaches of Yellow River started retaining water, the resulting landslide body began to show creep deformation toward the bed. The total volume of the dam slope on the right bank is tens of millions m3. Moreover, the measured maximum deformation rate is as high as 20–30 mm/day, and the cumulative deformation has reached 40 m, resulting in a high safety risk. The Tanggudong side slope is located in the basin of Yalong River in southwest China and consists of weathered sand slate from the Trias. Long-term river-water-related erosion and creep deformation resulted in a landslide body of 68 million m3 sliding into the river. This resulted in a giant landslide dam that blocked Yalong River for 9 days. The flood peak after dam break was up to 5.7 × 104 m3/s [5]. Danba county in Sichuan province in China is located in an ancient landslide body. Owing to long-term creep deformation and the effects of artificial disturbances, a landslide occurred in the county on February 18, 2005, resulting in damage to a large number of houses and total losses of up to 10.66 million yuan. This event was a great threat to the safety of the entire county [6]. These cases highlight the fact that it is essential to elucidate the rheological properties of rock mass in order to be able to predict the deformation of engineering rock mass and perform long-term safety assessments.Under the effect of long-term loads, the rock block and the structural plane constituting the rock mass both undergo rheological deformations. The aging-related mechanical behavior of rock mass is the result of these interactions [7, 8]. Experimental studies and numerical simulations on the rheological properties of rock block have yielded useful results [1, 9–15]. Zhu et al. [16] studied the greensand rock of the diversion tunnel in Jinping Second Cascade Hydropower Station. They performed rheological tests in the laboratory based on the stress path for constant axial compression and progressive unloading confining pressure. Based on the experimental results, the damage evolution equation as well as the variable-parameter nonlinear Burgers model could be established. This model can also be used to determine the decay creep stage and steady rheological behavior of green sandstone under unloading. Previous studies [17] have modified the anisotropic damage model to describe brittle rock. The aging deformation of brittle rock is believed to be the result of microcrack propagation caused by stress-related corrosion. Thus, the creep deformation of brittle rock can be simulated using the extended model. A three-axis shear creep test was performed by Li et al. under a high confining pressure and reunloading [18]. The transformation regularity of the stress state of soft rock in deep mines during tunnel excavation has been studied. The Mogi-Coulomb yield surface rheological constitutive model considering the volume deformation is proposed. This model has been used in a finite element simulation of the supporting process for tunnel excavation of a deep mine in Huainan. Moghadam et al. [19] proposed an elastic and plastic constitutive model that simulates the phenomena of expansion and strength evolution during salt rock creep and used the model to analyze aging deformation related to underground gas storage. Further, using finite element simulations, the effects of cavern excavation and the nonsalt rock interlayer on the large-scale displacement resulting from underground gas storage have been studied. Yang [20] performed numerical simulations of the flow deformation during the three-dimensional excavation for a dam slope at Dagangshan Hydropower Station. The simulations were based on the results of indoor and field rheological tests and used the nonlinear viscoelastic rheological model for hard brittle rock and a rheological model of variable-parameter damage for the weak rock zone. The rheological excavation behavior of rock slope under complex stress has also been investigated.The deformation and strength characteristics of rock mass that is rich in joints, fissures, and weak structural planes are generally determined by the structural plane. The rheological characteristics of the structural plane are more complex. In addition to external factors, such as the stress level, water content, and temperature, the production form of the plane, its filling material, and degree of weathering also play a critical role. However, there have been few studies on the rheological properties of the structural plane. The time-history curves of the shear stress-shear displacement of the weak interlayer in rock mass have been determined based on a laboratory-based shear rheological test by Zhu et al. [21]. Further, the long-term strength parameters of rock mass have also been determined. The standard linear rheological model has been obtained through model fitting. The creep characteristics of structural planes with different angles were studied under different stress conditions using the artificial tooth profile by Zhang et al. of Tongji University [22]. An improved Burgers model has been proposed based on the characteristics of the shear creep of the structural plane. A generalized Kelvin model and the Burgers model were used by Ding et al. [23] to describe the rheological behaviors of rock and weak structural planes, respectively. Numerical simulations of the uniaxial and triaxial compression of different types of rock mass have been performed using the software FLAC3D. The effects of the yield and number of structural planes on the creep behavior of rock mass have also been explored. The results showed that the yield of the structural planes not only changes the creep strength and displacement of rock mass but also controls its failure mode and condition. There have been studies [24] involving shear rheological tests of rock specimens with a weak structure surface. The existing Burgers model has been modified based on the brittle damage characteristics of the greenschist structure surface and the curves of shear creep tests. Generally speaking, most existing numerical simulation methods are based on the finite element method and the finite difference method. However, owing to the limitations of the methods themselves, the contact transformation between rock blocks cannot be simulated. Further, it is difficult to accurately simulate the large deformations and displacements caused by the shear creep of the structural plane.In the 1980s, Dr. Shi and Pei [25] proposed a numerical method based on discontinuous media mechanics, namely, the discontinuous deformation analysis (DDA) method. The method considers an arbitrary block as the basic unit and the rigid body displacement and deformation of the block as the basic unknown quantities. Based on a set of effective contact searching method, the DDA method can be used to determine arbitrary contact forms of the block system. At the same time, the method takes into account the interaction between the blocks. Based on the principle of minimum potential energy, a general equilibrium equation is established and an implicit solution is used to solve the problem at hand. The DDA method has an inherent advantage with respect to simulations of the contact problem of block systems [26–29]. In this study, we combine the DDA method and the Burgers rheological model to establish the recurrence formula for the increment in creep deformation for each calculation time step during numerical simulations. The general equilibrium equation for shear creep deformation is derived and the corresponding numerical program is developed. The effectiveness of the program is verified based on curves obtained from creep tests performed under different stress levels. Finally, the program is used to analyze the mechanism responsible for the effects that the structural plane creep characteristics have on the rock slope.


## 2. Rheological Model

### 2.1. Burgers Rheological Model

The rheological constitutive model of rock mass is a mathematical and physical model for describing the stress-strain-time relationship for the rock mass. Based on extensive experimental and theoretical studies, researchers have put forward a variety of rheological models [30]. The combination model is often used to simulate the rheological behavior of rock mass numerically. The basic principle of the combination model is to set the basic components, such as the elastic element, plastic element, and viscous element, based on the properties of the rock. These components can be combined into various constitutive models that reflect the rheological properties of different types of rock mass through series and parallel connections. Each component in the series connection shares the same load, and the sum of all the elements is equal to the total strain rate. In the case of parallel connections, the sum of all elements is equal to the total load, and the strain rate of each element is equal to one another.The Burgers model [30] is formed by connecting the Maxwell model and the Kelvin model in series. The model betters the description of the viscoelastic rheological characteristics of the structural plane. In this study, the shear creep characteristics of the structural plane as described by the model were adopted. The component diagram and creep curves are shown in Figures 1 and 2, respectively.Figure 1
Schematic diagram of the Burgers model.Figure 2
Creep curves for the Burgers model.The constitutive equation of the Burgers model is as follows:(1)η1ε˙+η1η2E2·ε¨=σ0+η1E1+η2E2+η1E2σ˙0+E1+E2E1E2·σ¨0.In this equation,E1, E2, η1, and η2 are the elastic moduli and viscosity coefficients of the elements shown in Figure 1, while σ0 is the normal stress and ε is the total strain.The creep equation for the Burgers model is given by(2)εt=σ0E1+σ0E21-e-E2/η2t+σ0η1t=σ0E1+σ0CtCt=σ0E21-e-E2/η2t+σ0η1t.In this equation,C(t) is the degree of creep and t is the time of operation of the stress, σ0.By considering Figure2 and (2) together, it can be seen that the creep equation for the Burgers model is composed of three parts: (3)ε1=σ0E1ε2=σ0E21-exp-E2η2tε3=σ0η1t,where ε1 is the instantaneous elastic strain; ε2 is the recoverable strain decaying with a negative exponent, meaning that it corresponds to the first section of the creep curve; and ε3 is the nonrecoverable strain corresponding to the steady strain rate, meaning that it is related to the second section of the creep curve.


### 2.2. Increment Expression of Rheology

In DDA method, the time,t, is usually divided into a series of time periods. The group of equilibrium equations is set up for each time interval and the stress increment in this period is determined. In addition, the creep deformation of the structural plane is related to the entire stress history. Therefore, in order to calculate the shear creep deformation for time t, it is necessary to establish the recurrence relation for the creep strain increment for each time step. As shown in Figure 3, the time, t, is divided into a series of moments (t0, t1, t2,…,tn-1, and tn), and it is assumed that the stress increments at these moments are Δσ0, Δσ1, Δσ2,…,Δσn-1, and Δσn, respectively. The three adjacent moments tn-1, tn, and tn+1 are selected as the target moments. Their time intervals are Δtn=tn-tn-1 and Δtn+1=tn+1-tn. Therefore, the creep deformations for the three adjacent moments can be given by (4)εctn-1=Δσ0Ctn-1-t0+Δσ1Ctn-1-t1+⋯+Δσn-2Ctn-1-tn-2(5)εctn=Δσ0Ctn-t0+Δσ1Ctn-t1+⋯+Δσn-1Ctn-tn-1(6)εctn+1=Δσ0Ctn+1-t0+Δσ1Ctn+1-t1+⋯+Δσn-1Ctn+1-tn-1+ΔσnCtn+1-tn.Figure 3
Stress increment during each time step of the calculation.If (5) is subtracted from (6), we get (7)Δεtn+1c=εctn+1-εctn=Δσ0Ctn+1-t0-Ctn-t0+Δσ1Ctn+1-t1-Ctn-t1+⋯+Δσn-1Ctn+1-tn-1-Ctn-tn-1+ΔσnCtn+1-tn.From (2), we get (8)Ctn+1-ti-Ctn-ti=1E21-e-E2/η2tn+1-ti-1E21-e-E2/η2tn-ti+1η1tn+1-ti-1η1tn-ti=1E2e-E2/η2tn-ti-e-E2/η2tn+1-ti+1η1Δtn+1=1E2e-E2/η2tn-ti1-e-E2/η2Δtn+1+1η1Δtn+1.By inserting (8) into (7), we get (9)Δεn+1c=1-e-E2/η2Δtn+1ωn+1+νn+1+ΔσnCtn+1-tn,where(10)ωn+1=Δσ0E2e-E2/η2tn-t0+Δσ1E2e-E2/η2tn-t1+⋯+Δσn-1E2e-E2/η2tn-tn-1(11)υn+1=σ0η1Δtn+1+σ1η1Δtn+1+⋯+σn-1η1Δtn+1.Using the same principle, (4) is subtracted from (5): (12)Δεnc=1-e-E2/η2Δtnωn+νn+Δσn-1Ctn-tn-1(13)ωn=Δσ0E2e-E2/η2tn-1-t0+Δσ1E2e-E2/η2tn-1-t1+⋯+Δσn-2E2e-E2/η2tn-1-tn-2(14)υn=σ0η1Δtn+σ1η1Δtn+⋯+σn-2η1Δtn.By comparing (10) and (13), we get (15)ωn+1=ωne-E2/η2Δtn+Δσn-1E2e-E2/η2Δtnω1=Δσ0E2.Similarly, by comparing (11) and (14), we get(16)νn+1=vnΔtn+1Δtn+σn-1η1Δtn+1ν1=σ0η1Δt1.Equations (12)–(16) form a group of recurrence formulas that can be used to determine the rheological deformation increment for each time step of the calculation. In this manner, the rheological deformation at any moment can be determined.


## 2.1. Burgers Rheological Model

The rheological constitutive model of rock mass is a mathematical and physical model for describing the stress-strain-time relationship for the rock mass. Based on extensive experimental and theoretical studies, researchers have put forward a variety of rheological models [30]. The combination model is often used to simulate the rheological behavior of rock mass numerically. The basic principle of the combination model is to set the basic components, such as the elastic element, plastic element, and viscous element, based on the properties of the rock. These components can be combined into various constitutive models that reflect the rheological properties of different types of rock mass through series and parallel connections. Each component in the series connection shares the same load, and the sum of all the elements is equal to the total strain rate. In the case of parallel connections, the sum of all elements is equal to the total load, and the strain rate of each element is equal to one another.The Burgers model [30] is formed by connecting the Maxwell model and the Kelvin model in series. The model betters the description of the viscoelastic rheological characteristics of the structural plane. In this study, the shear creep characteristics of the structural plane as described by the model were adopted. The component diagram and creep curves are shown in Figures 1 and 2, respectively.Figure 1
Schematic diagram of the Burgers model.Figure 2
Creep curves for the Burgers model.The constitutive equation of the Burgers model is as follows:(1)η1ε˙+η1η2E2·ε¨=σ0+η1E1+η2E2+η1E2σ˙0+E1+E2E1E2·σ¨0.In this equation,E1, E2, η1, and η2 are the elastic moduli and viscosity coefficients of the elements shown in Figure 1, while σ0 is the normal stress and ε is the total strain.The creep equation for the Burgers model is given by(2)εt=σ0E1+σ0E21-e-E2/η2t+σ0η1t=σ0E1+σ0CtCt=σ0E21-e-E2/η2t+σ0η1t.In this equation,C(t) is the degree of creep and t is the time of operation of the stress, σ0.By considering Figure2 and (2) together, it can be seen that the creep equation for the Burgers model is composed of three parts: (3)ε1=σ0E1ε2=σ0E21-exp-E2η2tε3=σ0η1t,where ε1 is the instantaneous elastic strain; ε2 is the recoverable strain decaying with a negative exponent, meaning that it corresponds to the first section of the creep curve; and ε3 is the nonrecoverable strain corresponding to the steady strain rate, meaning that it is related to the second section of the creep curve.


## 2.2. Increment Expression of Rheology

In DDA method, the time,t, is usually divided into a series of time periods. The group of equilibrium equations is set up for each time interval and the stress increment in this period is determined. In addition, the creep deformation of the structural plane is related to the entire stress history. Therefore, in order to calculate the shear creep deformation for time t, it is necessary to establish the recurrence relation for the creep strain increment for each time step. As shown in Figure 3, the time, t, is divided into a series of moments (t0, t1, t2,…,tn-1, and tn), and it is assumed that the stress increments at these moments are Δσ0, Δσ1, Δσ2,…,Δσn-1, and Δσn, respectively. The three adjacent moments tn-1, tn, and tn+1 are selected as the target moments. Their time intervals are Δtn=tn-tn-1 and Δtn+1=tn+1-tn. Therefore, the creep deformations for the three adjacent moments can be given by (4)εctn-1=Δσ0Ctn-1-t0+Δσ1Ctn-1-t1+⋯+Δσn-2Ctn-1-tn-2(5)εctn=Δσ0Ctn-t0+Δσ1Ctn-t1+⋯+Δσn-1Ctn-tn-1(6)εctn+1=Δσ0Ctn+1-t0+Δσ1Ctn+1-t1+⋯+Δσn-1Ctn+1-tn-1+ΔσnCtn+1-tn.Figure 3
Stress increment during each time step of the calculation.If (5) is subtracted from (6), we get (7)Δεtn+1c=εctn+1-εctn=Δσ0Ctn+1-t0-Ctn-t0+Δσ1Ctn+1-t1-Ctn-t1+⋯+Δσn-1Ctn+1-tn-1-Ctn-tn-1+ΔσnCtn+1-tn.From (2), we get (8)Ctn+1-ti-Ctn-ti=1E21-e-E2/η2tn+1-ti-1E21-e-E2/η2tn-ti+1η1tn+1-ti-1η1tn-ti=1E2e-E2/η2tn-ti-e-E2/η2tn+1-ti+1η1Δtn+1=1E2e-E2/η2tn-ti1-e-E2/η2Δtn+1+1η1Δtn+1.By inserting (8) into (7), we get (9)Δεn+1c=1-e-E2/η2Δtn+1ωn+1+νn+1+ΔσnCtn+1-tn,where(10)ωn+1=Δσ0E2e-E2/η2tn-t0+Δσ1E2e-E2/η2tn-t1+⋯+Δσn-1E2e-E2/η2tn-tn-1(11)υn+1=σ0η1Δtn+1+σ1η1Δtn+1+⋯+σn-1η1Δtn+1.Using the same principle, (4) is subtracted from (5): (12)Δεnc=1-e-E2/η2Δtnωn+νn+Δσn-1Ctn-tn-1(13)ωn=Δσ0E2e-E2/η2tn-1-t0+Δσ1E2e-E2/η2tn-1-t1+⋯+Δσn-2E2e-E2/η2tn-1-tn-2(14)υn=σ0η1Δtn+σ1η1Δtn+⋯+σn-2η1Δtn.By comparing (10) and (13), we get (15)ωn+1=ωne-E2/η2Δtn+Δσn-1E2e-E2/η2Δtnω1=Δσ0E2.Similarly, by comparing (11) and (14), we get(16)νn+1=vnΔtn+1Δtn+σn-1η1Δtn+1ν1=σ0η1Δt1.Equations (12)–(16) form a group of recurrence formulas that can be used to determine the rheological deformation increment for each time step of the calculation. In this manner, the rheological deformation at any moment can be determined.


## 3. Contact Rheological Simulation Using DDA

### 3.1. Contact Simulation Using DDA

The focus of the DDA method is a system which consists of blocks with arbitrary shapes. In the DDA method, Genhua put forward a set of strict and efficient contact searching method. Based on this method, all possible contact forms in the system can be determined, including that between two angles, an angle and a line, two lines, and other forms. The contact length between the two angles and that between the angle and the line is zero, while the contact length between two lines is the length of the overlapping part between the two lines. The contact between the two lines can be converted into two contacts. The form of each contact is a line and an angle. The contact length of each contact point is half of the total contact length of the two lines. The contact forms used in DDA are shown in Figure4.Figure 4
Contact simulation in DDA.
(a)
Contact between an angle and a line (b)
Contact between two linesThe DDA method uses the penalty function to ensure that there is no embedding between the contact surfaces of the two blocks. Further, the stiffness of the contact point can be determined using the following equations.Contact between two angles or an angle and a line is(17)PN=k1E,PT=k1E21+μ.Contact between two lines is(18)PN=12k2El,PT=k2El/221+μ.Here,PN is the normal contact stiffness, PT is the tangential contact stiffness, E is the bulk modulus of elasticity, μ is Poisson’s ratio of the block, l is the contact length, and k1 and k2 are the penalty factors (based on the problems simulated, they are usually 10–100).The tangential contact stiffness is usually set based on the relationship between the shear modulus and elastic modulus of the material. In the DDA method, the contact forces can be calculated from the following equation:(19)FN=PNdnFT=PTdt,where FN is the normal contact force, FT is the tangential contact force, and dn, and dt are the normal and tangential deformations, respectively, in the contact position.


### 3.2. Contact Rheological Model for DDA

The rheological deformation corresponding to contact between an angle and a line as well as contact between two lines can only occur in the shear direction. By using (3) and (19), the contact deformation in the direction of the upper-line shear under the action of the shear force, FT, can be obtained: (20)dte=FTPT(21)dtc=FTE21-e-E2/η2t+FTη1t.In (20), dte is the instantaneous elastic shear deformation caused by FT. In (21), dtc is the shear rheological deformation generated after time t by the action of FT. Using the same principle, based on (12)–(16), the rheological deformation of the contact surface along the shear direction can be deduced. That is to say, Δdtnc is substituted for Δεnc and ΔFtn is substituted for Δσn.


### 3.3. DDA Equation That Takes into Account Shear Rheological Deformation

A description of the tangential contact deformation in the DDA method is shown in Figure5. The vertex, P1, of block i is in contact with line P2P3, which is the entrance line. P0 is on line P2P3 and is also the locking point of P1 on line P2P3. That is to say, P0 is the contact point when the shear force between P1 and P2P3 is “0.” When there is a shear force between P1 and P2P3 and it is smaller than the friction force, a tangential spring should be set up between blocks i and j. The spring stiffness would be PT. After the spring has been stretched, the projection point of P1 on P2P3 will be P1′. The distance, dt, between P0 and P1′ satisfies (19).Figure 5
Schematic diagram of the shear deformation and rheology.
(a)
No rheology (b)
RheologyThe shear rheology of the structural plane can be of two kinds:(1) the shear force, FT, remains the same and contact point P1′ exhibits aging displacement along the action direction of FT. The distance between P0 and P1′ increases steadily; this phenomenon is called creep. (2) The position of contact point P1′ remains the same, while the shear force, FT, decreases steadily; this is called stress relaxation. According to these forms, the simulation can be performed in two ways: (1) locking point P0 will exhibit creep displacement along the action direction of FT and will move to P0′. At this time, if the distance between P0′ and P1′ is equal to that of P0P1′, that is, dt′=dt, then it can be used to simulate the creep deformation (the shear force is the same, while the tangential deformation increases). If the distance between P0′ and P1′ is smaller than that of the original P0P1′, that is, dt′=dt-dtc, then the stress relaxation can be simulated (the displacement is the same, while the force is reduced). (2) With an increase in the action time of loading, the rigidity of the spring, PT, which connects P0P1′ is reduced. In this study, the first method was adopted to simulate the shear rheology of the structural plane.It was assumed that whenn-1 calculation time step is over, the locking point is at P0′; the corresponding coordinates will be (x0′,y0′). Further, P0′, P2, and P3 are all points on block j and P1 is the vertex of block i. The coordinates of P1, P2, and P3 are (xk,yk), k=1,2, and 3, respectively. According to (12), the rheological deformation of P1 within the nth time interval and along P2P3 is given as follows: (22)Δdtnc=1-e-E2/η2Δtnωn+νn+ΔFtn-1Ctn-tn-1.It is assumed that the displacements ofP0, P1, P2, and P3 within the time interval are (uk,vk), k=0,1,2,3, respectively. After the rheological deformation has been taken into account, the shear displacement of vertex P1 along the entrance line P2P3 is(23)d=1lP0′P1⃑·P2P3⃑+Δdtnc=S0l+1lx1-x2,y1-y2u1-u0′v1-v0′+Δdnc,where(24)l=x2-x32+y2-y32S0=x3-x2,y3-y2x1-x0′y1-y0′.At this time, the strain energy of the shear spring is(25)∏s=PT2d2=PT2HTDi+GTDj+S0l+Δdtnc2.In this equation,Di and Dj are the six basic deformation components of i and j, respectively: (26)Di=Di1,Di2,Di3,Di4,Di5,Di6Dj=Dj1,Dj2,Dj3,Dj4,Dj5,Dj6.The matricesH and G are given by (27)H=1lTix1,y1Tx3-x2y3-y2G=1lTjx0,y0Tx2-x3y2-y3,where T(i)x1,y1 and T(j)x0,y0 are the displacement functions of P1 and P0, respectively.Therefore, the shear spring matrix for the creep deformation and the equivalent load vector triggered are(28)PTHHT⟶KiiPTHGT⟶KijPTGHT⟶KjiPTGGT⟶Kjj(29)-PTlS0-PTΔdtncH⟶Fi-PTlS0-PTΔdtncG⟶Fj.Equations (28) and (29) are integrated into the tangential stiffness and load array, respectively, of the DDA global equilibrium equation (the other submatrix and the integration methods for the integral stiffness equation are given in the literature [25]). This yields(30)K11K12⋯K1nK21K22⋯K2n⋯⋯⋯⋯Kn1Kn2⋯KnnD1D2⋯Dn=F1F2⋯Fn.By solving (30), the displacement and deformation of each block after the nth calculation time step can be obtained. Thus, the stress of each block and the contact stress between the blocks can also be calculated.


## 3.1. Contact Simulation Using DDA

The focus of the DDA method is a system which consists of blocks with arbitrary shapes. In the DDA method, Genhua put forward a set of strict and efficient contact searching method. Based on this method, all possible contact forms in the system can be determined, including that between two angles, an angle and a line, two lines, and other forms. The contact length between the two angles and that between the angle and the line is zero, while the contact length between two lines is the length of the overlapping part between the two lines. The contact between the two lines can be converted into two contacts. The form of each contact is a line and an angle. The contact length of each contact point is half of the total contact length of the two lines. The contact forms used in DDA are shown in Figure4.Figure 4
Contact simulation in DDA.
(a)
Contact between an angle and a line (b)
Contact between two linesThe DDA method uses the penalty function to ensure that there is no embedding between the contact surfaces of the two blocks. Further, the stiffness of the contact point can be determined using the following equations.Contact between two angles or an angle and a line is(17)PN=k1E,PT=k1E21+μ.Contact between two lines is(18)PN=12k2El,PT=k2El/221+μ.Here,PN is the normal contact stiffness, PT is the tangential contact stiffness, E is the bulk modulus of elasticity, μ is Poisson’s ratio of the block, l is the contact length, and k1 and k2 are the penalty factors (based on the problems simulated, they are usually 10–100).The tangential contact stiffness is usually set based on the relationship between the shear modulus and elastic modulus of the material. In the DDA method, the contact forces can be calculated from the following equation:(19)FN=PNdnFT=PTdt,where FN is the normal contact force, FT is the tangential contact force, and dn, and dt are the normal and tangential deformations, respectively, in the contact position.


## 3.2. Contact Rheological Model for DDA

The rheological deformation corresponding to contact between an angle and a line as well as contact between two lines can only occur in the shear direction. By using (3) and (19), the contact deformation in the direction of the upper-line shear under the action of the shear force, FT, can be obtained: (20)dte=FTPT(21)dtc=FTE21-e-E2/η2t+FTη1t.In (20), dte is the instantaneous elastic shear deformation caused by FT. In (21), dtc is the shear rheological deformation generated after time t by the action of FT. Using the same principle, based on (12)–(16), the rheological deformation of the contact surface along the shear direction can be deduced. That is to say, Δdtnc is substituted for Δεnc and ΔFtn is substituted for Δσn.


## 3.3. DDA Equation That Takes into Account Shear Rheological Deformation

A description of the tangential contact deformation in the DDA method is shown in Figure5. The vertex, P1, of block i is in contact with line P2P3, which is the entrance line. P0 is on line P2P3 and is also the locking point of P1 on line P2P3. That is to say, P0 is the contact point when the shear force between P1 and P2P3 is “0.” When there is a shear force between P1 and P2P3 and it is smaller than the friction force, a tangential spring should be set up between blocks i and j. The spring stiffness would be PT. After the spring has been stretched, the projection point of P1 on P2P3 will be P1′. The distance, dt, between P0 and P1′ satisfies (19).Figure 5
Schematic diagram of the shear deformation and rheology.
(a)
No rheology (b)
RheologyThe shear rheology of the structural plane can be of two kinds:(1) the shear force, FT, remains the same and contact point P1′ exhibits aging displacement along the action direction of FT. The distance between P0 and P1′ increases steadily; this phenomenon is called creep. (2) The position of contact point P1′ remains the same, while the shear force, FT, decreases steadily; this is called stress relaxation. According to these forms, the simulation can be performed in two ways: (1) locking point P0 will exhibit creep displacement along the action direction of FT and will move to P0′. At this time, if the distance between P0′ and P1′ is equal to that of P0P1′, that is, dt′=dt, then it can be used to simulate the creep deformation (the shear force is the same, while the tangential deformation increases). If the distance between P0′ and P1′ is smaller than that of the original P0P1′, that is, dt′=dt-dtc, then the stress relaxation can be simulated (the displacement is the same, while the force is reduced). (2) With an increase in the action time of loading, the rigidity of the spring, PT, which connects P0P1′ is reduced. In this study, the first method was adopted to simulate the shear rheology of the structural plane.It was assumed that whenn-1 calculation time step is over, the locking point is at P0′; the corresponding coordinates will be (x0′,y0′). Further, P0′, P2, and P3 are all points on block j and P1 is the vertex of block i. The coordinates of P1, P2, and P3 are (xk,yk), k=1,2, and 3, respectively. According to (12), the rheological deformation of P1 within the nth time interval and along P2P3 is given as follows: (22)Δdtnc=1-e-E2/η2Δtnωn+νn+ΔFtn-1Ctn-tn-1.It is assumed that the displacements ofP0, P1, P2, and P3 within the time interval are (uk,vk), k=0,1,2,3, respectively. After the rheological deformation has been taken into account, the shear displacement of vertex P1 along the entrance line P2P3 is(23)d=1lP0′P1⃑·P2P3⃑+Δdtnc=S0l+1lx1-x2,y1-y2u1-u0′v1-v0′+Δdnc,where(24)l=x2-x32+y2-y32S0=x3-x2,y3-y2x1-x0′y1-y0′.At this time, the strain energy of the shear spring is(25)∏s=PT2d2=PT2HTDi+GTDj+S0l+Δdtnc2.In this equation,Di and Dj are the six basic deformation components of i and j, respectively: (26)Di=Di1,Di2,Di3,Di4,Di5,Di6Dj=Dj1,Dj2,Dj3,Dj4,Dj5,Dj6.The matricesH and G are given by (27)H=1lTix1,y1Tx3-x2y3-y2G=1lTjx0,y0Tx2-x3y2-y3,where T(i)x1,y1 and T(j)x0,y0 are the displacement functions of P1 and P0, respectively.Therefore, the shear spring matrix for the creep deformation and the equivalent load vector triggered are(28)PTHHT⟶KiiPTHGT⟶KijPTGHT⟶KjiPTGGT⟶Kjj(29)-PTlS0-PTΔdtncH⟶Fi-PTlS0-PTΔdtncG⟶Fj.Equations (28) and (29) are integrated into the tangential stiffness and load array, respectively, of the DDA global equilibrium equation (the other submatrix and the integration methods for the integral stiffness equation are given in the literature [25]). This yields(30)K11K12⋯K1nK21K22⋯K2n⋯⋯⋯⋯Kn1Kn2⋯KnnD1D2⋯Dn=F1F2⋯Fn.By solving (30), the displacement and deformation of each block after the nth calculation time step can be obtained. Thus, the stress of each block and the contact stress between the blocks can also be calculated.


## 4. Verification Using Examples

There has been a previous study on the shear rheological behavior of the weak structural plane on the left bank of Jinping I Hydropower Station [31]. In this study, the time-history curves of the tangential displacement of the structural plane under different stress states were obtained. Then, the rheological parameters were determined by using the Burgers model for fitting the experimental curves. In this section, we used the extended DDA method to numerically validate the results obtained in the above-mentioned study. At the same time, using (20) and (21), the analytical solution for the shear displacement of the structural plane during each test was obtained. Finally, the numerical results obtained using the DDA method and the analytical results were compared with the experimental results obtained in the previous study, in order to validate the proposed method.The model used for the numerical simulations is shown in Figure6. Here, A and B are two rock blocks with the same properties. Along the side of rock block B, a normal constraint is exerted, while the full constraint is exerted at the bottom. Rock block A lies on top of B and is in contact with B. The normal force F1 and the tangential thrust F2 are exerted on rock block A; F2 is smaller than the shear force on the contact surface of rocks A and B. The tangential displacement on the contact surface of rock block A is recorded.Figure 6
Numerical model for testing the shear rheology of the structural plane.The rheological parameters of the left-bank weak structural plane in Jinping I Hydropower Station under different stress states are listed in Table1.Table 1
Rheological parameters of the weak structural plane under different stress states.
Normal stressσ (MPa) Tangential stressτ (MPa) Elastic shear modulusE1 (MPa) Viscoelastic shear modulusE2 (MPa) Viscoelastic coefficientη1 (MPa∗h) Viscoelastic coefficientη2 (MPa∗h) 1.408 0.232 26.87 192.43 34807.72 2141.23 1.156 0.424 40.13 249.33 101760.29 3236.41 0.942 0.439 29.45 260.76 52869.39 3771.63 0.751 0.512 33.65 305.53 23051.02 5736.91A comparison of the DDA numerical results, the analytical results, and the laboratory test results for the shear creep deformation of the structural plane is shown in Figure7.Figure 7
Comparison of experimental, analytical, and numerical results for the shear creep of the structural plane.It can be seen from Figure7 that the shear creep of the structural plane results in different deformation characteristics under the different stress states. When the ratio of the shear stress and normal stress is relatively small, the tangential instantaneous deformation of the structural plane is smaller, and the decay of the shear creep rate is relatively rapid. When the ratio of the shear stress and the normal stress is larger than the normal one, the tangential instantaneous deformation is larger. Further, the shear creep rate is the largest at the beginning, then decays slowly, and finally increases at a smaller but constant rate. The comparison of the numerical, analytical, and experimental results shows that the extended DDA method can accurately simulate the aging deformation of structures under various stress states based on the selection of the appropriate rheological parameters.


## 5. Analysis of Toppling Deformation While Taking Aging Factors into Consideration

The toppling deformation of antitoppling layered rock slopes is a common problem in geotechnical engineering. In order to overcome this issue, researchers have analyzed the underlying mechanism of the deformation and have proposed mechanical models for the phenomenon as well as various factors that affect the phenomenon [32–35]. However, the traditional analysis methods, including the limit equilibrium method for rigid bodies, can only analyze the instantaneous state before toppling and do not take into account the time factor of the toppling deformation. Therefore, they are unable to correctly reflect the slow development of the toppling deformation of the rock mass with the passing of time. Furthermore, they are also unable to predict the long-term deformation characteristics of the slope correctly. Based on the extended DDA method and using the classic Goodman toppling model, the effects of the aging factors on the toppling deformation of rock slopes were studied, while taking into account the shear creep characteristics of the contact interfaces between the blocks.The simulation model is shown in Figure8: the model length is 215 m and height is 125 m and it is divided into 16 blocks that are numbered from 1 to 16 along the foot of the slope to the top. The slope of the slip surface at the bottom is 30°. The basic parameters used in the analysis were as follows: volume weight: 2.5 × 104 N/m3, Poisson’s ratio: 0.25, cohesive force at the interface between the blocks: 0, friction angle of the structural planes between block 1 to block 16: 25°, and friction angle of the slip surface of the slope: φ. First, the values of φ were varied, and, without taking into account the aging factors, the deformation characteristics and failure modes of the slope blocks were observed. All mechanical parameters of the blocks and structural planes used in the paper are listed in Table 2.Table 2
Parameters used in the example.
Block Structural plane Slope angleα (°) Volume weightγ (N/m3) Elasticity modulusE (GPa) Poisson’s ratioν Cohesion forcec (MPa) Friction angle between blocksφ′ (°) Friction angle between block and slopeφ(°) 2.5 × 104 10 0.25 0 25 33~55 30Figure 8
Goodman toppling model.Figure9 shows the final deformation of blocks #1 and #10 for different φ values. Further, Figure 10 shows sketches of the slope deformation for different φ values. It can be seen that the smaller the φ value is, the more obvious the toppling deformation will be and the bigger the displacement of each block will be. With an increase in the φ value, the block displacement reduces gradually, with the block toppling characteristics becoming less obvious. When φ ≥ 50°, the deformation of the slope is very small and is mainly related to elastic deformation. Considering the slope from an instantaneous perspective, it can be said that it is in an absolutely safe condition. In order to further study the effects of the aging factors on the toppling deformation, the shear creep characteristics of the contact interface between the blocks were considered for φ = 50°, and the toppling deformation simulation was performed once again. The rheological parameters used in this simulation are shown in Table 3.Table 3
Rheological parameters of the structural plane.
Structural plane category Cohesive force (MPa) Internal friction angle (°) Elastic shear modulusE1 (MPa) Viscoelastic shear modulusE2 (MPa) Coefficient of viscosityη1 (MPa∗d) Viscoelastic coefficientη2 (MPa∗d) Structural plane along the slope 0 33~55 100 250 200 6.2 ∗ 10 3 Antitoppling structural plane 0 25 100 100 100 1.8 ∗ 10 3Figure 9
Block displacement for different friction angles of the slip surface at the bottom.Figure 10
Deformation form (rheology not considered) for differentφ values.
(a)
φ = 50° (b)
φ = 43° (c)
φ = 38° (d)
φ = 33°Figure11 shows the displacement time-history curves for blocks #1 and #10 for φ = 50°, with the rheology of the structural plane being taken into consideration. Further, Table 4 shows the displacement and deformation modes of the different blocks at different times, with the rheological properties of the structural plane being considered. Figure 12 shows the sketches of the slope deformation at different time for φ = 50°. It can be seen that only a small elastic deformation is generated by the slope for φ = 50° when the aging factors are not considered. Hence, both the safety and stability are high. If the rheological behavior of the contact interface between the layered blocks is considered, it is seen that the deformation of the slope occurs slowly and that the displacement of each block increases. This increases the possibility of instability. Thus, when the safety of the slope taken into account from an instantaneous perspective is satisfied, there is still the possibility of losing stability even after considering the aging factors. Therefore, its long-term safety cannot be guaranteed.Table 4
Displacement and deformation modes of the block system at different times forφ = 50°.
Number of blocks Without considering the rheology Considering the rheology t = 10 days t = 20 days t = 100 days Displacement/m Deformation mode Displacement/m Deformation mode Displacement/m Deformation mode Displacement/m Deformation mode 1 0.17 Sliding 8.12 Sliding 10.82 Sliding 15.15 Sliding 2 0.12 Toppling 8.28 Toppling and sliding 11.04 Toppling and sliding 15.46 Toppling and sliding 3 0.17 Toppling 8.39 Toppling and sliding 11.19 Toppling and sliding 15.67 Toppling and sliding 4 0.21 Toppling 8.51 Toppling and sliding 11.35 Toppling and sliding 15.89 Toppling and sliding 5 0.27 Toppling 8.63 Toppling and sliding 11.51 Toppling and sliding 16.11 Toppling and sliding 6 0.32 Toppling 8.75 Toppling and sliding 11.67 Toppling and sliding 16.34 Toppling and sliding 7 0.37 Toppling 8.91 Toppling and sliding 11.88 Toppling and sliding 16.63 Toppling and sliding 8 0.43 Toppling 8.99 Toppling and sliding 11.99 Toppling and sliding 16.79 Toppling and sliding 9 0.48 Toppling 8.37 Toppling and sliding 11.16 Toppling and sliding 17.02 Toppling and sliding 10 0.57 Toppling 6.02 Toppling and sliding 8.03 Toppling and sliding 17.24 Toppling and sliding 11 0.45 Toppling 5.57 Toppling and sliding 7.43 Toppling and sliding 14.60 Toppling and sliding 12 0.38 Toppling 5.64 Toppling and sliding 7.52 Toppling and sliding 11.93 Toppling and sliding 13 0.30 Toppling 5.00 Toppling and sliding 6.66 Toppling and sliding 9.32 Toppling and sliding 14 0.00 Stability 0.79 Sliding 1.05 Sliding 1.47 Sliding 15 0.00 Stability 0.29 Sliding 0.39 Sliding 0.55 Sliding 16 0.00 Stability 0.05 Sliding 0.07 Sliding 0.10 SlidingFigure 11
Displacement and time-history curve of the blocks forφ = 50°.Figure 12
Aging-related deformation process whenφ = 50°.
(a)
Time = 0 days (b)
Time = 10 days (c)
Time = 20 days (d)
Time = 100 days


## 6. Conclusions

This work studies the rheological properties of rock mass. The viscoelastic rheological model of Burgers is combined with the discontinuous deformation analysis (DDA) method to simulate the shear creep deformation of structural plane of rock mass. A DDA numerical program was also developed and its effectiveness was confirmed based on the curves obtained by the creep test of the structural plane of a rock mass under different stress levels. A further numerical simulation is then carried out to analyze the mechanism responsible for the creep features of the structural plane in the case of the toppling deformation of the rock slope. The conclusions of the study can be summarized as follows:(i)
The main form of slope deformation is creep. That is to say, when the slope is triggered by various factors, its deformation increases with the passing of time. The creep deformation of rock slopes and, in particular, that of hard rock slopes, occurs mainly in the structural plane. Further, numerical simulations of the creep deformation of rock slopes mainly focus on the shear creep on the structure plane.(ii)
The DDA method proposed by Shi Genhua mainly focuses on the block system cut by faults, joints, fractures, and other structural planes. The relative motion between the blocks, the opening and closing of the contact surface, and the sliding process could be simulated accurately. Moreover, the deformation of the block itself can also be simulated. The DDA method is a powerful tool for studying the deformation of rock slopes and whether they will lose stability.(iii)
In this study, the Burgers viscoelastic rheological model was combined with the DDA method; this will extend the DDA method so that it can be used to simulate the shear creep deformation of the structural plane. The extended DDA method was verified by comparing the results obtained with laboratory-based test results; it was found that the two were consistent.(iv)
The extended DDA method was used to simulate the Goodman toppling model. The results showed that the toppling deformation of the slope was composed of two parts: an instantaneous deformation and a creep deformation. The instantaneous deformation increased with a decrease in the friction angle. Forφ = 50°, the maximum instantaneous deformation was 2.17 m. Further, the slope deformation was not significant, and the entire structure was stable. When considering the creep characteristics of the structural plane, the slope deformation increased with the passing of time. The typical characteristics of the toppling deformation were observed after 100 days.(v)
The Burgers viscoelastic model was adopted to describe the creep characteristics of the shear structure in this study. In the same manner, the Maxwell model, the Kelvin model, Nishihara’s model, and other rheological models can also be integrated with the existing DDA method. When performing analyses, based on specific nature of the rock mass, the appropriate rheological model should be used for the simulation.



---

*Source: 1582825-2017-02-15.xml*