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.
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.
The constitutive equation of the Burgers model is as follows:
In this equation, , , , and are the elastic moduli and viscosity coefficients of the elements shown in Figure 1, while is the normal stress and is the total strain.
The creep equation for the Burgers model is given by
In this equation, is the degree of creep and is the time of operation of the stress, .
By considering Figure 2 and (2) together, it can be seen that the creep equation for the Burgers model is composed of three parts: where is the instantaneous elastic strain; is the recoverable strain decaying with a negative exponent, meaning that it corresponds to the first section of the creep curve; and 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, , 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 , it is necessary to establish the recurrence relation for the creep strain increment for each time step. As shown in Figure 3, the time, , is divided into a series of moments (, , , and ), and it is assumed that the stress increments at these moments are , , , and , respectively. The three adjacent moments , , and are selected as the target moments. Their time intervals are and . Therefore, the creep deformations for the three adjacent moments can be given by
If (5) is subtracted from (6), we get
From (2), we get
By inserting (8) into (7), we get where
Using the same principle, (4) is subtracted from (5):
By comparing (10) and (13), we get
Similarly, by comparing (11) and (14), we get
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 Figure 4.
(a) Contact between an angle and a line
(b) Contact between two lines
The 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
Contact between two lines is
Here, is the normal contact stiffness, is the tangential contact stiffness, is the bulk modulus of elasticity, is Poisson’s ratio of the block, is the contact length, and and 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: where is the normal contact force, is the tangential contact force, and , and 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, , can be obtained:
In (20), is the instantaneous elastic shear deformation caused by . In (21), is the shear rheological deformation generated after time by the action of . 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, is substituted for and is substituted for .
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 Figure 5. The vertex, , of block is in contact with line , which is the entrance line. is on line and is also the locking point of on line . That is to say, is the contact point when the shear force between and is “0.” When there is a shear force between and and it is smaller than the friction force, a tangential spring should be set up between blocks and . The spring stiffness would be . After the spring has been stretched, the projection point of on will be . The distance, , between and satisfies (19).
(a) No rheology
(b) Rheology
The shear rheology of the structural plane can be of two kinds: the shear force, , remains the same and contact point exhibits aging displacement along the action direction of . The distance between and increases steadily; this phenomenon is called creep. The position of contact point remains the same, while the shear force, , decreases steadily; this is called stress relaxation. According to these forms, the simulation can be performed in two ways: locking point will exhibit creep displacement along the action direction of and will move to . At this time, if the distance between and is equal to that of , that is, , 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 and is smaller than that of the original , that is, , then the stress relaxation can be simulated (the displacement is the same, while the force is reduced). With an increase in the action time of loading, the rigidity of the spring, , which connects is reduced. In this study, the first method was adopted to simulate the shear rheology of the structural plane.
It was assumed that when calculation time step is over, the locking point is at ; the corresponding coordinates will be (). Further, , , and are all points on block and is the vertex of block . The coordinates of , , and are (), , and 3, respectively. According to (12), the rheological deformation of within the th time interval and along is given as follows:
It is assumed that the displacements of , , , and within the time interval are , , respectively. After the rheological deformation has been taken into account, the shear displacement of vertex along the entrance line iswhere
At this time, the strain energy of the shear spring is
In this equation, and are the six basic deformation components of and , respectively:
The matrices and are given by where and are the displacement functions of and , respectively.
Therefore, the shear spring matrix for the creep deformation and the equivalent load vector triggered are
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
By solving (30), the displacement and deformation of each block after the th 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 Figure 6. 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 and the tangential thrust are exerted on rock block A; 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.
The rheological parameters of the left-bank weak structural plane in Jinping I Hydropower Station under different stress states are listed in Table 1.
A 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 Figure 7.
It can be seen from Figure 7 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 Figure 8: 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.
Figure 9 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.
(a) φ = 50°
(b) φ = 43°
(c) φ = 38°
(d) φ = 33°
Figure 11 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.
(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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was financially sponsored by the National Natural Science Foundation of China (no. 51579252 and no. 51439005), the National Basic Research Program of China (no. 2016YFB020100), and the Independent Research Subject of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (no. 2015TS05).