In this paper, the 3D continuum-discrete coupling analysis method combined with the strength reduction theory is used to study the macrodeformation, stress law, and microfracture mechanism of a high and steep cutting slope of Zhongkai Expressway in Guangdong Province. The 3D continuum-discrete coupling slope model is established based on the finite difference software FLAC and the discrete element program PFC. Then, based on the correlations between stress tensor and force chain, displacement and shear strain increment, plastic deformation, and microfracture, the mechanism of slope instability is analyzed from macro- and microscales. The research shows the following: (1) the continuum-discrete coupling model is reasonable and effective in the three-dimensional slope stability analysis. (2) The coupling domain was penetrated by the shear band, and the upper part displacement of the coupling domain was large, whose direction was basically consistent with the shear band direction; thus, the slope slip surface was formed. (3) The microfracture can represent the macroscopic failure phenomenon, and the slope landslide can be interpreted as the result of shear fracture in slip direction, accompanied by the extension and penetration of tension crack.

1. Introduction

In geotechnical engineering, any change in conditions such as seepage conditions, excavation and unloading, and temperature changes can cause changes in the soil state, resulting in soil deformation and even instability [15]. Especially, during the excavation of high-steep road cutting slope, the stresses in the slope will be redistributed to reach a new equilibrium. Meanwhile, the fractures within the slope are constantly breeding, evolving, expanding, and penetrating each other during the excavation, and thus, the stability of slope decreases gradually [3], and it is important to improve the stability of soils and reduce deformations via reinforcement of soils [6]. However, the slope failure process in practical engineering is difficult to be obtained by physical test. In recent years, with the rapid development of computer technology, the CEM (continuous element method), such as FEM (finite element method) or FDM (finite difference method), is used for the numerical simulation of geotechnical engineering stability evaluation [710]. The CEM as a common macroscopic analysis method can reproduce the process of stress and displacement evolutions during slope failure, but the evolution of microstructures in soil assembly cannot be directly simulated due to the phenomenality of constitutive model and sensitivity of element grid [1113]. In fact, the macroscopic instability failure of rock mass is the result of the development of bonding particle fracture, and the micromechanical analysis of rock deformation properties cannot be carried out only from the traditional CEM. Considering this limitation, Felippa and Park [14] proposed the continuous-discrete coupling method, which can better analyze the evolutions of forces and displacements of soil assembly.

At present, the continuous-discrete coupling method is widely applied in geotechnical engineering. Zhou et al. [15] calculated the interaction between particles at the continuous-discrete interface using DEM (discrete element method). In the continuous element, the velocity at contact interpolated to the node velocity at the interface according to the type function, and the contact force of interparticles is distributed to the interface node according to the weight function. The two-dimensional continuous-discrete coupling analysis was carried out by the Socket I/O interface between PFC and FLAC, and the feasibility of this method was further verified according to the agreement between the results of continuous-discrete coupling and those of FLAC. Chen et al. [16, 17] carried out continuous-discrete coupling numerical simulation of geogrid-wrapped vertical reinforced retaining wall on soft soil foundation through the Socket I/O interface between PFC and FLAC. FDM and DEM were adopted for the simulation analysis of soft soil foundation and retaining wall, respectively. The deformation, particle displacement field, and earth pressure distribution of the retaining wall under static load were analyzed and compared with the experimental results of similar scale model. The results of both two models are consistent, indicating that the method is reasonable to be used in the engineering field of reinforced soil retaining wall on soft soil foundation. Hu et al. [18] characterize the NPR anchor cable with continuous element and the surrounding rock and the anchoring material with discrete element. The NPR anchorage system was simulated with continuous-discrete coupling method. Then, the interaction mechanism between NPR anchor cable and surrounding rock under the static stretching was analyzed to reveal the special anchorage mechanism of NPR anchor cable. The numerical results show that this simulation method can predict and improve the properties of anchor cable anchorage system. Yan et al. [19] adopted the continuous-discrete coupling method combined with the strength reduction theory to realize the transmission of all kinds of data on the coupling boundary. The results show that the data of continuous domain coincide with that of discrete domain, and the microscopic fracture mechanism of coupled domain is consistent with the macroscopic plastic deformation of slope. In addition, the investigation proposed that the instability of slope is mainly caused by shear failure, and the variation of shear strain increment is mainly controlled by horizontal displacement of soil, while the particle displacement mainly results from the variation of contact force chain direction. Ma et al. [20, 21] analyze the influence of fracture dislocation movement on water conveyance tunnels of large hydraulic projects in western strong earthquake region using three-dimensional continuous-discrete coupling method. The fracture process of brittle materials such as rock is controlled by strain localization. And the continuous-discrete coupling method can simulate the formation process of local damage zone and internal fractures [22]. Zhang et al. [23] analyzed the complex fracture phenomenon of concrete under large strain rate based on the continuous-discrete coupling method combined with CT scanning imaging technology. Wang et al. [24] applied the continuous-discrete method to high-speed railway and studied the dynamic contact characteristics between track structure and ballast.

At present, the effectiveness of continuous-discrete coupling method has been proved in many investigations, and it is mostly used for the stability analysis of the interaction between foundation structure and soil [25], while the application of slope analysis is limited to two-dimensional condition, which cannot reflect the real stress state of three-dimensional slope. In this study, a coupled model for slope stability analysis is established based on the continuous-discrete coupling theory, which can be used for numerical simulation of slope with large size. This method not only takes the advantages of modeling and calculation in continuum method but also obtains the microscopic analysis indexes for local domain by discrete element method. This investigation expands the application of 3D continuous-discrete coupling analysis method and provides a new approach for slope stability analysis and mesoscopic failure mechanism.

2. Continuous-Discrete Coupling Method

FDM (finite difference method) is suitable for solving large deformation and nonlinear problems, which assumes the material as continuous medium and uses constitutive relation to calculate and solve. It is widely used in geotechnical engineering due to the high computational efficiency and strong practicability. However, a single macroscale continuous element cannot directly reflect the law of force and motion interparticles in soil, while the discrete element method can directly reproduce the mechanism of deformation development [26]. In actual engineering, the calculation for a site slope model is too large due to the millions of particles for DEM simulation, which requires high configuration of the computing machine. Besides, there is no need to analyze the whole slope model using DEM, but the region of interest can be analyzed with DEM, while other regions use continuous element method. In this sense, the continuous-discrete coupling method can not only obtain the macrodeformation failure of slope but also reveal the microfailure mechanism of slope in essence. FLAC3D 6.0 shows good compatibility with PFC3D 6.0 and does not need to be coupled through the Socket I/O interface in the traditional coupling method. Therefore, the continuous-discrete coupling method of version 6.0 is adopted for slope rupture mechanism analysis in this study.

2.1. Coupling Method

FLAC3D is a finite difference program whose forces act on the nodes and are computed by the constitutive equations, i.e., the force-displacement law. PFC3D is a discrete particle flow software, and the forces acting on the particles can be computed by the Newtonian equations of motion. The coupling schematic diagram of FLAC3D and PFC3D is shown in Figure 1. The FLAC3D-PFC3D coupling includes the mechanical and spatial variables, while the walls are the key to successful coupling.

2.2. Coupling Equation

The PFC walls overlap with the FLAC zone or the surface of the shell elements, transferring the forces and moments acting on the walls from the PFC to the FLAC through an equivalent mechanical system. The forces and moments acting on the wall by the PFC3D elements are successfully transferred into the FLAC3D solid cell by an equivalent mechanical system. Then, the coupling is realized successfully based on the spatial information exchange; i.e., the wall provides the spatial location of particles in PFC3D, and the deformation of FLAC3D is reflected to the zone surfaces (the wall of PFC3D). The coupling solution will be described by gravity interpolation (extrapolation) in Figure 2.

As shown in Figure 2, the wall consists of triangular faces connected by edges, where the vertex velocity and position can be specified as a function of time. According to the coupling method, the system of equivalent forces at the face vertices is determined by obtaining the contact forces and moments of wall, which are transferred to the grid points along with the stiffness.

Note that the vector pointing from the CP points to the vertices of each triangle is , where , ; the force applied at each grid point or node is ; the contact force applied at contact point is ; and the contact moment due to bonding is . Then, the total moment acting on the small plane is

According to the equivalence of forces, we have

The unit vector pointing in the normal to the triangle is denoted as , and the force vector along the surface of triangle is and the tangential vector is

Here, we consider a local coordinate system , where is consistent with the direction and is consistent with the direction . Because is located on the surface of the triangle, the -component of each is zero, i.e., . Then, the -component of the force at vertex in local system is

And the equilibrium equations in the - and -directions can be written as

Due to the weighting term, the -direction distribution of the force is . Thus, equations (7)–(9) can be reduced to where is the weighting factor.

In the absence of additional constraints, these underdetermined equations have infinitely many sets of solutions. Thus, in this paper, the additional constraint is given by

The forces obtained from the above calculation process are converted to forces in the overall coordinate system and then applied to the appropriate grid points or structural element nodes.

3. The Establishment of Coupling Model

Generally, a continuous-discrete coupling model includes a continuous model and a discrete model. The construction of the coupling model includes the establishment of initial geometric model, the selection of constitutive model and contact model, the assignment of macroscopic and mesoscopic parameters, the boundary condition, and the loading method.

3.1. Initial Geometric Model

The initial model includes a continuous element geometry model and a discrete element geometry model, and the continuous element geometry model is constructed according to the original dimensions of the slope in the engineering. The choice of the discrete domain is important, as both the location and size need to be considered. In order to observe clearly the microscopic phenomena such as contact forces and fractures in the discrete domain and thereby analyze the microscopic mechanism in the process of slope instability, the coupling domain should be selected in the slope failure area. In other words, the discrete domain should be traversed by the slope slip zone [27]. At the same time, the discrete domain has to be at the place where the slope displacement is large. So the discrete domain is placed under the fourth step of the slope. On the other hand, the size of discrete domain should not be too large; otherwise, it will affect the speed of the calculation. Moreover, as an embedded body, the discrete domain will have a strong boundary effect on the continuous domain if the size is too large. Therefore, the size of the discrete domain used for the micromechanical analysis is determined as , as shown in Figure 3.

3.2. The Selection of Microconstitutive Model and Microcontact Model

The continuous-discrete coupling model needs to determine the constitutive model of the continuous domain and the contact model of the discrete domain. The constitutive relation of the continuous domain model is the basis of the model calculation, which is used to describe the mechanical characteristic of materials [28, 29]. In this paper, the Mohr-Coulomb model is selected as the constitutive model of continuous domain, which is often used in slope stability analysis. The contact model of discrete model is considered to be parallel bond model, which is widely used in geotechnical engineering [30, 31].

3.3. The Identification of Macroscale and Microscale Parameters

The parameters of the continuous-discrete coupling model include macroscopic parameters in the continuous domain and microscopic parameters in the discrete domain. The macroscopic parameters are selected according to the constitutive model, and the microscale parameters are selected according to the contact model. The values of macroscopic parameters are obtained from engineering geological survey reports and laboratory rock experiments in the study, and the values of macroscopic parameters are shown in Table 1. Besides, microscopic parameters are calibrated based on DEM numerical simulation, and the values of microscopic parameters are shown in Table 2.

The Mohr-Coulomb strength curve of the sample is shown in Figure 4, the cohesion is 59.69 kPa, and the internal friction angle is 35.03°. The curves of axial stress-strain and radial strain-axial strain of discrete element samples under uniaxial compression compared with those of ideal elastoplastic materials are plotted in Figures 5 and 6, respectively. Clearly, the elastic modulus and Poisson’s ratio of the continuous element sample are 1.5 GPa and 0.3, respectively. Moreover, it can be seen from Figures 5 and 6 that the mechanical characteristic of discrete element sample and continuous element sample are basically the same, which indicates that the parameters are reasonable.

3.4. The Initial Conditions

By using the surface constraint method, the initial conditions in the continuous domain are normal displacement constraints() on the surrounding surfaces of the continuous domain model, including left, right, front, and back surfaces, and the bottom surface of the model is fully fixed (). The initial conditions of the discrete domain are zero contact force (force-contact multiply 0.0), zero linear force (), zero displacement (displacement multiply 0.0), free linear velocity, and free rotation velocity of the particles.

3.5. The Loading Method

Based on the reduction of strength parameters, the continuous-discrete coupling model is loaded with gravity stress, and the excavation process is divided into eight steps from top to bottom. The reduction factor corresponds to the safety factor calculated by the finite difference method, which is taken as 1.6. The equilibrium condition of the model is that the local unbalance force ratio less than 10-4, the model reaches equilibrium when both continuous and discrete domains converge, and the calculation terminates at 62450 steps.

4. Results and Analysis

4.1. Consistency of Continuous and Discrete Results

During the continuous-discrete simulations, the forces and displacements of continuum nodes and discrete particles must be consistent when the model reaches equilibrium. Figure 7 shows the stress condition and contact force distribution within the coupling domain.

It can be seen that the stress near the lower right of the discrete domain is larger while smaller in the upper left area, which corresponds to a thicker force chain in the lower part of the discrete domain and a thinner force chain in the upper part. Besides, the direction of the strong chain is practically consistent with the direction of the major principal stress. It can be also seen from Figure 8 that the particle displacement in the discrete domain is consistent with the displacement obtained from continuous domain. This consistency explains the accuracy of the coupling model to a certain extent.

4.2. Correlation Analysis of Shear Strain Increment and Particle Displacement

At macroscale level, the shear band is the significant characterization of the slope slip [3234]. Correspondingly, the analysis of particle displacement field of discrete domain, at microscale level, can improve the understanding of the movement of granular materials and allow better judgment in relation to the trends associated with displacements [35]. Figure 9 shows the development process of the shear band after slope excavation. When the excavation is completed initially (Figure 9(a)), the shear strain increment and displacement are almost zero. As the calculation reaches 5000 steps (Figure 9(b)), the shear strain increment begins to appear at the foot of the slope, and the coupling discrete domain undergoes a slight displacement. Thereafter, obvious shear strain increment extends upward along an arc from the toe of slope, and the displacement of the coupling discrete domain increases obviously and the displacement direction is consistent with the arc direction when calculating 10,000 steps (Figure 9(c)). As the calculation progresses, a large number of shear bands appear in the lower part of the slope, but they have not yet penetrated. At this point, the slope appears locally unstable, but the whole is in a stable state, with shear bands and large displacements in the coupled discrete domain (Figure 9(d)). Finally, when the model calculates the equilibrium (Figure 9(e)), the shear band basically penetrates and traverses the coupling domain. Slip surfaces appear within the slope and large displacements occur in the part of the coupled domain that lies on the slip body.

4.3. Correlation Analysis of Macroscale Plastic Deformation and Microscale Cracks

The penetration of the plastic zone is an important criterion for judging the instability of the slope [36, 37], and the generation and development of internal microcracks are the internal reasons for the gradual penetration of the plastic zone of the slope. Detailly, the excavation of the slope will produce elastic-plastic deformation, accompanied by the occurrence of internal microscale cracks occur. When the cracks develop to a certain extent, the plastic zone penetrates and the slope becomes unstable. Figure 10 shows the development process of the macroscale plastic zone and the microscale cracks within the slope. In discrete regions, green dots represent shear cracks, red dots represent tensile cracks, and red areas in continuous domains represent plastic zones. In the early stage after the excavation of the slope, the plastic zone is mainly distributed in the toe of slope and the deep part of the slope. Meanwhile, there are few cracks in the discrete region, and the crack distribution is relatively uniform. As the step increases, a top-down through plastic region was formed in the slope. The plastic zone crosses the middle and lower part of the coupled domain, in which a large number of cracks appear. Furthermore, in the discrete domain, the discrete domain is dominated by shear damage, supplemented by tensile damage, which is consistent with the dominance of shear damage in slope failure.

4.4. Correlation Analysis of Crack Development and Slope Failure

Crack development is a microscale characterization of the progressive failure of slopes [38, 39]. Figure 11 shows the whole process of crack development during slope instability. When the model computes 1000 time steps, the cracks are mainly concentrated at the bottom of the coupling discrete domain, and a total of 6220 cracks are generated, including 6148 shear cracks and 72 tension cracks. When the calculation reaches 5000 time steps, the bottom cracks further increase and expand to the middle of the coupling discrete domain, resulting in 9308 cracks, including 8855 shear cracks and 453 tension cracks. As the calculation proceeds, there are many cracks in the middle area, with a total of 12,375 cracks, including 11,393 shear cracks and 982 tension cracks at 10000 time steps. After that, the cracks begin to develop upward, and the number of cracks continued to gradually increase. Overall, it can be seen that the cracks increase rapidly from 0 to 10,000 time steps and then increase more slowly and gradually level off, which implies that the area where the discrete domain is located has been failure in the early stage of the slope instability. The crack-time step curve shows that the slope instability is dominated by shear failure, supplemented by tension failure [4042].

Figure 12 gives the polar isodensity map of cracks after slope instability. It can be seen from Figure 12 that the crack inclination is mainly concentrated in the range of 0-30°, which is generally consistent with the direction of the shear band and plastic zone near the coupled domain. Hence, the slip of slope is microscopically the result of the extension and penetration of shear and tension cracks in the shear band.

5. Discussion

Since the change of contact force between particles is an important reason for the macroscale displacement and has a close correlation to slope instability, the contact force and particle displacement in the coupled domain were carefully analyzed and discussed. Figure 13 shows the development of contact force and displacement in the coupled domain during slope instability and failure process. At the early stage of the completion of the slope excavation, the internal force chain is relatively evenly distributed. Under gravity loading, the vertical force chains inside the slope are strong force chains, and the horizontal force chains are less distributed, and the overall contact force is weak and the displacement is small, with the contact force distributed in 20 kN~40 kN, and the displacement is about 0~0.02 m. As the time step increases, the contact force is enhanced and the force chain distribution gradually changes, the upper force chains in the discrete domain are thin, and the lower force chains are thick and dense.

Since the lower part is located in the shear zone where cracks are in concentrated distribution, and for the stress at the tip of the tension cracks is concentrated, the lower part has a local contact force concentration phenomenon, with the magnitude reaching 200 kN. With the formation of the slip surface, the force chains in the upper part of discrete domain are located inside the sliding mass break, and the displacement in the upper part of discrete domain is significantly larger than the lower displacement, with the former magnitude reaching 0.1 m.

6. Conclusions

The force characteristics and deformation rule of high and steep cutting slopes during the instability process are analyzed. The correlations between stress and contact force, strain and particle displacement, shear strain increment and particle displacement, plastic failure and microfracture, crack development and slope instability, and coupled domain force chain and displacement are focused in this paper, and the slope instability mechanism is explored from the view of macro and microscales. The specific findings are summarized below:

In the established 3D continuous-discrete coupled slope model, a good transition of forces and displacements between the continuous domain and the discrete domain was realized, which indicates the reasonable validity of the continuous-discrete coupling in the 3D slope stability analysis.

The penetration of the shear zone is one of the criteria for slope instability. With the increase of equilibrium time step, the increment of shear strain increases, and the down-top penetrated shear band gradually forms. As the shear band crosses the coupling domain, the displacement of the upper part of the coupling domain is larger, the displacement direction is basically the same as the direction of the shear zone, and the slope slip surface is formed.

The slope instability is the result of the expansion and penetration of shear and tension microcracks in the shear zone. The slope failure is dominated by shear failure accompanied by tension cracks. The instability of the slope is a process from local failure to overall instability, and the coupling domain has been destroyed in the early stage of overall instability.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Han Dayong performed the data analyses and wrote the manuscript; Yang Chao established the numerical model and analyzed the results; Wang Liang contributed significantly to data analysis and manuscript preparation; Li Dong reviewed and edited the manuscript; Liu Yang contributed to the conception and methodology of the study. All authors have read and agreed to the published version of the manuscript.


The authors gratefully acknowledge the financial support provided by the PowerChina Roadbridge Group Co., Ltd. (No. LQKY2017-03).