Abstract
Rock slope instability by earthquakes results in substantial economic and property losses. The calculation method of interlayer load and stability coefficient of horizontal complex layered rock slopes in high-intensity areas is established from material mechanics, fracture mechanics, and dynamics. The stability of horizontal layered dangerous rock is calculated after combining it with PFC simulation technology to verify the rationality of the calculation in the Wenchuan area of Sichuan Province. The dynamic response characteristics of dangerous rocks under different weathering degrees are also analyzed. The results show that both methods have an excellent early warning effect on earthquake dangerous rocks. Among the PGA amplification factors, Model 1 has a relatively uniform distribution, Model 2 has a zigzag distribution, Models 3 and 4 have a “U”-shaped distribution, and the most severe acceleration dynamic responses are 4-1 and 4-2 rock blocks. The dynamic acceleration response of mudstone is affected by the crack propagation process of the upper sandstone and exhibits a particular elevation amplification effect. The peak stress gradually decreases with the increase in weathering and elevation. The stress change of the inner chain No. 2 in the horizontal x and y directions is severe, and the stress response of the outer chain No. 1 in the vertical z-direction is severe. It recommends that earthquake disaster protection projects should pay attention to the impact of low-frequency (0–10 Hz) and high-frequency (250 Hz) earthquakes on slope stability.
1. Introduction
According to the survey, sedimentary rocks with layered structures account for 66.7% of the total land area and around 77% in China. Meanwhile, many metamorphic rocks also have layered structure characteristics [1]. In the layered rock mass engineering, the horizontal layered dangerous rock earthquake disaster prevention and treatment is a project often encountered in the construction of highways, railways, water conservancy, and other infrastructures. Earthquake disasters frequently occur in high-intensity areas in seismic zones. There are numerous examples of earthquake-induced instability of dangerous rocks. It often causes substantial economic and property losses. The horizontal layered dangerous rock is the primary carrier of the earthquake-type slope instability disaster among them. Different structural planes cut the rock mass when the main control structural plane at the rear of the dangerous rock block gradually penetrates until it breaks under the action of multiple factors, and the slope will be unstable and destroyed. The existing studies mostly focus on the instability and failure of the rock slope in the nonseismic area under the action of gravity; nevertheless, the dynamic response of the horizontal layered rock slope under the earthquake action needs to be solved urgently. Utah is a rock cavity horizontal layered dangerous rock under the action of earthquakes. The study of its collapse deduction process, failure mechanism, and dynamic response characteristics under different weathering degrees has gradually become one of the research hotspots and technical problems in this field [2, 3].
Chen [4] deduced the chain development process of dangerous rocks from the geomorphological perspective, revealed the influence of rock cave weathering on the stability of dangerous rock slopes, and proposed the chain law of dangerous rock development in the Three Gorges Reservoir area. Tang et al. [5, 6] discussed the sequence of dangerous rock caving on the cliff with a weak base and established a calculation method for the time of dangerous rock caving based on fracture mechanics and damage mechanics. In contrast, Wang et al. [7] established a calculation method for interlayer load and stability coefficient of a complex gently inclined rock slope and revealed its caving sequence with finite element method using fracture mechanics and material mechanics. Zhou et al. [8] used QUIVER code to realize the response analysis of the gentle slope of an earthquake and revealed the slope response and sliding displacement law under the action of an earthquake. Huang et al. [9] used the discontinuous deformation analysis (DDA) method to simulate the movement characteristics of dangerous rock slopes under different seismic conditions, revealing the influence of seismic loads on the movement characteristics of dangerous rock mass collapse blocks. In the discrete element numerical analysis, Varnes classification is a new classification of landslide type. Hungr et al. [10] made a complete summary based on Varnes classification and generalized the slip, dump, fall, flow, expansion, and slope six significant categories of deformation. Dowing and Gilbert [11, 12] used the coupling system of discrete element and finite element to analyze the dynamic interaction between the underground chamber and surrounding rock medium. The discrete element simulated the jointed rock mass near the chamber, while the finite element analyzed the complete rock in the far-field. Also, Meguro and Hakuno [13] used the discrete element model to analyze the seismic failure process of concrete structures and made particular progress. With the rapid development of micromechanics and computer hardware, some numerical simulation methods for investigating the microstructure and properties of materials have emerged in recent years. For instance, the particle flow code (PFC) [14] derived from discrete element theory is a numerical analysis method in a thermotropic state. After the proposal by Potyondy and Cundall in 2004 of a bonded particle model (BPM) suitable for rock properties [15], many scholars have conducted particle flow related to rock fracture and crack propagation. Numerical test simulation research has been universally recognized by the international rock mechanics community [16–23]. Chen et al. [24] realized the particle flow simulation of the fracture propagation of the single dangerous rock main control structure surface and the stability evaluation of the dangerous rock. Tang et al. [25] used a two-dimensional discrete element simulation to study the formation mechanism of Caoling landslide in 1941.
Studies mentioned earlier consider the effects of a single factor and cannot directly reveal the failure process of the slope. As a result, not only the applicability is low but also the understanding of the slope seismic dynamic response law and its failure process is insufficient. In this paper, we take the dangerous rock slope in Wenchuan, Sichuan, as the research object and derive its dynamic stability coefficient and failure law. Moreover, we use discrete element PFC3D software to simulate its dynamic destruction process. At the same time, the dynamic response law under the action of an earthquake is revealed. Finally, the time-frequency characteristics of dangerous rocks are analyzed by HHT. This study has essential scientific guiding significance and practical value for the dynamic stability evaluation and prediction of such dangerous rock slopes under the action of an earthquake.
2. Materials and Methods
2.1. Artificial Synthetic Rock Mass Technology and Model Establishment
Complex geological bodies composed of structural plane networks and rock blocks often exhibit features such as nonuniformity, discontinuity, and anisotropy due to long-term complex geological structures. Structural planes are prone to damage because of the poor mechanical properties on the soft surface of the rock mass [26]. Due to the proposed technology of synthetic rock mass, it is possibly better to simulate rock mass and its mechanical properties in PFC. In PFC 6.0, the artificial synthetic rock mass is mainly composed of two parts: a parallel bonding model between particles and a discrete fracture network [27], as shown in Figure 1. In the dangerous rock, due to weathering, the central control structure plane penetrates to form a fissure, and the interlayer structure plane constitutes the inducing factor of the dangerous rock caving. In order to avoid the influence of the complex structural surface network inside the rock block on the main control structural surface, this paper makes the following assumptions for the PFC simulation of dangerous rocks:(1)Each dangerous rock block in the model does not contain a structural plane network, and each dangerous rock block needs to be set as a whole.(2)Arrange the discrete fracture network on each central control structure surface and interlayer structure surface.(3)The contact characteristics between particles are determined by the contact particles and the constitutive contact model.
The crack propagation of the main control structure and the weathering of the base rock cavity are the key factors affecting the stability of the dangerous rock. This paper uses the formula for calculating the stability of dangerous rock proposed by Chen [28, 29] and divides the slope into four types according to the degree of weathering of the rock cavity and the degree of penetration of the main control structural surface. The layout of the slope model and the particle flow model is shown in Figures 2(a) and 2(b), respectively. Table 1 shows the weathering of rock cavities and fissures in each model.
(a)
(b)
There are two or more joints in the horizontal layered dangerous rock, which are formed by the interaction of multiple microchains and multiple macrochains. There are mutual influences between each microchain and between each macrochain. Under the cutting of the main control structure and the layer, the slope is broken down into many rock blocks, contributing to failure mode differences between horizontal layered dangerous rock and single dangerous rock. In Figure 3, the vertical rows of rock blocks form a microchain, numbered from the outside to the inside. This paper proposes four models to apply seismic waves to dangerous rock slopes, respectively. According to the literature [28, 29], Models 1–3 are loaded by transverse waves, while Model 4 is loaded by longitudinal waves. When analyzing the dynamic response of a slope under earthquake action, it is necessary to input the seismic wave acceleration time-history curve at the bottom of the model. A viscous boundary is required to prevent seismic waves from reflecting at the model boundary. In PFC3D, it is realized by setting dampers in the standard and tangential directions of the two sides and the bottom boundary of the model. The specific seismic waves are shown in Figure 3. During the loading process, seismic acceleration is transformed into boundary acceleration, which drives the model to vibrate.
(a)
(b)
In PFC data monitoring and analysis, it is necessary to use the history command to track the data of the researched object, and the tracked content includes data such as acceleration and stress state. This article uses the PGA amplification factor (peak ground acceleration) to analyze the peak acceleration. The ratio of the peak acceleration of the rock block and the monitoring point to the peak acceleration of the input seismic wave (9.58 m/s2) is defined as the PGA amplification factor [26, 27]. The PGA magnification factor can reflect the size of rock dynamic response (see Figure 4).
In order to monitor the movement and destruction laws of rock slopes and the development laws of structural planes under the action of earthquakes, three monitoring methods are drawn up: monitoring points, monitoring blocks, and measuring balls. The monitoring points can accurately reflect the influence of the dynamic response inside the foundation on the upper rock mass. The monitoring block can precisely monitor the motion trajectory and dynamic response of each dangerous rock. And the measuring ball specifically completes the stress change monitoring of the main control structure surface. A total of 10 monitoring points are arranged in the mudstone particle model in a horizontal and vertical arrangement. The interval between the horizontal and vertical adjacent points is 1.5 m, and they are all distributed at the zero points of the space x-axis (where the space coordinate of No. 1 point is (0, 9, 7)), which mainly realizes real-time monitoring of the dynamic acceleration of mudstone foundation. There are eight monitoring blocks arranged, corresponding to the eight rock blocks cut by the main control structural plane and the interlayer structural plane in the model. In the PFC, the rock blocks are grouped to realize real-time monitoring of the blocks by the group. The measuring balls are mainly arranged at the central axis of the two vertical main control structural planes, and the stress changes when the structural planes fracture in the area is mainly monitored. The main control structure surface of each rock block corresponds to two measuring balls, a total of sixteen are arranged, and each measuring ball is tangent to each other in the vertical direction. In order to evenly arrange the measuring balls of each layer on the structural surface, the diameter of the measuring balls shall be half of the corresponding rock formation height. The overall model monitoring program is shown in Figure 5.
Parallel bonding contact model and smooth joint model are generally used to evaluate rock mesoscopic materials in particle flow simulation. Given the results in the literature [26, 27], significant microscopic parameters of the final sandstone mudstone were obtained through laboratory mechanical tests and sensitivity analysis, as shown in Tables 2 and 3. The microscopic parameter calibration process is not described in detail in this paper.
In Table 2, d is the particle size, P is the particle density, Ec is the particle contact modulus, kc∗ is the stiffness ratio of the contact model, μ is the friction coefficient, λ is the radius multiplier, k∗ is the stiffness ratio of the particles, Eb is the parallel Bond modulus, σc is the tensile strength, c is the cohesion, βn is the normal critical damping ratio, and βs is the shear critical damping ratio. In Table 3, sj_kn is the joint normal stiffness, sj_ks is the joint shear stiffness, sj_μ is the joint friction coefficient, sj_c is the joint cohesion, and sj_φ is the joint friction angle.
2.2. Establishment of Dynamic Stability Calculation Formula
This paper analyzes the dynamic stability of the most severely weathered Model 4, using longitudinal seismic waves [28, 29].
Supposedly, a horizontal layered dangerous rock contains n macroscopic chains and m layers of rock. First, the rock block is chosen in random for analysis and defined as rock block m-n, where m is the rock layer number and n is the macrochain number [7]. In order to improve the force of the rock block, the mechanical model of the dangerous layered rock is divided into the case of longitudinal wave action. The force diagram of the rock block is shown in Figure 6, where Gmn is the weight of the rock block, Pymn is the vertical seismic force of the longitudinal wave, Pwmn is the fracture water pressure, the interaction between the upper and lower rock layers is q′(m + 1)n and q′mn, the adjacent rock block on the right side of the m-n rock block acts as bending moment Mm(n−1), tensile force Nm(n−1), and shear force Tm(n−1), Hmn represents the vertical height of the rock block, emn represents the depth of the fracture, e1mn represents the height of the fracture water, and the thickness of the rock block is Lmn. When m > 1, the m-n rock block is subjected to shear wave seismic force. According to the calculation formula of material mechanics and seismic inertial force, the force at the main control structure surface of the m-n rock block can be obtained as follows:where is the bending moment at the structural plane of the m-n rock block under the action of longitudinal waves, is the shear force at the structural plane of the m-n rock block under the action of longitudinal waves, is the tensile force at the structural plane of the m-n rock block under the action of longitudinal waves, ay(t) is the seismic acceleration, and and represent the friction between layers.
Based on the above results, combined with the literature [29], the fracture intensity factor under the action of earthquake time history can be obtained as follows:where KI1mn is the fracture intensity factor generated by the fracture water pressure, KI2mn is the fracture intensity factor generated by the bending moment under the earthquake, KI3mn is the fracture intensity factor generated by the tensile stress under the earthquake, is the fracture water pressure, is the bulk density of fracture water, F(emn) is the crack shape parameter, σmnmax is the maximum tensile stress, and σmn is the tensile stress. Combining formulas (1)–(7), the type I stress intensity factor under the action of longitudinal waves can be obtained as
In the same way, the type II stress intensity factor under longitudinal seismic waves can be obtained as follows:
From the joint stress intensity factor and fracture angle at the main control structural plane of the rock block m-n [30], the joint stress intensity factor and fracture angle under seismic longitudinal waves can be obtained as follows:
Finally, from the ratio of the fracture toughness of the main control structure surface KIC(t) to the formula (12), the time-history formula of the stability coefficient of the m-n rock block can be obtained as follows:
The interlayer load under the earthquake needs to be determined for calculating the stability of the horizontal-layered dangerous rock. This article assumes that adjacent rock formations are in contact with each other, and considering that each layer of rock interacts with each other, the m-th layer of the rock mass is now taken for analysis. The geometric model of rock contact is shown in Figure 6. Among them, the layer height of the m-th rock mass is Hm, and the layer thickness is Lm. The adjacent rock layers are the m−1 and m + 1 layers. The layer heights are Hm−1 and Hm+1, respectively. The value is the same as that of the m-th layer, and the deepest potential failure contact points of each rock layer are Am and Am−1, respectively.
In Figure 7, the deflection of m layer rock under its own weight and m + 1-layer pressure under the action of longitudinal waves iswhere γ is the bulk density of rock, E is the modulus of elasticity, and Im is the moment of inertia.
When supported by the lower rock formation, its deflection y”m2 is
The friction between the m layer and the m + 1 layer is
Assuming that the layered dangerous rock is in contact with i layers at the same time, the rock layers are numbered from bottom to top. Suppose the contact points of adjacent rock formations are, respectively, A2, ..., Ai, based on the principle of equal deflection of the contact sections of each adjacent rock formation, the following equation is given:
The interlayer load can be obtained by formulas (15)∼(19).
2.3. HHT Method
When analyzing the seismic dynamic response of dangerous rocks, scholars consider the rock mass vibration acceleration in the time domain that is more critical than the frequency domain. From the perspective of structural dynamics, structures have their natural frequency. Vibration frequency and vibration acceleration are equally crucial for the dynamic response of rock slopes. Therefore, besides the acceleration value of the vibration signal, the analysis of the dynamic response of the slope should also be concerned with the frequency change law of the vibration. Accordingly, this article uses the HHT signal processing method to further process and analyze the PFC acceleration data.
HHT (Hilbert–Huang transform) is a nonlinear and nonstationary signal analysis method newly proposed by Huang E in 1998. It is a breakthrough in the linear and steady-state spectral analysis based on Fourier transform in the past 100 years. At present, it is widely used in seismic wave spectrum energy analysis. The method consists of empirical mode decomposition (EMD) and Hilbert transform (HT) [31]. The Hilbert energy spectrum provides the energy calculation formula for each frequency, which expresses the energy accumulated by each frequency over the entire length of time. The specific formula can be found in the literature [32]. Compared with Fourier transform and other signal wave processing methods, HHT can adaptively generate “base,” which has the characteristics of complete adaptability. Besides, the instantaneous frequency obtained by HHT has locality. It does not require the entire wave to define the local frequency. Nevertheless, it can accurately make a three-dimensional time-frequency-amplitude map, which is difficult to achieve by other signal analysis methods such as wavelet transform.
3. Results
3.1. Analysis of Mechanics Formula Examples
A particular horizontal layered dangerous rock on Duwen Highway is taken as an example, as shown in Figure 8(a). The dangerous upper rock of natural lithology is sandstone, and the lower bedrock is mudstone. The natural gravity of sandstone is 24.7 kN/m3. The fracture toughness is calculated as a parameter in Table 4. In this example, there are four layers in total, including two main control structural planes. According to the chain law of dangerous rock obtained in literature [28], the first chain is called one chain and the second chain is called two chains, as shown in Figure 8(b). For the seismic wave, the first 25 s of the real UD wave monitored in Wolong in the 2008 Wenchuan earthquake in Sichuan was selected as the longitudinal wave calculation mode, and the acceleration time-history curve is shown in Figure 3(b).
(a)
(b)
The deflection equations of the first and second layers at point A2 are as follows [7]:
The simplified formula (20) can be obtained:
From the above-derived dangerous rock formula calculation, the stability coefficient of the dangerous rock block under the interlayer load can be obtained. To visually exhibit the attenuation law and caving sequence of the stability coefficient of the rock block brought by the earthquake, this paper deals with the stability coefficient of the rock block and draws the attenuation curve of the stability coefficient under the earthquake action (Figure 9). Since the stability coefficient under earthquake fluctuates continuously, this paper defines the peak stability coefficient in Figure 9 as the peak of each lower envelope, which is the maximum value of the stability coefficient when it fluctuates downward. At the same time, according to the intersection of each curve and the critical line, the caving time of each rock block can be obtained (Table 5).
The results of the calculation represent that the chain rock block 1-1 in the bottom layer is lower than the critical line value of 1.0 at 2.10 s, and the main control structure near the rock block 1-1 is destroyed first. When the earthquake action time reaches 2.92 s, the stability coefficient of the outer chain rock block 1-2 on the bottom layer reaches the critical line, and it is successively destroyed. Reason analysis for the priority failure of the two dangerous rocks at the bottom is that the interlayer load of the dangerous rocks at the bottom is the largest, the bottom of the rocks is unsupported, and the internal rocks are subjected to more significant bending moments of the external rocks. Subsequently, the second layer of rock blocks 2-1 and 2-2 is collapsed and destroyed in 4.75 s and 4.84 s, respectively. It can be seen that the first and second layers of rock blocks are all caving preferentially by the external chain. Due to the successive destruction of the internal rock masses and the continuous action of the earthquake, the 4-2 stability coefficient of the top internal rock mass reached the critical line at 6.07 s, and the 3-2 stability coefficient of the rock mass reached the critical line at 7.98 s. It can be judged that the failure mode of the top two layers of rock blocks is that the internal structural plane penetrates from top to bottom, and the rock blocks 4-1 and 3-1 outside the top are destroyed with the failure of 4-2 and 3-2. The failure times are 9.50 s and 10.90 s, respectively. Among them, the stability coefficient of the rock block decays fastest between 5 and 6 s, and the decay value is between 0.12 and 0.25; then, it tends to decay slowly. Through the above phenomena, we can get the law that horizontal layered dangerous rock blocks successively collapse under the seismic force. Furthermore, before the seismic acceleration has reached the maximum value, the dangerous rock blocks on the two main control structures have collapsed and destroyed. The calculation method for the stability coefficient of the dangerous rock under the influence of interlayer load and seismic dynamics established this time can reflect the failure law and caving time of flat layered dangerous rock.
3.2. PFC Dynamic Response Analysis
3.2.1. Failure Mode and Formula Verification
Under the unidirectional action of seismic waves, the failure characteristics of horizontal layered dangerous rocks are affected by different cavern depths and the crack penetration rates of the main control structure. The primary failure modes are divided into tensile cracking-horizontal sliding failure (Model 1), tensile cracking-caving-toppling failure (Model 2), caving-toppling-rotation failure (Model 3), and avalanching-slip-slip-rotation failure (Model 4) [33]. Figure 10 shows the failure shape and displacement cloud diagram of the four models after 25 s. It can be seen that Models 2 and 3 have the most massive displacement of the rock block on the top of the slope, up to 6.32 m, and their rotation failure tendency is particularly significant. Model 4 takes the second place, and Model 1 has the smallest particle displacement and failure tendency. Therefore, the rock block is more prone to collapse at a faster speed with a greater depth of the rock cavity and the fissure. The dynamic stability and failure mode are jointly controlled by the depth of the rock cavity and the fracture. When the rock block collapses, the adjacent rock blocks between the horizontal and interlayer will squeeze and collide with each other, which promotes the fracture of the smooth joint model to intensify. The pull-cracking failure of a single rock block leads to the overall dumping (slipping) rotation failure.
(a)
(b)
(c)
(d)
This paper compares the rock avalanche time calculated by Model 4 with the calculated value of the formula shown in Figure 11. It should be noted that, in the simulation calculation, the 4-1 rock mass was not destroyed alone but collapsed together with the 4-2 rock mass. According to reasonable assumptions, this article makes the avalanche time of about 9 s. The rock collapse of Model 4 is a destruction process that starts from the bottom and the top and gradually develops toward the middle. The calculation time of instability of each rock block considering the interlayer load and longitudinal wave is the same as the caving time calculated by the PFC simulation, which proves the rationality of the theory and simulation method and provides a basis for the dynamic response analysis.
3.2.2. PGA Amplification Factor
Under the action of an earthquake, the distribution of PGA amplification coefficients of horizontal layered rock blocks under different working conditions is shown in Figure 12. Both the depths of the rock cavity and the cracks have an important influence on the dynamic failure mechanism of horizontal layered dangerous rocks. Moreover, each model shows different PGA dynamic response characteristics. The PGA amplification coefficient value of the rock block increases gradually from Model 1 to Model 4. It can be seen that the deeper the rock cavity and fissure, the more prominent the PGA response to each rock block.
As for the Model 1 with the smallest cave and crack depth, the tendency of rock mass movement is smaller due to the superior integrity of the slope. The PGA magnification coefficient is distributed in a zigzag pattern with the increase in the rock layer. Furthermore, the peak is basically distributed in the bottom rock block and the internal macrochain No. 2. It shows that the dynamic acceleration response of Model 1 dangerous rock is relatively uniform under the action of an earthquake, and the bottom and internal chain rock block 2 are more prone to damage. In Model 2, the acceleration response under the action of an earthquake is the most significant around the No. 1-2 rock mass and the top 4-2 rock mass, the possibility of priority failure is the largest, and there is a noticeable elevation amplification effect. As the rock layers increase, the PGA amplification coefficients of the rock blocks of Models 3 and 4 show a “U”-shaped distribution with the change in rock formations. The reason is that, in Models 3 to 4, weathering is intensified, the bottom rock block has a preferential caving trend, and the dynamic response is more prominent. As the top rock mass develops from top to bottom cracks, the dynamic response is enlarged. Therefore, after the weathering of rock cavities and fissures deepens to Models 3∼4, the acceleration response of the bottom rock block and the top rock block without base support under earthquake action is more prominent and more destructive.
In general, the dynamic response degree of rock PGA is related to the degree of weathering. The acceleration response of the first rock block and the top rock block of the Model 3 and Model 4 under earthquake action is more prominent. It can be seen that the greater is the degree of weathering, the greater is the dynamic response of the rock block PGA, which is affected by the change in the rock layer height.
This paper analyzes the PGA magnification factor of the horizontal and vertical monitoring points of the mudstone on the foundation, as shown in Figure 13. The horizontal and vertical distances of the monitoring point in Figure 13 refer to the distance between the monitoring point and the coordinate origin of the lower left of the model. It can be seen from Figure 13(a) that, in the vertical z-direction, the PGA amplification coefficient of each model increases linearly with the elevation change, and there is an elevation amplification effect. Furthermore, as the rock cavity and fissure deepen, the dynamic response of mudstone along the vertical direction gradually intensifies. It shows that the mudstone particles near the upper sandstone are greatly affected by the dynamic damage, and the dynamic response is significant that it gradually decreases with the elevation downward. In Figure 13(b), it can be seen that, as the caverns and fissures deepen, the law of the gradual intensification of mudstone dynamic response is consistent with the vertical direction, and with the development of the horizontal direction, there are large fluctuations. Through comparative analysis, the PGA amplification factor of mudstone near the slope surface is the largest, whereas the dynamic acceleration response is the most sensitive. Here, the PGA magnification coefficient curves of each model showed a “V”-shaped distribution within 0∼4 m from the horizontal distance of the slope. This is because the range corresponds to the two main control structural planes of the upper dangerous rock slope. Due to the seismic action, the dynamic response of the two horizontal ends of the structural plane is relatively large. It can be seen that, during the propagation process of seismic waves, it acts on the mudstone slope at the bottom first and then propagates the power to dangerous upper rock to destroy it. Subsequently, the upper rock mass reacts on the mudstone slope at the bottom. The PGA response of the mudstone slope interacts with the upper sandstone rock mass and is affected by the degree of weathering, and its PGA response has a particular elevation amplification effect.
(a)
(b)
3.2.3. Stress Analysis
In order to clarify the stress distribution of the main control structural plane of horizontal layered dangerous rock in the process of instability and failure, the stress changes in the two main control structural planes in the x, y, and z directions were monitored by measuring balls. Since PFC can only monitor the time-stress curve, the peak stress in this article is defined as the maximum stress value monitored by the model measuring ball within 25 s of the earthquake. In order to simplify the analysis, the peak stress method is adopted, as shown in Figure 14. The height of the measurement sphere center refers to the vertical distance from the measurement sphere center to the coordinate origin in the lower-left corner of the model. It can be seen from Figures 14(a) and 14(b) that, in the x-direction, the peak stress of each model is more massive in the bottom layer. As the height increases, the peak stress gradually decreases. The stress on the inner No. 2 chain is more massive than that on the outer No. 1 chain. This is mainly due to the large shear failure of the internal main control structure, which intensifies the dynamic response. Prior failure of the dangerous rock mass produces more significant stress, the contact force chain between the rock blocks is weakened, and the stress of the upper rock mass is successively weakened. Due to the continuous development and deepening of cracks (Models 1–4), the peak stress in the x-direction decreases, so Model 4 has the weakest stress effect under earthquake action. Therefore, in addition to the z and y directions, large shear stresses also occur in the longitudinal direction x of the slope, up to 18900 kPa, distributed in Model 3. There is a large stress fluctuation with the height of the slope, and both attenuate to the minimum at the top.
(a)
(b)
(c)
(d)
(e)
(f)
In the horizontal y direction, due to the displacement of the rock strata caused by the seismic dynamics and the collision and extrusion of the rocks, the peak stress changes intensify. The rock stress of Model 1 and Model 2 on the outer No. 1 chain first decreases with height and then increases and finally decreases after reaching the top. The maximum stress of the external chain of Model 1 can reach 45800 kPa. The above results are because when the external chain rock mass is broken, it will produce a considerable rock dislocation, and the joint contact model will undergo shear deformation. The shear stress in the y direction is too large, especially in the 1-1 rock block (elevation 11 m) and the 3-1 rock block (elevation 21.5 m), and the peak stress is the most obvious. On the inner No. 2 chain, the peak stress value range of each model is about 30,000 kPa larger than that of the outer chain. Among them, the stress changes in Model 2 and Model 3 are the most significant. In Model 2, the stress on the No. 2 chain from the first layer to the second layer increased to 61,800 kPa and then quickly decreased to 12,000 kPa. The reason for the analysis is that because Model 2 is a caving-toppling failure, a large amount of compressive stress will be generated in the bottom rock mass during the dumping process. After the earthquake, the third and fourth layers of rock have collapsed; therefore, the stress on the structural surface of the 2-2 rock mass of the second layer is significantly reduced. In the failure process of Model 3, the stress on the primary control structure of the bottom 1-2 rock blocks decreased from 73,600 kPa to 2,500 kPa, which showed an enormous change. The reason for the analysis is that the internal chain failure mode of Model 3 is a dumping-rotation failure, and the upper rock mass is an avalanche. Most of the weight is concentrated in the measuring ball area at the lower-left corner of the 1-2 rock block. Shearing between particles continuously occurs. The squeezing action produces a large amount of horizontal stress.
In the vertical z-direction, the peak stress is mainly caused by the vertical shear failure of the force chain at the main control structure surface. It can be seen from Figures 14(e) and 14(f) that, due to the smaller fracture depth and more structural surface force chains in Model 1, the shear stress generated during earthquake damage is more extensive. Moreover, there is a more obvious elevation amplification effect on the outer No. 1 chain; the maximum vertical shear stress can reach 52,000 kPa, which is distributed on the third layer No. 3-1 rock block (elevation 21.5 m). Subsequently, in the fourth layer, the No. 4–1 rock mass (elevation 24.7 m) was reduced to 5000 kPa. The reason is that the top rock block contacts the model with less fracture. For other Models 2 and 3, the vertical peak stress gradually decreases due to the increase in the crack depth of the structural plane.
Furthermore, in the outer No. 1 chain, the peak stress fluctuation gradually stabilizes as the crack deepens. The variation range of the vertical peak stress on the structural surface of chain 2 is 37000 kPa less than chain 1. Among them, the vertical peak stress of Model 1 and Model 2 showed a “U” distribution characteristic with height change, Model 3 gradually decreased from 10500 kPa to 1100 kPa with height change, and Model 4 showed a “sawtooth” distribution.
In summary, the development of cracks on the main control structure has a more significant impact on the damage of rock slopes under earthquake. As the cracks on the main control structure deepen from 25% to 62.5%, the number of force chains decreases. Therefore, the intensity of the peak stress decreases gradually, and the stress curve becomes flat with the elevation subsequently, but the degree of structural plane instability increases accordingly. In the horizontal x and y directions, the stress changes reflected by the internal No. 2 chain are more severe, the contact force chain is more likely to break, and the internal structural surface is destructive. In the vertical z-direction, the stress response of the outer No. 1 chain is more robust.
3.2.4. HHT 3D Time-Frequency Analysis
The HHT method is composed of empirical mode decomposition (EMD) and Hilbert transform (HT). The basic steps are as follows:(1)First, use EMD to decompose complex signals into a finite number of intrinsic mode functions (IMF). Use cubic spline function to interpolate the extreme points on the original signal x(t), connect the envelopes fitted by the upper and lower extreme points, and define the mean value of the two envelopes as m(t). The result of subtracting from the mean is as follows [31, 32]: In the formula, x(t) is the original signal, m(t) is the mean value of the upper and lower envelopes, and y(t) is the difference. Using y(t) as the original signal, repeat the above steps until the stopping criterion is satisfied after n iterations. The criterion is as follows: Among them, the SD value is generally between 0.2 and 0.3 to stop the iteration. At this time, c1 = y1(t) obtained after n screening is the first IMF component, and then, the difference between x(t) and c1 (the value r(t) = x(t)−c1(t)) is the new signal, and the IMF components c2, c3, ..., cn can be obtained in sequence. So far, the original signal can be decomposed into the sum of n IMF components and the margin rn: where c(t) is the modal component and r(t) is the residual margin.(2)Perform Hilbert transformation on each intrinsic mode function component to obtain the instantaneous frequency and amplitude of each IMF component over time. First, perform the Hilbert transform on the component c(t) as follows:
In the formula: H is the Hilbert spectrum and P is the Cauchy principal value.
Construct the analytical signal z(t):
Then, the amplitude function can be obtained:
Phase function:And then, get the instantaneous frequency:
After applying Hilbert transform to each IMF component, where Re means take the real part. The Hilbert spectrum can be obtained by expressing the above equation as a function of time domain and frequency domain. Then, the Hilbert energy spectrum can be obtained by integrating the square of the amplitude against time:where E is the Hilbert energy spectrum.
When processing the dynamic response signal of the rock slope under earthquake, this paper performs HHT transformation on the PFC acceleration time-history curve of the rock block. This paper selects the internal macrochain rock blocks No. 2 (1-2 rock block, 2-2 rock block, 3-2 rock block, and 4–2 rock block) in each model as the vertical analysis object, and the bottom layer of each model rock blocks (1-1 rock block and 1-2 rock block) is used as horizontal analysis objects. Use the Hilbert-Huang signal processing method to obtain the HHT time-frequency diagram with time domain, frequency domain, and amplitude (take the bottom and top rock blocks), as shown in Figure 15. For easy analysis, the instantaneous Fourier dominant frequency corresponding to the obtained maximum amplitude is extracted, and the obtained curve is shown in Figure 16.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
It can be seen from Figure 16 that the instantaneous dominant frequencies of each model rock block under the action of an earthquake are all distributed between 0 and 10 Hz. Moreover, the dominant frequency fluctuations of each rock block in the same model are different. Since the stability of the outer chain 1-1 rock block decreases with the increase in the weathering depth of the rock cavity and the penetration rate of the fracture, its instantaneous Fourier frequency also decreases. In Model 1, the Fourier dominant frequency of each rock block fluctuates between 2 and 2.75 Hz. The small variation range indicates that this type of slope has better integrity in the frequency domain. In Model 2, the dominant frequency of rock blocks varies significantly with elevation, and the high frequencies are concentrated in the second and third layers of the rock mass. In Model 3 and Model 4 with a higher degree of weathering, the main frequency increases from the external macrochain to the internal macrochain, and as the elevation increases, the main frequency gradually increases.
It can be seen that the frequencies corresponding to the peaks of the four dangerous rock models under the earthquake are concentrated within the low frequency (0∼10 Hz). The distribution of most amplitudes in the frequency domain is relatively monotonous and fluctuates with time on their respective Fourier main frequencies. From the analysis of the PGA magnification factor, it can be seen that the deepening of weathering leads to an increase in the dynamic response of the slope so that each rock block amplitude in the three-dimensional time-frequency diagram also increases with the rock cavity and fissure deepening.
In Model 1, the prominent frequency bands of each rock block are distributed between 0 and 90 Hz. The shape of the signal wave is mainly a combination of a single main peak and multiple small peaks. In Model 2, the main frequency band of the rock block has been widened. At the same time, the main peak of the signal wave increases. In Model 3, the main frequency band of the rock block increases further. Except for rock blocks 1-2, the other rock blocks all fluctuate at a frequency of 250 Hz, with a fluctuation range of 0–8 m/s2. This is a local high-frequency effect. At the same time, the number of rock blocks of the double main peak signal wave increase. In Model 4, the main frequency band of the rock block is distributed between 0 and 160 Hz, with the widest frequency band. Among them, the frequency band from block 1-2 to block 4–2 increases from 0–30 Hz to 0–160 Hz due to the elevation amplification effect of the frequency bandwidth. At the same time, high-frequency-amplitude fluctuations occur in each rock mass, and the fluctuation interval increases to 0–12 m/s2. At this time, there are multiple main peak signal distributions in dangerous rocks 1-2.
It can be seen that, in the case of deepening weathering, the possibility of interaction increases due to the intensified dynamic response of the slope. At this time, the overall frequency band of the slope will be widened, and the elevation amplification effect is present in Model 4. At the same time, high-frequency effects will increase accordingly. Moreover, the signal wave will evolve with the weathering of the slope, and there will be multiple main peaks at the dominant frequency. In earthquake disaster protection projects, as the degree of dangerous rocks weathering intensifies, in addition to strengthening the stability of single dangerous rocks, more treatments should be made on structural surface cracks between rock blocks (including interlayers). This can avoid the amplitude energy fluctuation caused by the mutual vibration of the rock masses. At the same time, Model 1 (the fracture penetration rate of the structural plane is 25%, and the thickness of the rock cavity is 4.5 m) should be considered under the earthquake action of the resonance failure effect of 0∼90 Hz. Model 2 (the penetration rate of structural plane cracks is 37.5%, and the thickness of the rock cavity is 3 m) should consider the effect of resonance failure of 0–100 Hz. Model 3 (the fracture penetration rate of the structural plane is 50%, and the thickness of the rock cavity is 1.5 m) should consider the resonance effect of 0∼120 Hz. In Model 4 (the fracture penetration rate of the structural plane is 62.5%, and the thickness of the rock cavity is 0 m), the resonance effect of 0∼160 Hz should be considered. Among them, the slopes under the four working conditions all focus on the seismic low frequency (0∼10 Hz). Models 2 and 3 also need to pay attention to the impact of high-frequency 250 Hz on the damage of dangerous rocks.
4. Conclusions
(1)The calculation methods of the horizontal layered dangerous rock slope under earthquake action and of the stability coefficient are established based on the dynamics theory, taking the seismic timeliness into consideration. The results show good consistency with the PFC simulation results, verifying the rationality of the formula and the PFC method. This method can not only reveal the attenuation law of the horizontal layered dangerous rock stability under earthquake (including timeliness) but also reveal the time and sequence of each rock block caving under earthquake action.(2)Under the action of an earthquake, the PGA distribution of Model 1 is relatively uniform. Model 2 is distributed in a zigzag pattern, and both Models 3 and 4 have a “U”-shaped distribution. It shows that the greater is the degree of sandstone joints weathering, the more severe is the dynamic response of the rock mass acceleration, and the degree of dynamic acceleration response is affected by the crack propagation process in the upper sandstone and exhibits a specific elevation amplification effect.(3)The peak stress of the main control structure surface of the slope model under the action of the earthquake gradually decreases with the weathering and the vertical elevation of the rock block. In the horizontal x and y directions, the internal stress of chain two changes more drastically, and the contact force chain is more likely to break. The stress response of the outer No. 1 chain in the vertical z-direction is severe. Among them, the vertical peak stress of Models 1 and 2 shows a “U”-shaped distribution with the increase in rock block elevation. The vertical peak stress of Model 4 is distributed in a “sawtooth” shape.(4)Based on the HHT signal processing method, the analysis found that the main frequency of each slope model under the action of the earthquake is mainly concentrated in the low frequency. As the degree of slope weathering increases, the overall frequency band of the slope widens, and the high-frequency effect increases. The corresponding time of the main peak is the same as the rock block collapse time of the PFC slope model, and there are multiple main peaks at the main frequency. It shows that the dynamic response of rock mass acceleration under earthquake is positively correlated with the degree of slope weathering. It is recommended that earthquake disaster protection projects should pay attention to the impact of low-frequency (0–10 Hz) and high-frequency (250 Hz) earthquakes on slope stability.
Data Availability
The table and figure data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was financially supported by the National Key Research and Development Project (2016YFC0802203), National Natural Science Foundation (51678097), Science and Technology Research Project of Chongqing Municipal Education Commission (KJQN201800706), General Project of Chongqing Natural Science Foundation (cstc2020jcyj-msxmX0218), and Science and Technology Research Project of Chongqing Education Commission. The authors would like to express their gratitude to EditSprings (https://www.editsprings.com/) for the expert linguistic services provided.