# Failure Mechanism of Horizontal Layered Rock Slope under the Coupling of Earthquake and Weathering

Authors: Yun Tian; Lin-feng Wang; Hong-hua Jin; Yu-chao Xia; Biao Zeng
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

In order to further reveal the dynamic failure process of horizontal complex layered rock slope in the coupling action of earthquake and weathering, based on the three-dimensional particle flow software PFC3D, the rock triaxial test was first carried out to obtain the macroscopic rock mechanics parameters. Secondly, PFC3D was used to establish a three-dimensional discrete element model of slope under 4 modes. Finally, the Wenchuan earthquake wave dynamic input is used to simulate the coupling effect of each model slope on earthquake-weathering. Under the progressive failure process, the results show that, during the triaxial test, the mudstone particles and sandstone particles exhibited shear failure and tensile failure, respectively. The seismic failure mode of the horizontal complex layered rock slope model is divided into tensile fracture-horizontal sliding failure (cavity depth 1.5 m and fissure depth 25%), tensile cracking-caving-dumping failure (cavity depth 3 m and fissure depth 37.5%), caving-dumping-rotation failure (cavity depth 4.5 m and fissure depth 50%), and caving-slip-rotation failure (cavity depth 6 m and fissure depth 62.5%), and the maximum vertical displacement of Model 2 can reach 6.2 m, and the maximum positive horizontal displacement is 7.8 m. With the deepening of the weathering degree, the coordination number of the rock block shows the attenuation advance effect and the elevation attenuation effect during the seismic period of 2–15 s. The total energy dissipation of the slope can be divided into two stages: fracture collision and overall collapse. The energy dissipation is more concentrated in the rock fracture collision stage, accounting for 89% to 95% of the total energy.


---

## Body

## 1. Introduction

The interior of the slope is often unstable, understable, or limit-balanced due to multiple groups of rock horizontal structural planes. With the increase of the weathering depth of the rock cavity and the expansion of the cracks in the main control structure, the possibility of slope instability will increase; this kind of slope is very easy to cause instability and collapse under the action of frequent earthquakes, and it is very easy to cause instability, collapse, and destruction and has the characteristics of diversity, complexity, serious danger, great harm, and rapid disaster. The instability of this slope will cause huge economic and property losses in the construction of local infrastructure [1, 2].Taking the horizontal layered rock slope as an example, this type of slope failure is due to the continuous instability and fracture of the internal structural surface of the slope under load, resulting in the deterioration of the stability of each rock mass of the slope until collapse. In recent years, many scholars have established such slope stability calculation methods based on fracture mechanics, catastrophe theory, and energy, etc.; for example, Chen et al. [3, 4] consider different load conditions, aiming at sliding, dumping, and falling rock quality and slopes, combined with rock mass structure theory and limit equilibrium theory, established the stability calculation methods of three typical types of rock slopes, and established the calculation methods of fracture strength factors of the first and second types based on fracture mechanics; Tang [5] based on the mechanics mechanism of rock mass interaction established the equivalent bending moment calculation formula of each rock mass from the rotation angle and calculated the rock collapse according to the creep damage, fracture criterion, and soft rock weathering speed, based on which the sequence of rock collapse was obtained; Wang et al. [6] took complex slow-dip rock slopes as the research object and established interlayer loads and stability of complex slow-dip rock slopes based on fracture mechanics and calculation method of material mechanical stability coefficient. Aiming at the dynamic failure process of complex horizontal rock slopes, e.g., He et al. [7], the discrete element software PFC2D was used to establish a rock mass slope model with dense joints. The mechanical properties of the complete rock and joints were studied, and the joints were discussed. Salmi et al. studied the effect of connectivity on the failure mode of slope [8]. Based on empirical method and numerical method, the stability of practical slope is evaluated and studied jointly; Tang et al. [9] used 2D discrete element simulation to study the genetic mechanism of Caoling Landslide in 1941; Hu et al. [10] based on the two-dimensional particle flow software (PFC2D) artificially synthesized rock mass technology and revealed the horizontal and intermittent joint thickness under different combinations of rock bridge length and joint spacing. The failure modes, dynamic response laws, and stress evolution characteristics of layered rock slopes under earthquake action have found that the failure mode of thick layered rock slopes with horizontal intermittent joints under earthquake action can be classified as collapse type failure, split-slip-block dumping mixed failure, and split-horizontal slip mixed failure, and horizontal intermittent joints are the key factors controlling the dynamic stability of the slope; Bian et al. [11] based on the two-dimensional particle flow software PFC2D synthetic rock technology (SRM) reveals the failure mode and dynamic response law of rock slopes with intermittent joints and discontinuities with different combinations of rock bridge dip angle and joint spacing. It is found that the potential sliding surface composed of bedding intermittent joints and rock bridges is the key factor to control the dynamic stability of the slope. In order to promote the reliability analysis of practical engineering problems, Ji et al. [12] proposed a simplified forward and reverse first-order reliability iterative algorithm (FORM) based on reliability design (RBD). At the same time, based on Newmark’s slider theory, Ji et al. [13] also studied the seismic slope failure mechanism of rotating sliding mass and obtained a simplified relationship between the rotating motion and the horizontal motion of circular failure mass.Most of the above scholars have focused on rock slopes and analyzed the stability and failure mechanism of slopes from the perspective of mechanical theory and PFC2D. The existing research mostly stays in the two-dimensional slope model, but does not consider the stability evolution process and the instability failure process of this kind of slope under the earthquake action. In addition, the interaction between the longitudinal particles of the three-dimensional slope also directly affects the failure mechanism of the slope, so the research on the dynamic failure process of the slope is still in the initial stage under the coupling of earthquake and weathering. It is still in its infancy, and the research on the mechanism of seismic failure of horizontally complex layered rock slopes has gradually become one of the research hotspots and technical difficulties in this field. Therefore, this paper uses the particle flow PFC3D method to establish horizontally complex layered rock slopes with different degrees of weathering, and the dynamic instability and failure process of such slopes under the coupling of earthquake-weathering are simulated by inputting the Wenchuan seismic wave. The damage degree of each rock block of the slope under different weathering degrees of the rock cavity was explored by the earthquake, and the caving time of the rock blocks was obtained, which revealed the mechanism of such damage under the coupling effect of earthquake-cavity weathering. It has important scientific guiding significance and practical value for the evaluation and prediction of such slope stability.In this paper, a horizontal complex layered rock slope of Duwen Highway in Sichuan is taken as the research object. The slope is about 22.5 m high. The upper rock mass is composed of feldspar sandstone, with good integrity, and is about 16.5 m high and has a natural bulk density of 24.7 kN/m3, the fracture toughness is 26.0 MPa·m1/2; the lower rock mass is mudstone, with poor integrity, and is about 6.0 m high and has a natural bulk density of 25.3 kN/m3. According to the dangerous rock chain rule [3–6] and the structural surface layout on site, this paper sets up two macrochains 1and 2 from the slope to the interior direction of the slope body, of which chains 1 and 2 are about 3 m wide. The dangerous rock bodies in each macrochain are numbered from bottom to top (chain 1 is 1-1, 2-1, 3-1, and 4-1, and chain 2 is 1-2, 2-2, 3-2, and 4-2), each macrochain contains 4 microchains (see Figure 1). The height of the first layer of the slope is 3.5 m, the height of the second and third layers is 4 m, the height of the fourth layer is 4.5 m, and the height of 4‐2 is 5 m.Figure 1
The real scene of the slope. (a) Duwen Highway Slope. (b) Rock slope model.
(a)(b)The above rock slope body is composed of sand and mudstone. In order to obtain the mesoparameters of the slope body in the process of PFC3D numerical simulation of particle flow, this article conducts a mesoscale observation through two methods: indoor rock triaxial test and PFC numerical triaxial test. For parameter debugging, the common trial adjustment method uses indoor triaxial-numerical triaxial stress-strain curve comparative analysis [14, 15]. The disadvantage of this method is that the content of the comparison is single, and the mesoparameters obtained are not reasonable. This paper intends to conduct a comparative analysis through the indoor triaxial-numerical triaxial stress-strain curve, the degree of compressive strength attenuation, and the failure process in order to obtain more reasonable rock microscopic parameters.


## 2. Materials and Methods

### 2.1. Indoor Triaxial Test

For the laboratory mechanical test of rock, many scholars have influenced the rock strength parameters and considered a lot of external and internal factors. Based on the research results of the literature studies [16–20], the triaxial test is adopted and carried out on the RMT-150C rock mechanics test system. The specific parameters of the test device are shown in Table 1. The rock specimens are taken from the rock slope at the above Duwen Expressway and are transported back to the laboratory by on-site sealed packaging. In order to ensure the moisture content of sandstone and mudstone, the sample preparation process is dry and wet processing, after several processes of cutting and grinding. To ensure the accuracy requirements of the two rock specimens, the specimen is a cylinder with a diameter of 50 mm and a height of 100 mm, and the error range is controlled within ±2 mm. The mudstone specimen number A and the sandstone specimen number B are shown in Figure 2.Table 1
Parameters of RMT-150C rock press.
Maximum load (kN)FrequencyMaximum confining pressure (MPa)Confining pressure rateMeasurement range of displacementStrain measurement range10000.2∼2 Hz1000.001∼1 MPa/s±50 mm±0.01Figure 2
Real scene of rock specimen and test device. (a) Rock sample. (b) Laboratory apparatus.
(a)(b)


### 2.2. PFC Numerical Triaxial Test

For the simulation study of geotechnical particle flow, the key condition to determine its rationality and reliability lies in the selection of mesoparameters in the model. In this paper, by constantly changing the mesoparameter settings, the triaxial test is repeated in large numbers and determined by the trial-and-error method and it constantly calibrates the mesoparameters. Combined with the principle of PFC3D triaxial test servo, a numerical model of rock specimens was established (see Figure3), and microscopic damage analysis and mesoparameter calibration of rock were carried out. The model generates particles in a cylindrical wall with a height of 10 m and a diameter of 5 m according to a real triaxial test. Considering the effect of particle size on the test results [14], the particles are appropriately enlarged during the generation process, which improves the calculation efficiency. Without affecting the accuracy of the simulation results, the particle radius is taken from 0.2 to 0.25 m. The confining pressure is correspondingly loaded according to the data of 2 MPa, 4 MPa, 6 MPa, and 8 MPa in the indoor test. The upper and lower loading plates move toward each other at a speed of 10 mm/min until the specimen fails, and the stress-strain curve of the entire process is recorded.Figure 3
Schematic diagram of PFC triaxial test.


### 2.3. PFC Model Method

#### 2.3.1. Geological and Numerical Model

Based on the calculation formula of rock slope stability proposed by Chen et al [21, 22], according to the degree of rock cavity weathering and the degree of crack penetration of the main structural plane, this kind of slope is divided into 4 types. The crack penetration of particle flow slope and rock cavity weathering are shown in Figures 4(a) and 4(b), respectively. Table 2 shows the weathering of each model cavity and fracture.Figure 4
Layout of slope fissure and weathering of rock cavity (unit: m). (a) Fracture (DFN joint) and force chain grid. (b) Weathering layout of particle flow rock cavity.
(a)(b)Table 2
Weathering situation of each model cavity and fissure.
ModelCrack depth (m)Fracture penetration rate (%)Cavity thickness (m)M14.125254.5M26.237.53M38.25501.5M410.3162.50


#### 2.3.2. Input of Seismic Waves

In the process of seismic failure of rock slopes, the propagation speeds of shear waves and longitudinal waves are successive and have different effects on the stability of rock slopes with different degrees of weathering. Often only the influence of the horizontal seismic force of the transverse wave is considered [23]. Later, some scholars have shown that the coupling effect of the transverse wave and the longitudinal wave has an amplified effect on the slope stability during the influence of the seismic force on the stability of the rock slope, and the probability of the two acting at the same point at the same time is small, and the actual effect is not large [24], and according to the “Code for Seismic Engineering of Highway Engineering” [25], it is known that the slope in high-intensity areas is significantly affected by the vertical wave vertical seismic action. Therefore, this paper separates the seismic longitudinal and longitudinal waves separately, based on the rock slope of Duwen Highway in Sichuan and based on the calculation results of the slope stability coefficients of different weathering degrees under the seismic longitudinal and longitudinal waves [22].The calculation formulas are shown in formulas (1)∼(8), the calculation parameters are shown in Table 3, and the calculation result is shown in Table 4. The calculation formula of the slope stability coefficient under the shear waves of Models 1∼4 is as follows:(1)Fs1′=mgla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(2)Fs2′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(3)Fs3′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mgla+maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(4)Fs4′=mgcosβ−maxtsinβ−1/18γwe2tanφ+cH−e/sinβmgsinβ+maxtcosβ.Table 3
Calculation parameters.
m (kg)ax (m/s2)ay (m/s2)γw (kN/m3)flk (MPa)fok (MPa)c (MPa)Φ (°)β (°)2.45 × 1059.589.481012.51.120.743890Table 4
Calculation result of the stability factor.
ModelDirectionModel 1Model2Model 3Model 4Shear wave action3.322.151.151.80Longitudinal wave action560.96178.823.600.92The calculation formula of the slope stability coefficient under Model 1∼4 earthquake longitudinal wave is as follows:(5)Fs1″=mg+maytla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(6)Fs2″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(7)Fs3″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mg+maytla+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(8)Fs4″=mgcosβ−maytcosβ−1/18γwe2tanφ+cH−e/sinβmg+maytsinβ,where β(°) is the vertical structural plane dip angle; c is the rock cohesion (MPa); φ is the internal friction angle (°); Gw is the rock’s own weight (kN); Px is the horizontal seismic force (kN); Py is the vertical seismic force (kN); Q is the dangerous rock fissure water pressure (kN); la is the horizontal distance between the center of gravity of the rock block and the weathering point of the cavity; lb is the horizontal distance between the vertical structural surface and the rock cavity surface; e is the fracture penetration length; H is the height of the dangerous rock; h0 is the vertical distance from the center of gravity to the top surface of the base; γw is the bulk density of the fracture water (kN/m3); m is the mass of the rock (kg); ax(t) is the horizontal seismic acceleration of transverse wave (m/s2); flk (MPa) is the standard value of the tensile strength of sandstone, fok (MPa) is the standard value of the tensile strength of mudstone; ay(t) is the vertical seismic acceleration of longitudinal wave (m/s2).It can be seen that the stability coefficient of Model 1 and Model 2 under the action of transverse wave is much smaller than the result under the action of longitudinal wave. Therefore, the seismic shear wave is the main factor of seismic damage of these two slope models. Similarly, the seismic wave with the smaller stability factor is selected as the main factor, that is, the main factor affecting the damage of Model 3. For the shear wave, the main factor affecting the destruction of Model 4 is the longitudinal wave. This result further improves the calculation efficiency.In this paper, the E-W and U-D waves obtained at the Wolongtai station in the 2008 Wenchuan earthquake in Sichuan were used as the transverse and longitudinal waves, and the first 25 s of the seismic time history records was taken as the horizontal and vertical seismic wave input models. The curve is shown below. Model 1∼3 seismic wave action uses shear wave, and Model 4 seismic wave action uses longitudinal wave (see Figure5).Figure 5
Model input seismic wave acceleration time history curve. (a) E-W wave curve (M1, M2, and M3). (b) U-D wave curve (M4) (the dataset is provided by China Earthquake Networks Center, National Earthquake Data Center (http://data.earthquake.cn)).
(a)(b)


#### 2.3.3. Boundary Conditions and Monitoring Plan

When analyzing the dynamic response of a slope under earthquake action, you need to enter the seismic wave acceleration time history curve at the bottom of the model. In order to prevent the reflection of seismic waves at the model boundary, a viscous boundary needs to be used. The normal and tangential boundaries of the boundary are provided with dampers. During the loading process, the seismic acceleration is converted into boundary acceleration, which drives the model particles to vibrate.This article is to monitor the law of rock slope movement damage under the action of earthquake. In the PFC data monitoring and analysis, the history command is used to track the data of the object under study. The tracked content includes data such as displacement, coordination number, and energy. Each rock block is grouped, and the movement data of each rock block are monitored in groups. A total of 8 groups are arranged. The boundary conditions and monitoring scheme are shown in Figure6.Figure 6
Boundary conditions of slope dynamic calculation.


## 2.1. Indoor Triaxial Test

For the laboratory mechanical test of rock, many scholars have influenced the rock strength parameters and considered a lot of external and internal factors. Based on the research results of the literature studies [16–20], the triaxial test is adopted and carried out on the RMT-150C rock mechanics test system. The specific parameters of the test device are shown in Table 1. The rock specimens are taken from the rock slope at the above Duwen Expressway and are transported back to the laboratory by on-site sealed packaging. In order to ensure the moisture content of sandstone and mudstone, the sample preparation process is dry and wet processing, after several processes of cutting and grinding. To ensure the accuracy requirements of the two rock specimens, the specimen is a cylinder with a diameter of 50 mm and a height of 100 mm, and the error range is controlled within ±2 mm. The mudstone specimen number A and the sandstone specimen number B are shown in Figure 2.Table 1
Parameters of RMT-150C rock press.
Maximum load (kN)FrequencyMaximum confining pressure (MPa)Confining pressure rateMeasurement range of displacementStrain measurement range10000.2∼2 Hz1000.001∼1 MPa/s±50 mm±0.01Figure 2
Real scene of rock specimen and test device. (a) Rock sample. (b) Laboratory apparatus.
(a)(b)


## 2.2. PFC Numerical Triaxial Test

For the simulation study of geotechnical particle flow, the key condition to determine its rationality and reliability lies in the selection of mesoparameters in the model. In this paper, by constantly changing the mesoparameter settings, the triaxial test is repeated in large numbers and determined by the trial-and-error method and it constantly calibrates the mesoparameters. Combined with the principle of PFC3D triaxial test servo, a numerical model of rock specimens was established (see Figure3), and microscopic damage analysis and mesoparameter calibration of rock were carried out. The model generates particles in a cylindrical wall with a height of 10 m and a diameter of 5 m according to a real triaxial test. Considering the effect of particle size on the test results [14], the particles are appropriately enlarged during the generation process, which improves the calculation efficiency. Without affecting the accuracy of the simulation results, the particle radius is taken from 0.2 to 0.25 m. The confining pressure is correspondingly loaded according to the data of 2 MPa, 4 MPa, 6 MPa, and 8 MPa in the indoor test. The upper and lower loading plates move toward each other at a speed of 10 mm/min until the specimen fails, and the stress-strain curve of the entire process is recorded.Figure 3
Schematic diagram of PFC triaxial test.


## 2.3. PFC Model Method

### 2.3.1. Geological and Numerical Model

Based on the calculation formula of rock slope stability proposed by Chen et al [21, 22], according to the degree of rock cavity weathering and the degree of crack penetration of the main structural plane, this kind of slope is divided into 4 types. The crack penetration of particle flow slope and rock cavity weathering are shown in Figures 4(a) and 4(b), respectively. Table 2 shows the weathering of each model cavity and fracture.Figure 4
Layout of slope fissure and weathering of rock cavity (unit: m). (a) Fracture (DFN joint) and force chain grid. (b) Weathering layout of particle flow rock cavity.
(a)(b)Table 2
Weathering situation of each model cavity and fissure.
ModelCrack depth (m)Fracture penetration rate (%)Cavity thickness (m)M14.125254.5M26.237.53M38.25501.5M410.3162.50


### 2.3.2. Input of Seismic Waves

In the process of seismic failure of rock slopes, the propagation speeds of shear waves and longitudinal waves are successive and have different effects on the stability of rock slopes with different degrees of weathering. Often only the influence of the horizontal seismic force of the transverse wave is considered [23]. Later, some scholars have shown that the coupling effect of the transverse wave and the longitudinal wave has an amplified effect on the slope stability during the influence of the seismic force on the stability of the rock slope, and the probability of the two acting at the same point at the same time is small, and the actual effect is not large [24], and according to the “Code for Seismic Engineering of Highway Engineering” [25], it is known that the slope in high-intensity areas is significantly affected by the vertical wave vertical seismic action. Therefore, this paper separates the seismic longitudinal and longitudinal waves separately, based on the rock slope of Duwen Highway in Sichuan and based on the calculation results of the slope stability coefficients of different weathering degrees under the seismic longitudinal and longitudinal waves [22].The calculation formulas are shown in formulas (1)∼(8), the calculation parameters are shown in Table 3, and the calculation result is shown in Table 4. The calculation formula of the slope stability coefficient under the shear waves of Models 1∼4 is as follows:(1)Fs1′=mgla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(2)Fs2′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(3)Fs3′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mgla+maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(4)Fs4′=mgcosβ−maxtsinβ−1/18γwe2tanφ+cH−e/sinβmgsinβ+maxtcosβ.Table 3
Calculation parameters.
m (kg)ax (m/s2)ay (m/s2)γw (kN/m3)flk (MPa)fok (MPa)c (MPa)Φ (°)β (°)2.45 × 1059.589.481012.51.120.743890Table 4
Calculation result of the stability factor.
ModelDirectionModel 1Model2Model 3Model 4Shear wave action3.322.151.151.80Longitudinal wave action560.96178.823.600.92The calculation formula of the slope stability coefficient under Model 1∼4 earthquake longitudinal wave is as follows:(5)Fs1″=mg+maytla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(6)Fs2″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(7)Fs3″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mg+maytla+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(8)Fs4″=mgcosβ−maytcosβ−1/18γwe2tanφ+cH−e/sinβmg+maytsinβ,where β(°) is the vertical structural plane dip angle; c is the rock cohesion (MPa); φ is the internal friction angle (°); Gw is the rock’s own weight (kN); Px is the horizontal seismic force (kN); Py is the vertical seismic force (kN); Q is the dangerous rock fissure water pressure (kN); la is the horizontal distance between the center of gravity of the rock block and the weathering point of the cavity; lb is the horizontal distance between the vertical structural surface and the rock cavity surface; e is the fracture penetration length; H is the height of the dangerous rock; h0 is the vertical distance from the center of gravity to the top surface of the base; γw is the bulk density of the fracture water (kN/m3); m is the mass of the rock (kg); ax(t) is the horizontal seismic acceleration of transverse wave (m/s2); flk (MPa) is the standard value of the tensile strength of sandstone, fok (MPa) is the standard value of the tensile strength of mudstone; ay(t) is the vertical seismic acceleration of longitudinal wave (m/s2).It can be seen that the stability coefficient of Model 1 and Model 2 under the action of transverse wave is much smaller than the result under the action of longitudinal wave. Therefore, the seismic shear wave is the main factor of seismic damage of these two slope models. Similarly, the seismic wave with the smaller stability factor is selected as the main factor, that is, the main factor affecting the damage of Model 3. For the shear wave, the main factor affecting the destruction of Model 4 is the longitudinal wave. This result further improves the calculation efficiency.In this paper, the E-W and U-D waves obtained at the Wolongtai station in the 2008 Wenchuan earthquake in Sichuan were used as the transverse and longitudinal waves, and the first 25 s of the seismic time history records was taken as the horizontal and vertical seismic wave input models. The curve is shown below. Model 1∼3 seismic wave action uses shear wave, and Model 4 seismic wave action uses longitudinal wave (see Figure5).Figure 5
Model input seismic wave acceleration time history curve. (a) E-W wave curve (M1, M2, and M3). (b) U-D wave curve (M4) (the dataset is provided by China Earthquake Networks Center, National Earthquake Data Center (http://data.earthquake.cn)).
(a)(b)


### 2.3.3. Boundary Conditions and Monitoring Plan

When analyzing the dynamic response of a slope under earthquake action, you need to enter the seismic wave acceleration time history curve at the bottom of the model. In order to prevent the reflection of seismic waves at the model boundary, a viscous boundary needs to be used. The normal and tangential boundaries of the boundary are provided with dampers. During the loading process, the seismic acceleration is converted into boundary acceleration, which drives the model particles to vibrate.This article is to monitor the law of rock slope movement damage under the action of earthquake. In the PFC data monitoring and analysis, the history command is used to track the data of the object under study. The tracked content includes data such as displacement, coordination number, and energy. Each rock block is grouped, and the movement data of each rock block are monitored in groups. A total of 8 groups are arranged. The boundary conditions and monitoring scheme are shown in Figure6.Figure 6
Boundary conditions of slope dynamic calculation.


## 2.3.1. Geological and Numerical Model

Based on the calculation formula of rock slope stability proposed by Chen et al [21, 22], according to the degree of rock cavity weathering and the degree of crack penetration of the main structural plane, this kind of slope is divided into 4 types. The crack penetration of particle flow slope and rock cavity weathering are shown in Figures 4(a) and 4(b), respectively. Table 2 shows the weathering of each model cavity and fracture.Figure 4
Layout of slope fissure and weathering of rock cavity (unit: m). (a) Fracture (DFN joint) and force chain grid. (b) Weathering layout of particle flow rock cavity.
(a)(b)Table 2
Weathering situation of each model cavity and fissure.
ModelCrack depth (m)Fracture penetration rate (%)Cavity thickness (m)M14.125254.5M26.237.53M38.25501.5M410.3162.50


## 2.3.2. Input of Seismic Waves

In the process of seismic failure of rock slopes, the propagation speeds of shear waves and longitudinal waves are successive and have different effects on the stability of rock slopes with different degrees of weathering. Often only the influence of the horizontal seismic force of the transverse wave is considered [23]. Later, some scholars have shown that the coupling effect of the transverse wave and the longitudinal wave has an amplified effect on the slope stability during the influence of the seismic force on the stability of the rock slope, and the probability of the two acting at the same point at the same time is small, and the actual effect is not large [24], and according to the “Code for Seismic Engineering of Highway Engineering” [25], it is known that the slope in high-intensity areas is significantly affected by the vertical wave vertical seismic action. Therefore, this paper separates the seismic longitudinal and longitudinal waves separately, based on the rock slope of Duwen Highway in Sichuan and based on the calculation results of the slope stability coefficients of different weathering degrees under the seismic longitudinal and longitudinal waves [22].The calculation formulas are shown in formulas (1)∼(8), the calculation parameters are shown in Table 3, and the calculation result is shown in Table 4. The calculation formula of the slope stability coefficient under the shear waves of Models 1∼4 is as follows:(1)Fs1′=mgla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(2)Fs2′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb2/2maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(3)Fs3′=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mgla+maxth0+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(4)Fs4′=mgcosβ−maxtsinβ−1/18γwe2tanφ+cH−e/sinβmgsinβ+maxtcosβ.Table 3
Calculation parameters.
m (kg)ax (m/s2)ay (m/s2)γw (kN/m3)flk (MPa)fok (MPa)c (MPa)Φ (°)β (°)2.45 × 1059.589.481012.51.120.743890Table 4
Calculation result of the stability factor.
ModelDirectionModel 1Model2Model 3Model 4Shear wave action3.322.151.151.80Longitudinal wave action560.96178.823.600.92The calculation formula of the slope stability coefficient under Model 1∼4 earthquake longitudinal wave is as follows:(5)Fs1″=mg+maytla+flkH−e2/2sin2β+H−elb/tanβ+fok⋅lb2/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(6)Fs2″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/21/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(7)Fs3″=flkH−e2/2sin2β+H−elb/tanβ+lbfok⋅lb/2mg+maytla+1/18γwe2e/9sinβ+H−e/sinβ+lbcosβ,(8)Fs4″=mgcosβ−maytcosβ−1/18γwe2tanφ+cH−e/sinβmg+maytsinβ,where β(°) is the vertical structural plane dip angle; c is the rock cohesion (MPa); φ is the internal friction angle (°); Gw is the rock’s own weight (kN); Px is the horizontal seismic force (kN); Py is the vertical seismic force (kN); Q is the dangerous rock fissure water pressure (kN); la is the horizontal distance between the center of gravity of the rock block and the weathering point of the cavity; lb is the horizontal distance between the vertical structural surface and the rock cavity surface; e is the fracture penetration length; H is the height of the dangerous rock; h0 is the vertical distance from the center of gravity to the top surface of the base; γw is the bulk density of the fracture water (kN/m3); m is the mass of the rock (kg); ax(t) is the horizontal seismic acceleration of transverse wave (m/s2); flk (MPa) is the standard value of the tensile strength of sandstone, fok (MPa) is the standard value of the tensile strength of mudstone; ay(t) is the vertical seismic acceleration of longitudinal wave (m/s2).It can be seen that the stability coefficient of Model 1 and Model 2 under the action of transverse wave is much smaller than the result under the action of longitudinal wave. Therefore, the seismic shear wave is the main factor of seismic damage of these two slope models. Similarly, the seismic wave with the smaller stability factor is selected as the main factor, that is, the main factor affecting the damage of Model 3. For the shear wave, the main factor affecting the destruction of Model 4 is the longitudinal wave. This result further improves the calculation efficiency.In this paper, the E-W and U-D waves obtained at the Wolongtai station in the 2008 Wenchuan earthquake in Sichuan were used as the transverse and longitudinal waves, and the first 25 s of the seismic time history records was taken as the horizontal and vertical seismic wave input models. The curve is shown below. Model 1∼3 seismic wave action uses shear wave, and Model 4 seismic wave action uses longitudinal wave (see Figure5).Figure 5
Model input seismic wave acceleration time history curve. (a) E-W wave curve (M1, M2, and M3). (b) U-D wave curve (M4) (the dataset is provided by China Earthquake Networks Center, National Earthquake Data Center (http://data.earthquake.cn)).
(a)(b)


## 2.3.3. Boundary Conditions and Monitoring Plan

When analyzing the dynamic response of a slope under earthquake action, you need to enter the seismic wave acceleration time history curve at the bottom of the model. In order to prevent the reflection of seismic waves at the model boundary, a viscous boundary needs to be used. The normal and tangential boundaries of the boundary are provided with dampers. During the loading process, the seismic acceleration is converted into boundary acceleration, which drives the model particles to vibrate.This article is to monitor the law of rock slope movement damage under the action of earthquake. In the PFC data monitoring and analysis, the history command is used to track the data of the object under study. The tracked content includes data such as displacement, coordination number, and energy. Each rock block is grouped, and the movement data of each rock block are monitored in groups. A total of 8 groups are arranged. The boundary conditions and monitoring scheme are shown in Figure6.Figure 6
Boundary conditions of slope dynamic calculation.


## 3. Results and Discussion

### 3.1. Rock Triaxial Test Analysis

#### 3.1.1. Comparative Analysis of Rock Stress-Strain

The stress-strain curves under different confining pressures obtained from the rock chamber and the numerical triaxial test are shown in Figures7 and 8. Compared with the two stress-strain curves, the trend is basically similar. It can be seen that the peak bias stress error of sandstone is kept within 0∼1 MPa and the peak strain error is controlled within 0.1%, so the parameters are reasonable.Figure 7
Conventional triaxial stress-strain curves of two rocks. (a) Mudstone. (b) Sandstone.
(a)(b)Figure 8
Stress-strain curve of rock numerical triaxial test. (a) Mudstone. (b) Sandstone.
(a)(b)


#### 3.1.2. Comparative Analysis of Rock Strength Attenuation

The rock compressive strength attenuation value and strength attenuation coefficient reflect the brittleness degree of the rock. In this paper, we cite the mechanical indexes proposed by Peng et al. [15] to characterize the post-peak strength decay behavior of rock samples—rock strength attenuation coefficient:(9)Ds=σ1P−σrσ1P=Δσσ1P,(10)Δσ=σ1P−σr,where Ds is the intensity attenuation coefficient, value range [0, 1]; Δσ is the intensity attenuation value; σ1P is the peak intensity; σr is the residual intensity.During the trial adjustment of mesoparameters, the comparison curves of the relationship between the compressive strength attenuation value and attenuation coefficient of sandstone and mudstone calculated according to equations (9) and (10) are shown in Figures 9(a) and 9(b). It can be found that with the increase of confining pressure, the strength attenuation value of sandstone shows a gradual decreasing trend, which is consistent with the conclusion obtained in [26], and verifies the rationality of this conventional triaxial test. The sandstone attenuation coefficient curve decreased sharply after 2 MPa confining pressure, from 0.67 to 0.37, and then decreased slowly, indicating that the brittleness of sandstone gradually weakened with the increase of confining pressure, and the brittleness weakening ability increased after 2 MPa. In Figures 9(a), the curve of mudstone strength attenuation value and attenuation coefficient changes in a V-shape. When the curve drops to 4 MPa confining pressure, the attenuation value and attenuation coefficient increase simultaneously, rising 7.05 and 0.15, respectively, and the rising curve slows down after 6 MPa. It shows that at a low confining pressure of 2 MPa, the brittleness of the mudstone first weakened due to the large initial porosity, and the ductility increased. At 4 MPa, the brittleness increased and the ductility decreased. When the confining pressure is small, the two rock samples have a large attenuation coefficient and high brittleness. They are basically one-time shear failure, and there are few local cracks, and some tensile cracks will appear; when the confining pressure is large, the attenuation coefficient of the rock samples is relatively reduced, brittleness is weakened, and there are many local cracks. The literature [27] also confirmed this law.Figure 9
Comparison of rock attenuation value and attenuation coefficient. (a) Mudstone. (b) Sandstone.
(a)(b)Comparing the results of the laboratory test and numerical test, it is found that the attenuation value and attenuation coefficient of compressive strength of mudstone are relatively low in the numerical test and those of sandstone are relatively high. The parameters mainly depend on the parallel bonding modulus and the critical damping ratio, while the mudstone modulus is relatively low and the critical damping is relatively low, resulting in low brittleness in the numerical experiment; the sandstone modulus is high and the critical damping is relatively high, resulting in high brittleness.


#### 3.1.3. Comparative Analysis of Rock Failure Process

During the numerical triaxial test, PFC3D was used to track the broken contact force chain in the parallel bonding contact model through the DFN joint system, and the broken contact link was expressed in the form of DFN_disk (disk joint), thereby forming a damage crack. In DFN, the internal cracks of rock particles are divided into tensile cracks and shear cracks. During the loading of the model, the type of microcracks can be identified and the number of them can be counted by programming the FISH function. The mesoscopic tensile cracks and shear cracks do not directly reflect the tensile and shear cracks of the macroscopic failure of the rock sample, but in the three-dimensional particle flow model, different mesoscopic crack spatial trajectories can observe the development of the macroscopic damage cracks of the rock. The disc joint reflects the damage characteristics of the mesoscopic particles of the rock. This method is especially important in the analysis of the failure characteristics of the rock sample and the precise trial adjustment of the mesoscopic parameters.The sandstone and mud numerical triaxial test failure shear cracks under different confining pressure conditions and loading methods are shown in Figures10 and 11. In this paper, only the failure test graph of the specimen under the minimum confining pressure of 2 MPa and the maximum confining pressure of 8 MPa is compared with the PFC numerical simulation failure map. From Figure 7, it can be seen that the mudstones mostly show shear cracks on the microjoints, and the fractures between the particles are mostly shear damage. At a confining pressure of 2 MPa, the whole mudstone exhibits a single shear failure, with fewer shear failure cracks, and brittleness characteristics are more obvious, accompanied by a small amount of tensile cracks from the bottom up; as the confining pressure increases, the mudstone exhibits mixed shear failure, multiple shear cracks appear in the sample, the brittleness is weakened, the ductility is enhanced, and the damage cracks increase. When the confining pressure continues to increase to 8 MPa, the microscopic cracks at the bottom and top gradually increase and develop toward the middle, and the shear failure cracks increase and become more obvious. The cracks increase and become more obvious. The development trend of mudstone cracks is similar to the indoor test, and the development characteristics are consistent with the results of the indoor test. It can be seen from Figure 8 that the sandstone shows more tensile cracks on the microjoints, and the tensile fracture failure is mainly among the particles. When the confining pressure is 2 MPa, the mesocracks are generated from the bottom of the sandstone and gradually develop upwards. Due to the low confining pressure, the macroshear crack did not penetrate to the top, and the top microcrack developed slowly; when the confining pressure increased to 8 MPa, the bottom crack gradually expanded to the top, forming a penetrating macroscopic shear crack. As the confining pressure increases, the number gradually increases, forming multiple shear cracks that alternately develop failure characteristics. The development trend of sandstone cracks is also similar to laboratory tests. Through the three-dimensional crack tracing observation of the numerical triaxial test of the particle flow simulation above, it can be seen that the macrofailure characteristics of the two rocks under the PFC numerical test are similar to the indoor test, which further verifies the reliability of the mesoparameter adjustment in this paper.Figure 10
Numerical triaxial failure joints (cracks) of mudstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)Figure 11
Triaxial failure joints (cracks) of sandstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)In summary, the macroscopic mechanical properties of the mesoparameters adjusted by the three-axis numerical test have certain reliability, and the materials are truly usable. The mesoparameters of the obtained sandstone and mudstone materials and the mesoparameters of structural plane are shown in Tables5 and 6, respectively.Table 5
Mesoscopic parameters of sandstone and mudstone materials.
Lithologyd (mm)Ρ (kg·m−3)Ec(GPa)kc∗μλk∗Eb(GPa)σc (MPa)c (MPa)βnβsMudstone0.24310021.00.81.001.02018280.20.2Sandstone0.26300031.00.61.011.03020350.30.6Table 6
Mesoscopic parameters of structural plane.
Structural plane typesj_kn (N·m−1)sj_ks (N·m−1)sj_μsj_c (MPa)sj_φ (°)Horizontal structural plane1010.70.310Vertical structural plane1010.308


### 3.2. Slope Particle Flow Simulation and Failure Mode

#### 3.2.1. Slope Failure Process

After the establishment of four slope models with different cavity depths and fissure depths, the earthquake action is applied to them. The failure process of each model after 25 s of earthquake duration is shown in Figure11, and the displacement of the rock block in the vertical z-direction and horizontal y-direction is shown in Figure 12. It can be seen that after 25 s of seismic waves, all the slope models have significant structural plane fracture expansion and penetration phenomena, of which Model 2 and Model 3 both have obvious rock mass collapse and overall failure phenomena, and the main control plane surface failure. The degree gradually increases with the deepening of the rock cavity and the initial fissures; the onset time of failure under the earthquake action of each model is different. 2∼2.5 s began to fail, indicating that as the weathering degree of the rock cavity deepened, the initial crack penetration degree increased, and the faster the rock slope collapse rate under earthquake.Figure 12
The destruction process diagram of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The failure process of Model 1 seismic wave within 25 seconds is shown in Figure11(a). It can be seen that there is no sign of rock collapse during the earthquake period of this model. The displacement values are relatively small. The negative vertical displacements are distributed between 0 and 0.4 m, and the positive horizontal displacements are distributed between 0 and 1.2 m. This indicates that the damage characteristics of Model 1 in the horizontal direction are relatively large. When the earthquake action time is between 0 and 8.5 seconds, there is no vertical displacement change, and each rock block has a horizontal y-direction displacement fluctuation with the shear wave action, and the amplitude of the fluctuation increases with the extension of the earthquake action time. At the same time, there is no obvious deformation and destruction of the two rock mass macrochains, and there is no crack development in the mudstone. At this time, the slope is relatively stable. When the earthquake acts for 9 seconds, the seismic acceleration curve gradually reaches the peak, and the displacement of each rock block in Model 1 has a small frequency fluctuation in the vertical direction. At this time, the contact model of the main control structure near the 1-1 rock block and the 2-2 rock block is locally fractured. Local fracture occurred in the contact model of the main control structure. The first and second layers of rock began to move, and the lateral force of the smooth joint model between the layers increased and gradually destroyed. After 10 seconds, the structure of the no. 1-1 rock block penetrated upward. Therefore, the structural surface of the no. 2-1 rock block fractured, and the crack propagation in the rock cavity intensified. The mudstone particles disintegrated severely. After 11.5 seconds, the vertical displacement curve of Model 1 fell slowly downward. The horizontal displacement to the right gradually increases, and the rock strata intensify. Rock no. 3-2 and no. 1-2 have large vertical displacement due to local fragmentation on the top, which can reach a maximum of 1.35 m. Due to the interaction of adjacent rock blocks at the bottom under earthquake, the structural surface of the no. 4-1 rock block at the top layer was partially broken. After 25 seconds, the upper rock mass was destroyed except for the no. 4-2 rock block, the slope as a whole had a tendency to collapse, and the maximum horizontal displacement of the no. 4-1 rock block was 1.2 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 25% and a mudstone cavity depth of 1.5 m will undergo slight instability and failure under the action of a shear wave of 25 seconds. The destruction mode is fracture-horizontal sliding failure.The failure process of Model 2 seismic wave action within 25 s is shown in Figure12(b). The displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(b) and 13(f). At 8 seconds after the earthquake, due to the deepening of the weathering of the rock cavity, there is no support at the bottom of the no. 1 macrochain. The local fracture of the no. 1-1 rock block preferentially occurs, and the horizontal and vertical displacements increase sharply. After 9 seconds, the internal chain no. 1-2 rock, no. 2-2 rock, and no. 3-2 rock structural surface broke. A partial penetration surface is formed, and then the top of the rock block is shifted outward. After 11 seconds, the no. 1-1 rock collided with the slope for the first time, and the displacement curve suddenly changed. Rock 3-1, rock 3-2, rock 4-1, and rock 4-2 were dumped-rotated and their horizontal and vertical curve slopes increased simultaneously. The maximum horizontal displacement can reach 7.8 m, and the maximum vertical displacement can reach 6.3 m. Rock 2-1 began to collapse due to the earthquake and the squeezing effect of the adjacent rock mass. At the same time, the lower mudstones have cracked and damaged. When the earthquake took place for 19 seconds, the no. 2-1 rock block collided with the no. 1-1 rock block for the first time. Under the action of the earthquake, the top four rock blocks were rotated and destroyed around the rock block no. 2-2. The earthquake lasted from 22 seconds to 25 seconds, and the rock mass collapsed and intensified. The top rock mass separated between layers. Rocks 3-1 and 3-2 simultaneously fell and collided with rock 2-1. Rocks 4-1 and 4-2 rolled outward. Until the end of the calculation, the maximum vertical displacement of rock no. 1-1 is 4.5 m and the maximum horizontal displacement is 2.5 m. The maximum vertical displacement of no. 2-1 rock block is 4 m, and the maximum horizontal displacement is 0.8 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 37.5% and a mudstone cavity depth of 3 m will undergo cracking-caving-dumping failure under the action of a 25-second shear wave.Figure 13
Displacement time history curve of each model. (a) Model 1,y-direction. (b) Model 2, y-direction. (c) Model 3, y-direction. (d) Model 4, y-direction. (e) Model 1, z-direction. (f) Model 2, z-direction. (g) Model 3, z-direction. (h) Model 4, z-direction.
(a)(b)(c)(d)(e)(f)(g)(h)The damage process of Model 3 seismic wave action within 25 s is shown in Figure12(c), and the horizontal y-direction and vertical z-direction displacements are shown in Figures 13(c) and 13(g). After 5 seconds, the crack of no. 1-1 rock’s main control surface quickly initiated and propagated, and local fracture and shedding occurred. After 6 seconds, rock 1-1 collapsed, followed by rocks 4-2 and 3-1. The horizontal displacement of the rock block under the action of earthquake has changed obviously. After 7 seconds, the crack at the no. 2 macrochain at the top of the slope developed downward. As a result, no. 3-2 and no. 2-2 rocks fractured and separated from the parent rock. At this time, the destruction of the no. 3-1 rock block at the top prompted the rapid destruction of the no. 2-1 rock structural plane. The two rocks collapsed at the same time and collided with the no. 1-1 rock. After 10 seconds, as the seismic force continued to increase, the bottom no. 1-2 rock collapsed and broke away from the parent rock. At this time , the internal main control structure has been completely broken and penetrated, and the rock cavity and mudstone particles were destroyed. After 15 seconds, the no. 4-1 rock dumped and destroyed along with the four internal rock blocks and collided with the no. 3-1 rock. Due to the excessive torque of the top rock block, at 20 seconds, the rock joint model between no. 3-2 rock and no. 4-2 rock was destroyed. And no. 4-1 and no. 4-2 rocks rotate and destroy around no. 3-2. At the end of 25 seconds, the top two rock masses turned and destroyed and collided with no. 3-1 rock mass, which has not yet completely collapsed. Throughout the process, the maximum vertical displacement of no. 1–1 rock collapse can reach 4 m, and the horizontal roll-off displacement can reach 1.7 m. The top no. 4-1 and no. 4-2 rocks have a large roll-off tendency. Therefore, the failure mode of a complex horizontal layered rock slope with a fracture penetration rate of 50% and a mudstone cavity depth of 4.5 m under the action of a shear wave of 25 s is caving-dumping-rotation failure. The development of fissures and the weathering of rock cavities are more serious.The weathering degree of Model 4 is the most serious. The destruction process of seismic wave action within 25 s is shown in Figure12(d), and the displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(d) and 13(h). Within 0∼2 seconds, the slope was in a stable state. After 2.5 seconds, the no. 1-1 rock began to fall, and the vertical displacement curve increased significantly, accompanied by horizontal displacement. After 3 seconds, the no. 1-2 rock collapsed and was squeezed by the upper rock mass. After 5 seconds, the no. 1-1 rock rolled down along the surface of the mudstone slope, with a vertical fall displacement of 3.8 m and a horizontal roll-off displacement of 1.45 m. Due to the continued action of earthquakes and interlayer loads, the no. 2-1 rock broke away from the parent rock and began to fall. The earthquake action gradually increased, and the squeezing effect of no. 2-2 rock and no. 1-2 rock on the upper rock mass also continued to increase. After 8 seconds, the slope began to slip overall. After 11 seconds, the no. 3-1 rock began to fall under the action of the earthquake. Within 15∼25 seconds, the joint model of the top rock layer began to break, and the no. 4-1 and no. 4-2 rocks were rotationally damaged. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. At last, Model 4 suffered the most severe damage and the fastest damage speed after 25 seconds of earthquake action. And there is no rolling movement of single rock. The overall trend of slope collapse is more significant, and the failure mode is caving-slip-rotation failure.In summary, under the action of seismic wave in one direction, the damage characteristics of complex horizontal layered rock slopes are affected by the depth of the cavity and the crack penetration rate of the main control surface. The basic failure mode (see Figure12) can be divided into fracture-horizontal slid failure (Model 1), tensile cracking-caving-dumping failure (Model 2), caving-dumping-rotation failure (Model 3), and caving-slip-rotation failure (Model 4). The greater the weathering degree of the rock cavity and fissures, the more obvious the tendency of rock avalanches and the faster the avalanche time. When the rock blocks are collapsing, the horizontal and interlayer adjacent rock blocks will squeeze and collide with each other, which will intensify the fracture of the smooth joint model. The slope has developed from the single rock block’s pull failure to the overall dumping (slippage) rotation failure. In the work of geological disasters, according to the failure modes and displacement characteristics of the slopes obtained in this paper, the reinforcement of Model 1 type slopes can be strengthened, and the rock masses of Model 2∼4 type slopes can be prevented.


#### 3.2.2. Coordination Number Analysis

In order to further reveal the seismic failure process and displacement characteristics of the slope, the specific strength attenuation of the contact model during the fracture process and the specific rock collapse time need to be quantitatively described. Therefore, the concept of coordination number is introduced in this paper. Coordination number, that is, the average number of particles in contact is equal to twice the number of particles in contact/number of particles. This parameter comprehensively reflects the physical properties of the particle contact model (friction coefficient, elastic modulus, Poisson’s ratio, etc.). The strength is closely related [28]. In the simulation process of PFC rock fracture failure, it can reflect the fracture time and strength attenuation of the contact force chain of the specific structural surface. In this paper, the particle coordination values in units of rocks in each model are monitored, and the time history curve is plotted (see Figure 14). It can be seen from the figure that the coordination values of Model 1 and Model 4 are distributed between 1–11 and 6–11 respectively, and the coordination values of Model 2 and Model 3 are distributed between 5 and 11, which can be seen as the depth of the cavity and the cracks. As it deepens, the change range of the coordination value gradually decreases. The reason is that the initial stability of the slope gradually decreases. The coordination number attenuation of each model slope in the period of anticavitation collapse is about 2∼15 s, and as the weathering deepens, the coordination number attenuation effect shows a certain advance effect, that is, Models 1∼4 correspond to 8.7∼11.4 s, 7.1∼10.7 s, 5.4∼8.3 s, and 2.3∼11.3 s respectively. In Model 1, the coordination number of rock block 1-2 is the most attenuated by earthquake, from 9.2 s from 10.5 to 1.8. The reason is that during this period, the rock has undergone large local fragmentation and destruction, and some particles break away. For the volume, the attenuation value includes not only the fracture of the structural plane, but also the reduction of the coordination number caused by the fracture of the particles inside the rock. The attenuation values of coordination number of other rock blocks are maintained between 2 and 4. As the No. 4-2 rock mass has not collapsed 25 s after the earthquake, the coordination number attenuation value is at least 2. In each model, as the elevation increases, the coordination number attenuation value decreases, of which the no. 4-1 and no. 4-2 rock bodies placed on the top of Model 4 have the lowest coordination number attenuation value of 0.6 and 1.2. It can be seen that the coordination number has an elevation attenuation effect in dangerous rock, and with the development of the cavity and the depth of the fissure, the elevation attenuation effect gradually weakens. Comparing the coordination number curves of each model, it can be seen that with the deepening of the cavity and cracks, the stability of the main structural surface of the rock slope decreases, the start time of the fluctuation of the coordination number gradually decreases, and the fluctuation range of the coordination value gradually increases. After reaching Model 5, due to the intensified collision motion of rocks, the particle contact behavior between different rocks is more significant, resulting in an increasing trend of the coordination number curve and further increase in fluctuations. It can be seen that the attenuation speed and attenuation degree of the structural plane under the action of an earthquake are determined by the initial stability of the slope (i.e., the depth of the rock cavity and the depth of the crack). The development trend of the coordination number of each model is consistent with the failure process diagram and displacement diagram, which further confirms the specific failure time of each rock block of each model, as shown in Table 7. According to this, it can provide a strong reference for the treatment of horizontal complex layered rock slopes under different working conditions (Models 1 to 4).Figure 14
Time history curve of coordination number of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)Table 7
Fracture failure timetable of each model rock block.
ModelRock block1-1 (s)1-2 (s)2-1 (s)2-2 (s)3-1 (s)3-2 (s)4-1 (s)4-2 (s)Model 18.79.29.68.79.38.911.4>25Model 27.18.410.78.99.18.89.89.7Model 35.48.36.26.35.96.15.85.8Model 42.33.04.64.911.38.09.06.1


#### 3.2.3. Energy Dissipation Analysis

The laws of thermodynamics show that energy conversion is an essential feature of the physical change process of matter, and material destruction is a state instability phenomenon driven by energy. The energy dissipation during the failure of the slope is mainly reflected in the failure of the structural plane and the damage and rupture of the intact rock block, reflecting the internal mechanism of the failure evolution. This article assumes that the rock slope system is a closed system that has no heat exchange with the outside world, and it can be obtained according to the first law of thermodynamics [29]:(11)U=Ud+Ue,where U is the total input energy generated by the work of the slope’s body force during the earthquake; Ud is the dissipated energy during the dynamic response of the slope; Ue is the elastic strain energy that can be released inside the slope particles. The total input energy U in the model of this paper is obtained by the following formula:(12)U=Upre+∑NpγpVpΔu,where Upre is the total input energy accumulated in the last calculation time step; γp and Δu are the particle weight and displacement in the current time step, respectively; Np is the number of particles. In the model, the sandstone and mudstone particle contact model is a parallel bonding model. The elastic strain energy Ue consists of the particle strain energy Uc and the parallel bond strain energy Upb, namely,(13)Ue=Uc+Upb,Uc=12∑NcFin2/kn+Fis2/ks,Upb=12∑NpbF¯in2/Ak¯n+F¯is2/Ak¯s+M¯2/Ik¯n,where Nc and Npb are the number of contact force chains and the number of parallel bonding; i represents the i-th contact chain; Fin and Fis are the normal contact force and tangential contact force; F¯in, F¯is, and M¯ are the normal parallel bonding force, tangential parallel bonding force, and parallel bonding moment; kn and ks are the contact normal stiffness and contact tangential stiffness; k¯n, k¯s, and k¯n are, respectively, the stiffness corresponding to each parallel bonding force; A and I are the cross-sectional area and moment of inertia of the parallel bonding.The energy dissipationUd of the final model can be obtained by substituting various equations into the solution. According to the first law of thermodynamics and the particle flow comes with history energy detection system, the instantaneous kinetic energy of the slope can be obtained.The kinetic energy response of each model under the action of earthquake is shown in Figure15. It can be seen that from Model 1 to Model 4, the maximum kinetic energy decreased from 155 MJ to 72 MJ, indicating that as the development of the cavity and cracks deepens, the energy generated by the failure of the slope gradually decreases. This is due to the deepening of the cavity and the upper and lower sandstones. The friction with the mudstone particles is reduced, and the degree of mudstone particle disintegration is reduced; the cracks are deepened, and the flat bond chain of the structural surface is reduced, that is, the energy reflection is reduced. In Models 1 to 3, the kinetic energy fluctuations gradually increase with the duration of the earthquake, and the energy gradually increases with the severity of the slope damage. Among them, Model 1 reached the first peak of 46 MJ at 6 s, indicating that the contact force chain began to break, and then the peak fluctuated between 75 and 90 MJ during the caving process, which was caused by the collision of the rocks. As a result, the maximum peak reached 155 MJ at 13 s; Model 2 produced a peak of 45 MJ kinetic energy at 2 s due to local particle collisions and subsequently produced peaks of 85 MJ and 89 MJ at 13 s and 15 s; similarly, in Model 3 at 15 s, the instantaneous energy of 104 MJ was generated due to collision of rock particles. In Model 4, the kinetic energy response increases rapidly to a maximum value of 73 MJ with the duration of the earthquake and then decreases to fluctuate between 0 and 43 MJ. It can be seen that due to the strong stability of Models 1 to 3, the dangerous rock gradually collapses and breaks down over time, the slope shows progressive failure, and the kinetic energy response gradually increases. Because of the weak stability of Model 4, the whole slope is in the early stage When the destruction is completed in a short time, the kinetic energy response reaches the maximum value at the initial stage, and then it collapses over time, which is in good synchronization with the failure mode of the slope.Figure 15
Variation of kinetic energy response of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The energy dissipation evolution process of each model is shown in Figure16. It can be seen that Model 1 and Model 4 produce the most total energy. The maximum energy dissipation of Models 1 to 4 is 105 MJ, 76 MJ, 63 MJ, and 82 MJ, respectively. Model 1 has the largest total energy dissipated due to more structural surface contact, which means the dissipated energy is the largest. The interaction between the blocks is large, so the dissipation energy is higher than that of Model 2 and Model 3. The energy dissipation of each slope model along with the development of earthquake damage and destruction can be divided into two stages, and the whole has the characteristics of first followed by slow. The first stage is the rock block fracture and collision stage, which corresponds to a period of 2.5∼15 s, the period when the dangerous rock blocks break away from the parent rock and collides with each other. At this stage, the change of dissipated energy of Model 1 is 0∼100 MJ, accounting for 95% of the total energy; the change of dissipated energy of Model 2 is 0∼68 MJ, accounting for 89% of the total energy; that of Model 3 is 0∼56 MJ, accounting for 89% of the total energy; that of Model 4 is 74 MJ, accounting for 90% of the total energy. At this stage, the dissipated energy curve develops faster and faster with the duration of the earthquake. Corresponding to the process of “cracking-caving” and “dumping-collision” in the failure mode, it is first “cracking-caving” and then “dump-collision.” And the energy dissipation caused by the dumping collision is greater than the energy generated by the fracture of the structural plane, and it also reflects the progressive cumulative energy dissipation mechanism of the development of dangerous rock failure. The second stage is the stage of the overall collapse of the slope. In this stage, the rock mass has completed local fracture and gradually fell to the foot of the slope. At this stage, the energy dissipation change of Model 1 is 100–105 MJ, accounting for 5% of the total energy. That of Model 2 is 68–76 MJ, of Model 3 is 56–63 MJ, and of Model 4 is 74–82 MJ. The change with time is slower and the change value is smaller, corresponding to the process of “overall dumping” and “overall seat sliding” in the failure mode. Comprehensive analysis shows that the trend of energy dissipation is closely related to the failure mode. The degree of overall energy dissipation affects the cavity and the depth of the fissures. There is a significant difference between the energy consumption of the first stage and the second stage. The energy consumption of horizontal complex layered rock slopes under earthquake action is more concentrated in the first stage of rock fracture collision.Figure 16
Comparison of energy dissipation of different models.


## 3.1. Rock Triaxial Test Analysis

### 3.1.1. Comparative Analysis of Rock Stress-Strain

The stress-strain curves under different confining pressures obtained from the rock chamber and the numerical triaxial test are shown in Figures7 and 8. Compared with the two stress-strain curves, the trend is basically similar. It can be seen that the peak bias stress error of sandstone is kept within 0∼1 MPa and the peak strain error is controlled within 0.1%, so the parameters are reasonable.Figure 7
Conventional triaxial stress-strain curves of two rocks. (a) Mudstone. (b) Sandstone.
(a)(b)Figure 8
Stress-strain curve of rock numerical triaxial test. (a) Mudstone. (b) Sandstone.
(a)(b)


### 3.1.2. Comparative Analysis of Rock Strength Attenuation

The rock compressive strength attenuation value and strength attenuation coefficient reflect the brittleness degree of the rock. In this paper, we cite the mechanical indexes proposed by Peng et al. [15] to characterize the post-peak strength decay behavior of rock samples—rock strength attenuation coefficient:(9)Ds=σ1P−σrσ1P=Δσσ1P,(10)Δσ=σ1P−σr,where Ds is the intensity attenuation coefficient, value range [0, 1]; Δσ is the intensity attenuation value; σ1P is the peak intensity; σr is the residual intensity.During the trial adjustment of mesoparameters, the comparison curves of the relationship between the compressive strength attenuation value and attenuation coefficient of sandstone and mudstone calculated according to equations (9) and (10) are shown in Figures 9(a) and 9(b). It can be found that with the increase of confining pressure, the strength attenuation value of sandstone shows a gradual decreasing trend, which is consistent with the conclusion obtained in [26], and verifies the rationality of this conventional triaxial test. The sandstone attenuation coefficient curve decreased sharply after 2 MPa confining pressure, from 0.67 to 0.37, and then decreased slowly, indicating that the brittleness of sandstone gradually weakened with the increase of confining pressure, and the brittleness weakening ability increased after 2 MPa. In Figures 9(a), the curve of mudstone strength attenuation value and attenuation coefficient changes in a V-shape. When the curve drops to 4 MPa confining pressure, the attenuation value and attenuation coefficient increase simultaneously, rising 7.05 and 0.15, respectively, and the rising curve slows down after 6 MPa. It shows that at a low confining pressure of 2 MPa, the brittleness of the mudstone first weakened due to the large initial porosity, and the ductility increased. At 4 MPa, the brittleness increased and the ductility decreased. When the confining pressure is small, the two rock samples have a large attenuation coefficient and high brittleness. They are basically one-time shear failure, and there are few local cracks, and some tensile cracks will appear; when the confining pressure is large, the attenuation coefficient of the rock samples is relatively reduced, brittleness is weakened, and there are many local cracks. The literature [27] also confirmed this law.Figure 9
Comparison of rock attenuation value and attenuation coefficient. (a) Mudstone. (b) Sandstone.
(a)(b)Comparing the results of the laboratory test and numerical test, it is found that the attenuation value and attenuation coefficient of compressive strength of mudstone are relatively low in the numerical test and those of sandstone are relatively high. The parameters mainly depend on the parallel bonding modulus and the critical damping ratio, while the mudstone modulus is relatively low and the critical damping is relatively low, resulting in low brittleness in the numerical experiment; the sandstone modulus is high and the critical damping is relatively high, resulting in high brittleness.


### 3.1.3. Comparative Analysis of Rock Failure Process

During the numerical triaxial test, PFC3D was used to track the broken contact force chain in the parallel bonding contact model through the DFN joint system, and the broken contact link was expressed in the form of DFN_disk (disk joint), thereby forming a damage crack. In DFN, the internal cracks of rock particles are divided into tensile cracks and shear cracks. During the loading of the model, the type of microcracks can be identified and the number of them can be counted by programming the FISH function. The mesoscopic tensile cracks and shear cracks do not directly reflect the tensile and shear cracks of the macroscopic failure of the rock sample, but in the three-dimensional particle flow model, different mesoscopic crack spatial trajectories can observe the development of the macroscopic damage cracks of the rock. The disc joint reflects the damage characteristics of the mesoscopic particles of the rock. This method is especially important in the analysis of the failure characteristics of the rock sample and the precise trial adjustment of the mesoscopic parameters.The sandstone and mud numerical triaxial test failure shear cracks under different confining pressure conditions and loading methods are shown in Figures10 and 11. In this paper, only the failure test graph of the specimen under the minimum confining pressure of 2 MPa and the maximum confining pressure of 8 MPa is compared with the PFC numerical simulation failure map. From Figure 7, it can be seen that the mudstones mostly show shear cracks on the microjoints, and the fractures between the particles are mostly shear damage. At a confining pressure of 2 MPa, the whole mudstone exhibits a single shear failure, with fewer shear failure cracks, and brittleness characteristics are more obvious, accompanied by a small amount of tensile cracks from the bottom up; as the confining pressure increases, the mudstone exhibits mixed shear failure, multiple shear cracks appear in the sample, the brittleness is weakened, the ductility is enhanced, and the damage cracks increase. When the confining pressure continues to increase to 8 MPa, the microscopic cracks at the bottom and top gradually increase and develop toward the middle, and the shear failure cracks increase and become more obvious. The cracks increase and become more obvious. The development trend of mudstone cracks is similar to the indoor test, and the development characteristics are consistent with the results of the indoor test. It can be seen from Figure 8 that the sandstone shows more tensile cracks on the microjoints, and the tensile fracture failure is mainly among the particles. When the confining pressure is 2 MPa, the mesocracks are generated from the bottom of the sandstone and gradually develop upwards. Due to the low confining pressure, the macroshear crack did not penetrate to the top, and the top microcrack developed slowly; when the confining pressure increased to 8 MPa, the bottom crack gradually expanded to the top, forming a penetrating macroscopic shear crack. As the confining pressure increases, the number gradually increases, forming multiple shear cracks that alternately develop failure characteristics. The development trend of sandstone cracks is also similar to laboratory tests. Through the three-dimensional crack tracing observation of the numerical triaxial test of the particle flow simulation above, it can be seen that the macrofailure characteristics of the two rocks under the PFC numerical test are similar to the indoor test, which further verifies the reliability of the mesoparameter adjustment in this paper.Figure 10
Numerical triaxial failure joints (cracks) of mudstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)Figure 11
Triaxial failure joints (cracks) of sandstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)In summary, the macroscopic mechanical properties of the mesoparameters adjusted by the three-axis numerical test have certain reliability, and the materials are truly usable. The mesoparameters of the obtained sandstone and mudstone materials and the mesoparameters of structural plane are shown in Tables5 and 6, respectively.Table 5
Mesoscopic parameters of sandstone and mudstone materials.
Lithologyd (mm)Ρ (kg·m−3)Ec(GPa)kc∗μλk∗Eb(GPa)σc (MPa)c (MPa)βnβsMudstone0.24310021.00.81.001.02018280.20.2Sandstone0.26300031.00.61.011.03020350.30.6Table 6
Mesoscopic parameters of structural plane.
Structural plane typesj_kn (N·m−1)sj_ks (N·m−1)sj_μsj_c (MPa)sj_φ (°)Horizontal structural plane1010.70.310Vertical structural plane1010.308


## 3.1.1. Comparative Analysis of Rock Stress-Strain

The stress-strain curves under different confining pressures obtained from the rock chamber and the numerical triaxial test are shown in Figures7 and 8. Compared with the two stress-strain curves, the trend is basically similar. It can be seen that the peak bias stress error of sandstone is kept within 0∼1 MPa and the peak strain error is controlled within 0.1%, so the parameters are reasonable.Figure 7
Conventional triaxial stress-strain curves of two rocks. (a) Mudstone. (b) Sandstone.
(a)(b)Figure 8
Stress-strain curve of rock numerical triaxial test. (a) Mudstone. (b) Sandstone.
(a)(b)


## 3.1.2. Comparative Analysis of Rock Strength Attenuation

The rock compressive strength attenuation value and strength attenuation coefficient reflect the brittleness degree of the rock. In this paper, we cite the mechanical indexes proposed by Peng et al. [15] to characterize the post-peak strength decay behavior of rock samples—rock strength attenuation coefficient:(9)Ds=σ1P−σrσ1P=Δσσ1P,(10)Δσ=σ1P−σr,where Ds is the intensity attenuation coefficient, value range [0, 1]; Δσ is the intensity attenuation value; σ1P is the peak intensity; σr is the residual intensity.During the trial adjustment of mesoparameters, the comparison curves of the relationship between the compressive strength attenuation value and attenuation coefficient of sandstone and mudstone calculated according to equations (9) and (10) are shown in Figures 9(a) and 9(b). It can be found that with the increase of confining pressure, the strength attenuation value of sandstone shows a gradual decreasing trend, which is consistent with the conclusion obtained in [26], and verifies the rationality of this conventional triaxial test. The sandstone attenuation coefficient curve decreased sharply after 2 MPa confining pressure, from 0.67 to 0.37, and then decreased slowly, indicating that the brittleness of sandstone gradually weakened with the increase of confining pressure, and the brittleness weakening ability increased after 2 MPa. In Figures 9(a), the curve of mudstone strength attenuation value and attenuation coefficient changes in a V-shape. When the curve drops to 4 MPa confining pressure, the attenuation value and attenuation coefficient increase simultaneously, rising 7.05 and 0.15, respectively, and the rising curve slows down after 6 MPa. It shows that at a low confining pressure of 2 MPa, the brittleness of the mudstone first weakened due to the large initial porosity, and the ductility increased. At 4 MPa, the brittleness increased and the ductility decreased. When the confining pressure is small, the two rock samples have a large attenuation coefficient and high brittleness. They are basically one-time shear failure, and there are few local cracks, and some tensile cracks will appear; when the confining pressure is large, the attenuation coefficient of the rock samples is relatively reduced, brittleness is weakened, and there are many local cracks. The literature [27] also confirmed this law.Figure 9
Comparison of rock attenuation value and attenuation coefficient. (a) Mudstone. (b) Sandstone.
(a)(b)Comparing the results of the laboratory test and numerical test, it is found that the attenuation value and attenuation coefficient of compressive strength of mudstone are relatively low in the numerical test and those of sandstone are relatively high. The parameters mainly depend on the parallel bonding modulus and the critical damping ratio, while the mudstone modulus is relatively low and the critical damping is relatively low, resulting in low brittleness in the numerical experiment; the sandstone modulus is high and the critical damping is relatively high, resulting in high brittleness.


## 3.1.3. Comparative Analysis of Rock Failure Process

During the numerical triaxial test, PFC3D was used to track the broken contact force chain in the parallel bonding contact model through the DFN joint system, and the broken contact link was expressed in the form of DFN_disk (disk joint), thereby forming a damage crack. In DFN, the internal cracks of rock particles are divided into tensile cracks and shear cracks. During the loading of the model, the type of microcracks can be identified and the number of them can be counted by programming the FISH function. The mesoscopic tensile cracks and shear cracks do not directly reflect the tensile and shear cracks of the macroscopic failure of the rock sample, but in the three-dimensional particle flow model, different mesoscopic crack spatial trajectories can observe the development of the macroscopic damage cracks of the rock. The disc joint reflects the damage characteristics of the mesoscopic particles of the rock. This method is especially important in the analysis of the failure characteristics of the rock sample and the precise trial adjustment of the mesoscopic parameters.The sandstone and mud numerical triaxial test failure shear cracks under different confining pressure conditions and loading methods are shown in Figures10 and 11. In this paper, only the failure test graph of the specimen under the minimum confining pressure of 2 MPa and the maximum confining pressure of 8 MPa is compared with the PFC numerical simulation failure map. From Figure 7, it can be seen that the mudstones mostly show shear cracks on the microjoints, and the fractures between the particles are mostly shear damage. At a confining pressure of 2 MPa, the whole mudstone exhibits a single shear failure, with fewer shear failure cracks, and brittleness characteristics are more obvious, accompanied by a small amount of tensile cracks from the bottom up; as the confining pressure increases, the mudstone exhibits mixed shear failure, multiple shear cracks appear in the sample, the brittleness is weakened, the ductility is enhanced, and the damage cracks increase. When the confining pressure continues to increase to 8 MPa, the microscopic cracks at the bottom and top gradually increase and develop toward the middle, and the shear failure cracks increase and become more obvious. The cracks increase and become more obvious. The development trend of mudstone cracks is similar to the indoor test, and the development characteristics are consistent with the results of the indoor test. It can be seen from Figure 8 that the sandstone shows more tensile cracks on the microjoints, and the tensile fracture failure is mainly among the particles. When the confining pressure is 2 MPa, the mesocracks are generated from the bottom of the sandstone and gradually develop upwards. Due to the low confining pressure, the macroshear crack did not penetrate to the top, and the top microcrack developed slowly; when the confining pressure increased to 8 MPa, the bottom crack gradually expanded to the top, forming a penetrating macroscopic shear crack. As the confining pressure increases, the number gradually increases, forming multiple shear cracks that alternately develop failure characteristics. The development trend of sandstone cracks is also similar to laboratory tests. Through the three-dimensional crack tracing observation of the numerical triaxial test of the particle flow simulation above, it can be seen that the macrofailure characteristics of the two rocks under the PFC numerical test are similar to the indoor test, which further verifies the reliability of the mesoparameter adjustment in this paper.Figure 10
Numerical triaxial failure joints (cracks) of mudstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)Figure 11
Triaxial failure joints (cracks) of sandstone under different confining pressures. (a) Confining pressure of 2 MPa. (b) Confining pressure of 8 MPa.
(a)(b)In summary, the macroscopic mechanical properties of the mesoparameters adjusted by the three-axis numerical test have certain reliability, and the materials are truly usable. The mesoparameters of the obtained sandstone and mudstone materials and the mesoparameters of structural plane are shown in Tables5 and 6, respectively.Table 5
Mesoscopic parameters of sandstone and mudstone materials.
Lithologyd (mm)Ρ (kg·m−3)Ec(GPa)kc∗μλk∗Eb(GPa)σc (MPa)c (MPa)βnβsMudstone0.24310021.00.81.001.02018280.20.2Sandstone0.26300031.00.61.011.03020350.30.6Table 6
Mesoscopic parameters of structural plane.
Structural plane typesj_kn (N·m−1)sj_ks (N·m−1)sj_μsj_c (MPa)sj_φ (°)Horizontal structural plane1010.70.310Vertical structural plane1010.308


## 3.2. Slope Particle Flow Simulation and Failure Mode

### 3.2.1. Slope Failure Process

After the establishment of four slope models with different cavity depths and fissure depths, the earthquake action is applied to them. The failure process of each model after 25 s of earthquake duration is shown in Figure11, and the displacement of the rock block in the vertical z-direction and horizontal y-direction is shown in Figure 12. It can be seen that after 25 s of seismic waves, all the slope models have significant structural plane fracture expansion and penetration phenomena, of which Model 2 and Model 3 both have obvious rock mass collapse and overall failure phenomena, and the main control plane surface failure. The degree gradually increases with the deepening of the rock cavity and the initial fissures; the onset time of failure under the earthquake action of each model is different. 2∼2.5 s began to fail, indicating that as the weathering degree of the rock cavity deepened, the initial crack penetration degree increased, and the faster the rock slope collapse rate under earthquake.Figure 12
The destruction process diagram of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The failure process of Model 1 seismic wave within 25 seconds is shown in Figure11(a). It can be seen that there is no sign of rock collapse during the earthquake period of this model. The displacement values are relatively small. The negative vertical displacements are distributed between 0 and 0.4 m, and the positive horizontal displacements are distributed between 0 and 1.2 m. This indicates that the damage characteristics of Model 1 in the horizontal direction are relatively large. When the earthquake action time is between 0 and 8.5 seconds, there is no vertical displacement change, and each rock block has a horizontal y-direction displacement fluctuation with the shear wave action, and the amplitude of the fluctuation increases with the extension of the earthquake action time. At the same time, there is no obvious deformation and destruction of the two rock mass macrochains, and there is no crack development in the mudstone. At this time, the slope is relatively stable. When the earthquake acts for 9 seconds, the seismic acceleration curve gradually reaches the peak, and the displacement of each rock block in Model 1 has a small frequency fluctuation in the vertical direction. At this time, the contact model of the main control structure near the 1-1 rock block and the 2-2 rock block is locally fractured. Local fracture occurred in the contact model of the main control structure. The first and second layers of rock began to move, and the lateral force of the smooth joint model between the layers increased and gradually destroyed. After 10 seconds, the structure of the no. 1-1 rock block penetrated upward. Therefore, the structural surface of the no. 2-1 rock block fractured, and the crack propagation in the rock cavity intensified. The mudstone particles disintegrated severely. After 11.5 seconds, the vertical displacement curve of Model 1 fell slowly downward. The horizontal displacement to the right gradually increases, and the rock strata intensify. Rock no. 3-2 and no. 1-2 have large vertical displacement due to local fragmentation on the top, which can reach a maximum of 1.35 m. Due to the interaction of adjacent rock blocks at the bottom under earthquake, the structural surface of the no. 4-1 rock block at the top layer was partially broken. After 25 seconds, the upper rock mass was destroyed except for the no. 4-2 rock block, the slope as a whole had a tendency to collapse, and the maximum horizontal displacement of the no. 4-1 rock block was 1.2 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 25% and a mudstone cavity depth of 1.5 m will undergo slight instability and failure under the action of a shear wave of 25 seconds. The destruction mode is fracture-horizontal sliding failure.The failure process of Model 2 seismic wave action within 25 s is shown in Figure12(b). The displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(b) and 13(f). At 8 seconds after the earthquake, due to the deepening of the weathering of the rock cavity, there is no support at the bottom of the no. 1 macrochain. The local fracture of the no. 1-1 rock block preferentially occurs, and the horizontal and vertical displacements increase sharply. After 9 seconds, the internal chain no. 1-2 rock, no. 2-2 rock, and no. 3-2 rock structural surface broke. A partial penetration surface is formed, and then the top of the rock block is shifted outward. After 11 seconds, the no. 1-1 rock collided with the slope for the first time, and the displacement curve suddenly changed. Rock 3-1, rock 3-2, rock 4-1, and rock 4-2 were dumped-rotated and their horizontal and vertical curve slopes increased simultaneously. The maximum horizontal displacement can reach 7.8 m, and the maximum vertical displacement can reach 6.3 m. Rock 2-1 began to collapse due to the earthquake and the squeezing effect of the adjacent rock mass. At the same time, the lower mudstones have cracked and damaged. When the earthquake took place for 19 seconds, the no. 2-1 rock block collided with the no. 1-1 rock block for the first time. Under the action of the earthquake, the top four rock blocks were rotated and destroyed around the rock block no. 2-2. The earthquake lasted from 22 seconds to 25 seconds, and the rock mass collapsed and intensified. The top rock mass separated between layers. Rocks 3-1 and 3-2 simultaneously fell and collided with rock 2-1. Rocks 4-1 and 4-2 rolled outward. Until the end of the calculation, the maximum vertical displacement of rock no. 1-1 is 4.5 m and the maximum horizontal displacement is 2.5 m. The maximum vertical displacement of no. 2-1 rock block is 4 m, and the maximum horizontal displacement is 0.8 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 37.5% and a mudstone cavity depth of 3 m will undergo cracking-caving-dumping failure under the action of a 25-second shear wave.Figure 13
Displacement time history curve of each model. (a) Model 1,y-direction. (b) Model 2, y-direction. (c) Model 3, y-direction. (d) Model 4, y-direction. (e) Model 1, z-direction. (f) Model 2, z-direction. (g) Model 3, z-direction. (h) Model 4, z-direction.
(a)(b)(c)(d)(e)(f)(g)(h)The damage process of Model 3 seismic wave action within 25 s is shown in Figure12(c), and the horizontal y-direction and vertical z-direction displacements are shown in Figures 13(c) and 13(g). After 5 seconds, the crack of no. 1-1 rock’s main control surface quickly initiated and propagated, and local fracture and shedding occurred. After 6 seconds, rock 1-1 collapsed, followed by rocks 4-2 and 3-1. The horizontal displacement of the rock block under the action of earthquake has changed obviously. After 7 seconds, the crack at the no. 2 macrochain at the top of the slope developed downward. As a result, no. 3-2 and no. 2-2 rocks fractured and separated from the parent rock. At this time, the destruction of the no. 3-1 rock block at the top prompted the rapid destruction of the no. 2-1 rock structural plane. The two rocks collapsed at the same time and collided with the no. 1-1 rock. After 10 seconds, as the seismic force continued to increase, the bottom no. 1-2 rock collapsed and broke away from the parent rock. At this time , the internal main control structure has been completely broken and penetrated, and the rock cavity and mudstone particles were destroyed. After 15 seconds, the no. 4-1 rock dumped and destroyed along with the four internal rock blocks and collided with the no. 3-1 rock. Due to the excessive torque of the top rock block, at 20 seconds, the rock joint model between no. 3-2 rock and no. 4-2 rock was destroyed. And no. 4-1 and no. 4-2 rocks rotate and destroy around no. 3-2. At the end of 25 seconds, the top two rock masses turned and destroyed and collided with no. 3-1 rock mass, which has not yet completely collapsed. Throughout the process, the maximum vertical displacement of no. 1–1 rock collapse can reach 4 m, and the horizontal roll-off displacement can reach 1.7 m. The top no. 4-1 and no. 4-2 rocks have a large roll-off tendency. Therefore, the failure mode of a complex horizontal layered rock slope with a fracture penetration rate of 50% and a mudstone cavity depth of 4.5 m under the action of a shear wave of 25 s is caving-dumping-rotation failure. The development of fissures and the weathering of rock cavities are more serious.The weathering degree of Model 4 is the most serious. The destruction process of seismic wave action within 25 s is shown in Figure12(d), and the displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(d) and 13(h). Within 0∼2 seconds, the slope was in a stable state. After 2.5 seconds, the no. 1-1 rock began to fall, and the vertical displacement curve increased significantly, accompanied by horizontal displacement. After 3 seconds, the no. 1-2 rock collapsed and was squeezed by the upper rock mass. After 5 seconds, the no. 1-1 rock rolled down along the surface of the mudstone slope, with a vertical fall displacement of 3.8 m and a horizontal roll-off displacement of 1.45 m. Due to the continued action of earthquakes and interlayer loads, the no. 2-1 rock broke away from the parent rock and began to fall. The earthquake action gradually increased, and the squeezing effect of no. 2-2 rock and no. 1-2 rock on the upper rock mass also continued to increase. After 8 seconds, the slope began to slip overall. After 11 seconds, the no. 3-1 rock began to fall under the action of the earthquake. Within 15∼25 seconds, the joint model of the top rock layer began to break, and the no. 4-1 and no. 4-2 rocks were rotationally damaged. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. At last, Model 4 suffered the most severe damage and the fastest damage speed after 25 seconds of earthquake action. And there is no rolling movement of single rock. The overall trend of slope collapse is more significant, and the failure mode is caving-slip-rotation failure.In summary, under the action of seismic wave in one direction, the damage characteristics of complex horizontal layered rock slopes are affected by the depth of the cavity and the crack penetration rate of the main control surface. The basic failure mode (see Figure12) can be divided into fracture-horizontal slid failure (Model 1), tensile cracking-caving-dumping failure (Model 2), caving-dumping-rotation failure (Model 3), and caving-slip-rotation failure (Model 4). The greater the weathering degree of the rock cavity and fissures, the more obvious the tendency of rock avalanches and the faster the avalanche time. When the rock blocks are collapsing, the horizontal and interlayer adjacent rock blocks will squeeze and collide with each other, which will intensify the fracture of the smooth joint model. The slope has developed from the single rock block’s pull failure to the overall dumping (slippage) rotation failure. In the work of geological disasters, according to the failure modes and displacement characteristics of the slopes obtained in this paper, the reinforcement of Model 1 type slopes can be strengthened, and the rock masses of Model 2∼4 type slopes can be prevented.


### 3.2.2. Coordination Number Analysis

In order to further reveal the seismic failure process and displacement characteristics of the slope, the specific strength attenuation of the contact model during the fracture process and the specific rock collapse time need to be quantitatively described. Therefore, the concept of coordination number is introduced in this paper. Coordination number, that is, the average number of particles in contact is equal to twice the number of particles in contact/number of particles. This parameter comprehensively reflects the physical properties of the particle contact model (friction coefficient, elastic modulus, Poisson’s ratio, etc.). The strength is closely related [28]. In the simulation process of PFC rock fracture failure, it can reflect the fracture time and strength attenuation of the contact force chain of the specific structural surface. In this paper, the particle coordination values in units of rocks in each model are monitored, and the time history curve is plotted (see Figure 14). It can be seen from the figure that the coordination values of Model 1 and Model 4 are distributed between 1–11 and 6–11 respectively, and the coordination values of Model 2 and Model 3 are distributed between 5 and 11, which can be seen as the depth of the cavity and the cracks. As it deepens, the change range of the coordination value gradually decreases. The reason is that the initial stability of the slope gradually decreases. The coordination number attenuation of each model slope in the period of anticavitation collapse is about 2∼15 s, and as the weathering deepens, the coordination number attenuation effect shows a certain advance effect, that is, Models 1∼4 correspond to 8.7∼11.4 s, 7.1∼10.7 s, 5.4∼8.3 s, and 2.3∼11.3 s respectively. In Model 1, the coordination number of rock block 1-2 is the most attenuated by earthquake, from 9.2 s from 10.5 to 1.8. The reason is that during this period, the rock has undergone large local fragmentation and destruction, and some particles break away. For the volume, the attenuation value includes not only the fracture of the structural plane, but also the reduction of the coordination number caused by the fracture of the particles inside the rock. The attenuation values of coordination number of other rock blocks are maintained between 2 and 4. As the No. 4-2 rock mass has not collapsed 25 s after the earthquake, the coordination number attenuation value is at least 2. In each model, as the elevation increases, the coordination number attenuation value decreases, of which the no. 4-1 and no. 4-2 rock bodies placed on the top of Model 4 have the lowest coordination number attenuation value of 0.6 and 1.2. It can be seen that the coordination number has an elevation attenuation effect in dangerous rock, and with the development of the cavity and the depth of the fissure, the elevation attenuation effect gradually weakens. Comparing the coordination number curves of each model, it can be seen that with the deepening of the cavity and cracks, the stability of the main structural surface of the rock slope decreases, the start time of the fluctuation of the coordination number gradually decreases, and the fluctuation range of the coordination value gradually increases. After reaching Model 5, due to the intensified collision motion of rocks, the particle contact behavior between different rocks is more significant, resulting in an increasing trend of the coordination number curve and further increase in fluctuations. It can be seen that the attenuation speed and attenuation degree of the structural plane under the action of an earthquake are determined by the initial stability of the slope (i.e., the depth of the rock cavity and the depth of the crack). The development trend of the coordination number of each model is consistent with the failure process diagram and displacement diagram, which further confirms the specific failure time of each rock block of each model, as shown in Table 7. According to this, it can provide a strong reference for the treatment of horizontal complex layered rock slopes under different working conditions (Models 1 to 4).Figure 14
Time history curve of coordination number of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)Table 7
Fracture failure timetable of each model rock block.
ModelRock block1-1 (s)1-2 (s)2-1 (s)2-2 (s)3-1 (s)3-2 (s)4-1 (s)4-2 (s)Model 18.79.29.68.79.38.911.4>25Model 27.18.410.78.99.18.89.89.7Model 35.48.36.26.35.96.15.85.8Model 42.33.04.64.911.38.09.06.1


### 3.2.3. Energy Dissipation Analysis

The laws of thermodynamics show that energy conversion is an essential feature of the physical change process of matter, and material destruction is a state instability phenomenon driven by energy. The energy dissipation during the failure of the slope is mainly reflected in the failure of the structural plane and the damage and rupture of the intact rock block, reflecting the internal mechanism of the failure evolution. This article assumes that the rock slope system is a closed system that has no heat exchange with the outside world, and it can be obtained according to the first law of thermodynamics [29]:(11)U=Ud+Ue,where U is the total input energy generated by the work of the slope’s body force during the earthquake; Ud is the dissipated energy during the dynamic response of the slope; Ue is the elastic strain energy that can be released inside the slope particles. The total input energy U in the model of this paper is obtained by the following formula:(12)U=Upre+∑NpγpVpΔu,where Upre is the total input energy accumulated in the last calculation time step; γp and Δu are the particle weight and displacement in the current time step, respectively; Np is the number of particles. In the model, the sandstone and mudstone particle contact model is a parallel bonding model. The elastic strain energy Ue consists of the particle strain energy Uc and the parallel bond strain energy Upb, namely,(13)Ue=Uc+Upb,Uc=12∑NcFin2/kn+Fis2/ks,Upb=12∑NpbF¯in2/Ak¯n+F¯is2/Ak¯s+M¯2/Ik¯n,where Nc and Npb are the number of contact force chains and the number of parallel bonding; i represents the i-th contact chain; Fin and Fis are the normal contact force and tangential contact force; F¯in, F¯is, and M¯ are the normal parallel bonding force, tangential parallel bonding force, and parallel bonding moment; kn and ks are the contact normal stiffness and contact tangential stiffness; k¯n, k¯s, and k¯n are, respectively, the stiffness corresponding to each parallel bonding force; A and I are the cross-sectional area and moment of inertia of the parallel bonding.The energy dissipationUd of the final model can be obtained by substituting various equations into the solution. According to the first law of thermodynamics and the particle flow comes with history energy detection system, the instantaneous kinetic energy of the slope can be obtained.The kinetic energy response of each model under the action of earthquake is shown in Figure15. It can be seen that from Model 1 to Model 4, the maximum kinetic energy decreased from 155 MJ to 72 MJ, indicating that as the development of the cavity and cracks deepens, the energy generated by the failure of the slope gradually decreases. This is due to the deepening of the cavity and the upper and lower sandstones. The friction with the mudstone particles is reduced, and the degree of mudstone particle disintegration is reduced; the cracks are deepened, and the flat bond chain of the structural surface is reduced, that is, the energy reflection is reduced. In Models 1 to 3, the kinetic energy fluctuations gradually increase with the duration of the earthquake, and the energy gradually increases with the severity of the slope damage. Among them, Model 1 reached the first peak of 46 MJ at 6 s, indicating that the contact force chain began to break, and then the peak fluctuated between 75 and 90 MJ during the caving process, which was caused by the collision of the rocks. As a result, the maximum peak reached 155 MJ at 13 s; Model 2 produced a peak of 45 MJ kinetic energy at 2 s due to local particle collisions and subsequently produced peaks of 85 MJ and 89 MJ at 13 s and 15 s; similarly, in Model 3 at 15 s, the instantaneous energy of 104 MJ was generated due to collision of rock particles. In Model 4, the kinetic energy response increases rapidly to a maximum value of 73 MJ with the duration of the earthquake and then decreases to fluctuate between 0 and 43 MJ. It can be seen that due to the strong stability of Models 1 to 3, the dangerous rock gradually collapses and breaks down over time, the slope shows progressive failure, and the kinetic energy response gradually increases. Because of the weak stability of Model 4, the whole slope is in the early stage When the destruction is completed in a short time, the kinetic energy response reaches the maximum value at the initial stage, and then it collapses over time, which is in good synchronization with the failure mode of the slope.Figure 15
Variation of kinetic energy response of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The energy dissipation evolution process of each model is shown in Figure16. It can be seen that Model 1 and Model 4 produce the most total energy. The maximum energy dissipation of Models 1 to 4 is 105 MJ, 76 MJ, 63 MJ, and 82 MJ, respectively. Model 1 has the largest total energy dissipated due to more structural surface contact, which means the dissipated energy is the largest. The interaction between the blocks is large, so the dissipation energy is higher than that of Model 2 and Model 3. The energy dissipation of each slope model along with the development of earthquake damage and destruction can be divided into two stages, and the whole has the characteristics of first followed by slow. The first stage is the rock block fracture and collision stage, which corresponds to a period of 2.5∼15 s, the period when the dangerous rock blocks break away from the parent rock and collides with each other. At this stage, the change of dissipated energy of Model 1 is 0∼100 MJ, accounting for 95% of the total energy; the change of dissipated energy of Model 2 is 0∼68 MJ, accounting for 89% of the total energy; that of Model 3 is 0∼56 MJ, accounting for 89% of the total energy; that of Model 4 is 74 MJ, accounting for 90% of the total energy. At this stage, the dissipated energy curve develops faster and faster with the duration of the earthquake. Corresponding to the process of “cracking-caving” and “dumping-collision” in the failure mode, it is first “cracking-caving” and then “dump-collision.” And the energy dissipation caused by the dumping collision is greater than the energy generated by the fracture of the structural plane, and it also reflects the progressive cumulative energy dissipation mechanism of the development of dangerous rock failure. The second stage is the stage of the overall collapse of the slope. In this stage, the rock mass has completed local fracture and gradually fell to the foot of the slope. At this stage, the energy dissipation change of Model 1 is 100–105 MJ, accounting for 5% of the total energy. That of Model 2 is 68–76 MJ, of Model 3 is 56–63 MJ, and of Model 4 is 74–82 MJ. The change with time is slower and the change value is smaller, corresponding to the process of “overall dumping” and “overall seat sliding” in the failure mode. Comprehensive analysis shows that the trend of energy dissipation is closely related to the failure mode. The degree of overall energy dissipation affects the cavity and the depth of the fissures. There is a significant difference between the energy consumption of the first stage and the second stage. The energy consumption of horizontal complex layered rock slopes under earthquake action is more concentrated in the first stage of rock fracture collision.Figure 16
Comparison of energy dissipation of different models.


## 3.2.1. Slope Failure Process

After the establishment of four slope models with different cavity depths and fissure depths, the earthquake action is applied to them. The failure process of each model after 25 s of earthquake duration is shown in Figure11, and the displacement of the rock block in the vertical z-direction and horizontal y-direction is shown in Figure 12. It can be seen that after 25 s of seismic waves, all the slope models have significant structural plane fracture expansion and penetration phenomena, of which Model 2 and Model 3 both have obvious rock mass collapse and overall failure phenomena, and the main control plane surface failure. The degree gradually increases with the deepening of the rock cavity and the initial fissures; the onset time of failure under the earthquake action of each model is different. 2∼2.5 s began to fail, indicating that as the weathering degree of the rock cavity deepened, the initial crack penetration degree increased, and the faster the rock slope collapse rate under earthquake.Figure 12
The destruction process diagram of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The failure process of Model 1 seismic wave within 25 seconds is shown in Figure11(a). It can be seen that there is no sign of rock collapse during the earthquake period of this model. The displacement values are relatively small. The negative vertical displacements are distributed between 0 and 0.4 m, and the positive horizontal displacements are distributed between 0 and 1.2 m. This indicates that the damage characteristics of Model 1 in the horizontal direction are relatively large. When the earthquake action time is between 0 and 8.5 seconds, there is no vertical displacement change, and each rock block has a horizontal y-direction displacement fluctuation with the shear wave action, and the amplitude of the fluctuation increases with the extension of the earthquake action time. At the same time, there is no obvious deformation and destruction of the two rock mass macrochains, and there is no crack development in the mudstone. At this time, the slope is relatively stable. When the earthquake acts for 9 seconds, the seismic acceleration curve gradually reaches the peak, and the displacement of each rock block in Model 1 has a small frequency fluctuation in the vertical direction. At this time, the contact model of the main control structure near the 1-1 rock block and the 2-2 rock block is locally fractured. Local fracture occurred in the contact model of the main control structure. The first and second layers of rock began to move, and the lateral force of the smooth joint model between the layers increased and gradually destroyed. After 10 seconds, the structure of the no. 1-1 rock block penetrated upward. Therefore, the structural surface of the no. 2-1 rock block fractured, and the crack propagation in the rock cavity intensified. The mudstone particles disintegrated severely. After 11.5 seconds, the vertical displacement curve of Model 1 fell slowly downward. The horizontal displacement to the right gradually increases, and the rock strata intensify. Rock no. 3-2 and no. 1-2 have large vertical displacement due to local fragmentation on the top, which can reach a maximum of 1.35 m. Due to the interaction of adjacent rock blocks at the bottom under earthquake, the structural surface of the no. 4-1 rock block at the top layer was partially broken. After 25 seconds, the upper rock mass was destroyed except for the no. 4-2 rock block, the slope as a whole had a tendency to collapse, and the maximum horizontal displacement of the no. 4-1 rock block was 1.2 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 25% and a mudstone cavity depth of 1.5 m will undergo slight instability and failure under the action of a shear wave of 25 seconds. The destruction mode is fracture-horizontal sliding failure.The failure process of Model 2 seismic wave action within 25 s is shown in Figure12(b). The displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(b) and 13(f). At 8 seconds after the earthquake, due to the deepening of the weathering of the rock cavity, there is no support at the bottom of the no. 1 macrochain. The local fracture of the no. 1-1 rock block preferentially occurs, and the horizontal and vertical displacements increase sharply. After 9 seconds, the internal chain no. 1-2 rock, no. 2-2 rock, and no. 3-2 rock structural surface broke. A partial penetration surface is formed, and then the top of the rock block is shifted outward. After 11 seconds, the no. 1-1 rock collided with the slope for the first time, and the displacement curve suddenly changed. Rock 3-1, rock 3-2, rock 4-1, and rock 4-2 were dumped-rotated and their horizontal and vertical curve slopes increased simultaneously. The maximum horizontal displacement can reach 7.8 m, and the maximum vertical displacement can reach 6.3 m. Rock 2-1 began to collapse due to the earthquake and the squeezing effect of the adjacent rock mass. At the same time, the lower mudstones have cracked and damaged. When the earthquake took place for 19 seconds, the no. 2-1 rock block collided with the no. 1-1 rock block for the first time. Under the action of the earthquake, the top four rock blocks were rotated and destroyed around the rock block no. 2-2. The earthquake lasted from 22 seconds to 25 seconds, and the rock mass collapsed and intensified. The top rock mass separated between layers. Rocks 3-1 and 3-2 simultaneously fell and collided with rock 2-1. Rocks 4-1 and 4-2 rolled outward. Until the end of the calculation, the maximum vertical displacement of rock no. 1-1 is 4.5 m and the maximum horizontal displacement is 2.5 m. The maximum vertical displacement of no. 2-1 rock block is 4 m, and the maximum horizontal displacement is 0.8 m. In summary, a complex horizontal layered rock slope with a fracture penetration rate of 37.5% and a mudstone cavity depth of 3 m will undergo cracking-caving-dumping failure under the action of a 25-second shear wave.Figure 13
Displacement time history curve of each model. (a) Model 1,y-direction. (b) Model 2, y-direction. (c) Model 3, y-direction. (d) Model 4, y-direction. (e) Model 1, z-direction. (f) Model 2, z-direction. (g) Model 3, z-direction. (h) Model 4, z-direction.
(a)(b)(c)(d)(e)(f)(g)(h)The damage process of Model 3 seismic wave action within 25 s is shown in Figure12(c), and the horizontal y-direction and vertical z-direction displacements are shown in Figures 13(c) and 13(g). After 5 seconds, the crack of no. 1-1 rock’s main control surface quickly initiated and propagated, and local fracture and shedding occurred. After 6 seconds, rock 1-1 collapsed, followed by rocks 4-2 and 3-1. The horizontal displacement of the rock block under the action of earthquake has changed obviously. After 7 seconds, the crack at the no. 2 macrochain at the top of the slope developed downward. As a result, no. 3-2 and no. 2-2 rocks fractured and separated from the parent rock. At this time, the destruction of the no. 3-1 rock block at the top prompted the rapid destruction of the no. 2-1 rock structural plane. The two rocks collapsed at the same time and collided with the no. 1-1 rock. After 10 seconds, as the seismic force continued to increase, the bottom no. 1-2 rock collapsed and broke away from the parent rock. At this time , the internal main control structure has been completely broken and penetrated, and the rock cavity and mudstone particles were destroyed. After 15 seconds, the no. 4-1 rock dumped and destroyed along with the four internal rock blocks and collided with the no. 3-1 rock. Due to the excessive torque of the top rock block, at 20 seconds, the rock joint model between no. 3-2 rock and no. 4-2 rock was destroyed. And no. 4-1 and no. 4-2 rocks rotate and destroy around no. 3-2. At the end of 25 seconds, the top two rock masses turned and destroyed and collided with no. 3-1 rock mass, which has not yet completely collapsed. Throughout the process, the maximum vertical displacement of no. 1–1 rock collapse can reach 4 m, and the horizontal roll-off displacement can reach 1.7 m. The top no. 4-1 and no. 4-2 rocks have a large roll-off tendency. Therefore, the failure mode of a complex horizontal layered rock slope with a fracture penetration rate of 50% and a mudstone cavity depth of 4.5 m under the action of a shear wave of 25 s is caving-dumping-rotation failure. The development of fissures and the weathering of rock cavities are more serious.The weathering degree of Model 4 is the most serious. The destruction process of seismic wave action within 25 s is shown in Figure12(d), and the displacements in the horizontal y-direction and vertical z-direction are shown in Figures 13(d) and 13(h). Within 0∼2 seconds, the slope was in a stable state. After 2.5 seconds, the no. 1-1 rock began to fall, and the vertical displacement curve increased significantly, accompanied by horizontal displacement. After 3 seconds, the no. 1-2 rock collapsed and was squeezed by the upper rock mass. After 5 seconds, the no. 1-1 rock rolled down along the surface of the mudstone slope, with a vertical fall displacement of 3.8 m and a horizontal roll-off displacement of 1.45 m. Due to the continued action of earthquakes and interlayer loads, the no. 2-1 rock broke away from the parent rock and began to fall. The earthquake action gradually increased, and the squeezing effect of no. 2-2 rock and no. 1-2 rock on the upper rock mass also continued to increase. After 8 seconds, the slope began to slip overall. After 11 seconds, the no. 3-1 rock began to fall under the action of the earthquake. Within 15∼25 seconds, the joint model of the top rock layer began to break, and the no. 4-1 and no. 4-2 rocks were rotationally damaged. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. The downward pressure caused the 1-2, 2-2, and 3-2 rocks to undergo an integrated seat-slide failure. At last, Model 4 suffered the most severe damage and the fastest damage speed after 25 seconds of earthquake action. And there is no rolling movement of single rock. The overall trend of slope collapse is more significant, and the failure mode is caving-slip-rotation failure.In summary, under the action of seismic wave in one direction, the damage characteristics of complex horizontal layered rock slopes are affected by the depth of the cavity and the crack penetration rate of the main control surface. The basic failure mode (see Figure12) can be divided into fracture-horizontal slid failure (Model 1), tensile cracking-caving-dumping failure (Model 2), caving-dumping-rotation failure (Model 3), and caving-slip-rotation failure (Model 4). The greater the weathering degree of the rock cavity and fissures, the more obvious the tendency of rock avalanches and the faster the avalanche time. When the rock blocks are collapsing, the horizontal and interlayer adjacent rock blocks will squeeze and collide with each other, which will intensify the fracture of the smooth joint model. The slope has developed from the single rock block’s pull failure to the overall dumping (slippage) rotation failure. In the work of geological disasters, according to the failure modes and displacement characteristics of the slopes obtained in this paper, the reinforcement of Model 1 type slopes can be strengthened, and the rock masses of Model 2∼4 type slopes can be prevented.


## 3.2.2. Coordination Number Analysis

In order to further reveal the seismic failure process and displacement characteristics of the slope, the specific strength attenuation of the contact model during the fracture process and the specific rock collapse time need to be quantitatively described. Therefore, the concept of coordination number is introduced in this paper. Coordination number, that is, the average number of particles in contact is equal to twice the number of particles in contact/number of particles. This parameter comprehensively reflects the physical properties of the particle contact model (friction coefficient, elastic modulus, Poisson’s ratio, etc.). The strength is closely related [28]. In the simulation process of PFC rock fracture failure, it can reflect the fracture time and strength attenuation of the contact force chain of the specific structural surface. In this paper, the particle coordination values in units of rocks in each model are monitored, and the time history curve is plotted (see Figure 14). It can be seen from the figure that the coordination values of Model 1 and Model 4 are distributed between 1–11 and 6–11 respectively, and the coordination values of Model 2 and Model 3 are distributed between 5 and 11, which can be seen as the depth of the cavity and the cracks. As it deepens, the change range of the coordination value gradually decreases. The reason is that the initial stability of the slope gradually decreases. The coordination number attenuation of each model slope in the period of anticavitation collapse is about 2∼15 s, and as the weathering deepens, the coordination number attenuation effect shows a certain advance effect, that is, Models 1∼4 correspond to 8.7∼11.4 s, 7.1∼10.7 s, 5.4∼8.3 s, and 2.3∼11.3 s respectively. In Model 1, the coordination number of rock block 1-2 is the most attenuated by earthquake, from 9.2 s from 10.5 to 1.8. The reason is that during this period, the rock has undergone large local fragmentation and destruction, and some particles break away. For the volume, the attenuation value includes not only the fracture of the structural plane, but also the reduction of the coordination number caused by the fracture of the particles inside the rock. The attenuation values of coordination number of other rock blocks are maintained between 2 and 4. As the No. 4-2 rock mass has not collapsed 25 s after the earthquake, the coordination number attenuation value is at least 2. In each model, as the elevation increases, the coordination number attenuation value decreases, of which the no. 4-1 and no. 4-2 rock bodies placed on the top of Model 4 have the lowest coordination number attenuation value of 0.6 and 1.2. It can be seen that the coordination number has an elevation attenuation effect in dangerous rock, and with the development of the cavity and the depth of the fissure, the elevation attenuation effect gradually weakens. Comparing the coordination number curves of each model, it can be seen that with the deepening of the cavity and cracks, the stability of the main structural surface of the rock slope decreases, the start time of the fluctuation of the coordination number gradually decreases, and the fluctuation range of the coordination value gradually increases. After reaching Model 5, due to the intensified collision motion of rocks, the particle contact behavior between different rocks is more significant, resulting in an increasing trend of the coordination number curve and further increase in fluctuations. It can be seen that the attenuation speed and attenuation degree of the structural plane under the action of an earthquake are determined by the initial stability of the slope (i.e., the depth of the rock cavity and the depth of the crack). The development trend of the coordination number of each model is consistent with the failure process diagram and displacement diagram, which further confirms the specific failure time of each rock block of each model, as shown in Table 7. According to this, it can provide a strong reference for the treatment of horizontal complex layered rock slopes under different working conditions (Models 1 to 4).Figure 14
Time history curve of coordination number of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)Table 7
Fracture failure timetable of each model rock block.
ModelRock block1-1 (s)1-2 (s)2-1 (s)2-2 (s)3-1 (s)3-2 (s)4-1 (s)4-2 (s)Model 18.79.29.68.79.38.911.4>25Model 27.18.410.78.99.18.89.89.7Model 35.48.36.26.35.96.15.85.8Model 42.33.04.64.911.38.09.06.1


## 3.2.3. Energy Dissipation Analysis

The laws of thermodynamics show that energy conversion is an essential feature of the physical change process of matter, and material destruction is a state instability phenomenon driven by energy. The energy dissipation during the failure of the slope is mainly reflected in the failure of the structural plane and the damage and rupture of the intact rock block, reflecting the internal mechanism of the failure evolution. This article assumes that the rock slope system is a closed system that has no heat exchange with the outside world, and it can be obtained according to the first law of thermodynamics [29]:(11)U=Ud+Ue,where U is the total input energy generated by the work of the slope’s body force during the earthquake; Ud is the dissipated energy during the dynamic response of the slope; Ue is the elastic strain energy that can be released inside the slope particles. The total input energy U in the model of this paper is obtained by the following formula:(12)U=Upre+∑NpγpVpΔu,where Upre is the total input energy accumulated in the last calculation time step; γp and Δu are the particle weight and displacement in the current time step, respectively; Np is the number of particles. In the model, the sandstone and mudstone particle contact model is a parallel bonding model. The elastic strain energy Ue consists of the particle strain energy Uc and the parallel bond strain energy Upb, namely,(13)Ue=Uc+Upb,Uc=12∑NcFin2/kn+Fis2/ks,Upb=12∑NpbF¯in2/Ak¯n+F¯is2/Ak¯s+M¯2/Ik¯n,where Nc and Npb are the number of contact force chains and the number of parallel bonding; i represents the i-th contact chain; Fin and Fis are the normal contact force and tangential contact force; F¯in, F¯is, and M¯ are the normal parallel bonding force, tangential parallel bonding force, and parallel bonding moment; kn and ks are the contact normal stiffness and contact tangential stiffness; k¯n, k¯s, and k¯n are, respectively, the stiffness corresponding to each parallel bonding force; A and I are the cross-sectional area and moment of inertia of the parallel bonding.The energy dissipationUd of the final model can be obtained by substituting various equations into the solution. According to the first law of thermodynamics and the particle flow comes with history energy detection system, the instantaneous kinetic energy of the slope can be obtained.The kinetic energy response of each model under the action of earthquake is shown in Figure15. It can be seen that from Model 1 to Model 4, the maximum kinetic energy decreased from 155 MJ to 72 MJ, indicating that as the development of the cavity and cracks deepens, the energy generated by the failure of the slope gradually decreases. This is due to the deepening of the cavity and the upper and lower sandstones. The friction with the mudstone particles is reduced, and the degree of mudstone particle disintegration is reduced; the cracks are deepened, and the flat bond chain of the structural surface is reduced, that is, the energy reflection is reduced. In Models 1 to 3, the kinetic energy fluctuations gradually increase with the duration of the earthquake, and the energy gradually increases with the severity of the slope damage. Among them, Model 1 reached the first peak of 46 MJ at 6 s, indicating that the contact force chain began to break, and then the peak fluctuated between 75 and 90 MJ during the caving process, which was caused by the collision of the rocks. As a result, the maximum peak reached 155 MJ at 13 s; Model 2 produced a peak of 45 MJ kinetic energy at 2 s due to local particle collisions and subsequently produced peaks of 85 MJ and 89 MJ at 13 s and 15 s; similarly, in Model 3 at 15 s, the instantaneous energy of 104 MJ was generated due to collision of rock particles. In Model 4, the kinetic energy response increases rapidly to a maximum value of 73 MJ with the duration of the earthquake and then decreases to fluctuate between 0 and 43 MJ. It can be seen that due to the strong stability of Models 1 to 3, the dangerous rock gradually collapses and breaks down over time, the slope shows progressive failure, and the kinetic energy response gradually increases. Because of the weak stability of Model 4, the whole slope is in the early stage When the destruction is completed in a short time, the kinetic energy response reaches the maximum value at the initial stage, and then it collapses over time, which is in good synchronization with the failure mode of the slope.Figure 15
Variation of kinetic energy response of each model. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4.
(a)(b)(c)(d)The energy dissipation evolution process of each model is shown in Figure16. It can be seen that Model 1 and Model 4 produce the most total energy. The maximum energy dissipation of Models 1 to 4 is 105 MJ, 76 MJ, 63 MJ, and 82 MJ, respectively. Model 1 has the largest total energy dissipated due to more structural surface contact, which means the dissipated energy is the largest. The interaction between the blocks is large, so the dissipation energy is higher than that of Model 2 and Model 3. The energy dissipation of each slope model along with the development of earthquake damage and destruction can be divided into two stages, and the whole has the characteristics of first followed by slow. The first stage is the rock block fracture and collision stage, which corresponds to a period of 2.5∼15 s, the period when the dangerous rock blocks break away from the parent rock and collides with each other. At this stage, the change of dissipated energy of Model 1 is 0∼100 MJ, accounting for 95% of the total energy; the change of dissipated energy of Model 2 is 0∼68 MJ, accounting for 89% of the total energy; that of Model 3 is 0∼56 MJ, accounting for 89% of the total energy; that of Model 4 is 74 MJ, accounting for 90% of the total energy. At this stage, the dissipated energy curve develops faster and faster with the duration of the earthquake. Corresponding to the process of “cracking-caving” and “dumping-collision” in the failure mode, it is first “cracking-caving” and then “dump-collision.” And the energy dissipation caused by the dumping collision is greater than the energy generated by the fracture of the structural plane, and it also reflects the progressive cumulative energy dissipation mechanism of the development of dangerous rock failure. The second stage is the stage of the overall collapse of the slope. In this stage, the rock mass has completed local fracture and gradually fell to the foot of the slope. At this stage, the energy dissipation change of Model 1 is 100–105 MJ, accounting for 5% of the total energy. That of Model 2 is 68–76 MJ, of Model 3 is 56–63 MJ, and of Model 4 is 74–82 MJ. The change with time is slower and the change value is smaller, corresponding to the process of “overall dumping” and “overall seat sliding” in the failure mode. Comprehensive analysis shows that the trend of energy dissipation is closely related to the failure mode. The degree of overall energy dissipation affects the cavity and the depth of the fissures. There is a significant difference between the energy consumption of the first stage and the second stage. The energy consumption of horizontal complex layered rock slopes under earthquake action is more concentrated in the first stage of rock fracture collision.Figure 16
Comparison of energy dissipation of different models.


## 4. Conclusion

In order to study the failure mechanism of horizontal complex layered rock slope under seismic-weathering coupling, triaxial test of sandstone and mudstone, and PFC numerical simulation of horizontal complex layered rock slope were carried out in this paper. The conclusions obtained with practical application value are as follows:(1)
In the conventional triaxial test, the brittleness of sandstone gradually weakened with the increase of confining pressure, and the brittleness weakened more sharply after the confining pressure of 2 MPa; the characteristic value of the brittleness attenuation of the mudstone showed a V-shaped change. At a low confining pressure of 2 MPa, the brittleness of mudstone is preferentially weakened and the ductility characteristics are enhanced. At a confining pressure of 4 MPa, the brittleness of mudstone is enhanced and ductility decreases. During the PFC numerical triaxial servo test, it was found that the mudstone particles and sandstone particles exhibited shear failure and tensile failure, respectively. The failure characteristics of the two rocks in the PFC numerical test are similar to those in the laboratory test.(2)
Under the action of the earthquake in one direction, the destruction characteristics of complex horizontal layered rock slopes change with the depth of the cavity and the deepening of the penetration rate of the main control structural surface. The basic failure modes are divided into tensile cracking-horizontal sliding failure (Model 1), tensile cracking-caving-dumping failure (Model 2), caving-dumping-rotating failure (Model 3), and caving-slip-rotation failure (Model 4). By analyzing they and z displacements of the horizontal complex layered rock slope after 25 s of earthquake action, we can see that Model 2 has the largest horizontal and vertical final displacement, the maximum vertical displacement can reach 6.2 m, and the maximum positive horizontal displacement is 7.8 m. And the maximum horizontal final displacements of models 1∼4 are all distributed on the top rock block.(3)
The coordination number of the slope rock blocks of each model increases with the degree of weathering during the earthquake action period of 2∼15 s, showing the attenuation advance effect and the elevation attenuation effect. Among them, the coordination number of the rock block of Model 1 is the most attenuated during the seismic action from 8.7 to 11.4 s, and the coordination number of the rock block 1-2 is attenuated from 10.5 to 1.8 from 9.2 s to the end of the seismic action, and the fragmentation and damage of the rock block are serious; the coordination number of block 4 in Model 4 has the largest attenuation in 2.3∼11.3 s of earthquake action, and the coordination number of top 4-1 and 4-2 rock blocks is reduced to 0.6 and 1.2, respectively, after the end of earthquake action. It shows that the decay rate of the coordination number of the slope rock block under the action of earthquake is determined by the depth of the slope cavity and the depth of the fissure, and as the depth of the cavity and fissure deepens, the effect of coordination number elevation attenuation gradually weakens.(4)
By analyzing the energy dissipation of complex layered rock slopes, it is found that the total energy dissipation of the slope under earthquake action can be divided into two stages: fracture collision and overall collapse, and the overall energy consumption curve changes first with time. After the slow characteristic, the energy dissipation is more concentrated in the rock fracture and collision stage, accounting for 89% to 95% of the total energy. Moreover, with the increasing weathering degree of the slope, the peak value of the overall kinetic energy of Models 1 to 4 decreased from 155 MJ to 72 MJ, indicating that the instantaneous energy generated by the failure of the slope was negatively related to the degree of slope weathering.



---

*Source: 8825440-2020-08-04.xml*