Numerical research method of rock slope with intermittent fractures based on parallel bond contact model

Within the domain of rock slopes characterized by intermittent fracturing, the fracture distribution and the mechanical properties of the interstitial rock bridges significantly influence the stability of the slope. The parallel bond contact model is advantageous in replicating the mechanical behavior of rock particles. This research introduces a numerical methodology for analyzing rock slopes with intermittent fractures using the parallel bond contact model. Initially, the model’s microscale parameters are refined through calibration with empirical data derived from macroscopic mechanical tests on rocks. Following this, the discrete element modeling software is employed to construct a detailed rock slope model. This model incorporates a smooth joint approach to define the intermittent fractures, enabling the creation of slope models with varying configurations of coplanar rock bridges and diverse rock types. The research methodologically investigates the mechanical properties and failure patterns of rock slopes under a spectrum of variable combinations. The findings reveal that slopes with multiple rock bridges demonstrate progressive failure and interlocking phenomena during their load-deformation cycles. These insights provide a foundational understanding for the analysis of catastrophic mechanisms and the stability assessment of rock slopes.


Introduction
Throughout extensive geological structural processes and environmental influences, rock masses evolve with a multitude of structural planes, each varying in scale, orientation, and mechanical characteristics.These planes typically exhibit discontinuity, segmented and interspersed by unfractured rock segments, known as rock bridges or lock segments, which are crucial in determining the stability of rock slopes.Numerous instances of landslides have illustrated that the failure of rock slopes with intermittent cracks is a progressive phenomenon, evolving from localized to comprehensive areas.This progression is significantly influenced by the interplay of various factors, including the integration of rock bridges, diversity in rock types, and the methodologies of loading [1,2].
When subjected to adverse external loadings, rock slopes harboring intermittent cracks experience a cascade of shear failures in their lock segments, diminishing to residual strength.This degradation leads to the downslope force surpassing the resisting force, culminating in landslides.At the microscale, this process entails rock particle flow and micro-fracturing, contributing to internal damage and strength reduction.Consequently, substantial research has been dedicated to understanding the initiation and coalescence mechanisms of rock bridges under diverse stress conditions [3][4][5].Wong et al. explored the impact of rock bridge length and inclination on crack propagation in artificially fractured rock-like material specimens [6].Park and Bobet developed gypsum specimens with three differently angled and connected cracks, conducting experimental investigations into the onset, progression, and merging of cracks under uniaxial compression [7].Building upon this methodology, Cao et al. discovered that the interplay of multiple cracks diminishes the relative sliding of coplanar cracks [8].Zhao et al. and Xue et al. performed uniaxial compression tests on rock-like specimens with various configurations of two parallel cracks, integrating digital image correlation technology (DIC).They observed that the bridge angle significantly influences the crack propagation pathway and the strength of the rock [9,10].These empirical studies have established a solid rock mechanics foundation for elucidating the disaster evolution mechanism in rock slopes.
In the realm of model experimental studies, Huang and colleagues carried out physical model experiments on rock landslides, scrutinizing the initiation, progression, and failure mechanisms of locktype landslides [11,12].Gong et al. employed analogous materials in their physical model experiments to probe the compound toppling failure mechanism and its evolution in layered rock slopes under varied stress conditions [13].Dong et al. explored the effects of seismic activities on lock-type rock slope disasters, conducting comparative physical model tests on a large vibration table and formulating instability criteria for lock-type rock slopes [14].Zhou and associates utilized extended finite element methods to simulate rock slopes with intermittent joints, numerically examining the brittle failure mechanisms of jointed rock slopes and presenting an effective analysis of joint crack mechanical behavior [15].While physical model experiments adeptly analyze the catastrophic process in rock slopes with intermittent cracks, they entail high experimental costs and extended durations, with challenges in thoroughly examining diverse conditions.Numerical simulation methods have gained widespread application in investigating geotechnical engineering disasters and instability mechanisms.
Accordingly, this paper introduces a rock slope model with intermittent fractures, developed using the parallel bond contact model in particle discrete element software.It analyzes the impact of varying numbers of rock bridges and rock property parameters on the slope model's mechanical characteristics and failure modes.

Parallel bond contact model
2.1.1.Conception.In the parallel bond contact model, the interaction between particles is replicated using normal and tangential springs, effectively representing the inter-granular cementation characteristic of rock minerals.When the cementation integrity is compromised, the model defaults to a conventional linear contact model.As a result, the parallel bond contact model plays a vital role in accurately simulating the mechanical property variations of rock materials in a collective assembly of numerous particles.It operates by combining parallel-bonded units with linear contacts in both normal and tangential directions.This model comprises two fundamental components: a linear-elastic, frictional unit that transmits force but cannot convey bending moments and is susceptible to slipping under excessive shear force, and a second unit designed to withstand both force and bending moments, which fails when its strength threshold is surpassed. .The total contact force in the parallel bond is given as follows: where n i F represents the normal component along the local coordinate system at the contact interface, and s i F signifies the tangential component.In computations, vector forces can be transformed into the product of corresponding scalars and unit vectors.
Before iterative calculations, contact forces are typically reset to zero (i.e., setting lin_mode=1), initializing the total contact force ( i F ) and total bending moment ( 3 M ) to null.uubsequently, incremental calculations are conducted at each timestep, cumulatively adding the changes in force and bending moment to the previous timestep's values.The increment in each timestep can be expressed as: The increments per timestep for normal and tangential displacements are represented as  represents the relative rotational angle, and I represents the particle's section polar moment of inertia.
Consequently, the maximum tensile and shear stresses acting on the bonded unit can be defined as: where A represents the cross-sectional area at the bond, and  is defined as the moment contribution parameter, ranging from 0 to 1.

Calibration.
In rock mechanics, rock attributes are conventionally characterized using macroscopic mechanical properties such as strength and modulus.However, for the Discrete Element Method (DEM), the parameters in the modeling process are on a microscopic scale, encompassing aspects like bond radius and inter-particle friction angle.Consequently, it becomes crucial to forge a link between the microscopic parameters used in PFC (Particle Flow Code) and the macroscopic mechanical properties of rock masses, a procedure known as parameter calibration.Within PFC, model parameters are bifurcated into modeling parameters and microscopic mechanical parameters.
During the calibration process, ensuring that modeling parameters mirror the actual model's specifications is imperative, followed by the meticulous determination of microscopic mechanical parameters via extensive numerical testing.Key modeling parameters include particle size distribution and sample preparation methods, with particle size distribution typically adopting a 1.66 size ratio and adhering to a normal distribution to maintain a balanced porosity ratio.This approach aids in realistically simulating rock mass characteristics while circumventing excessive suspended particles that might impede computational efficiency.
In this investigation, the sample preparation method employed is the servo sampling method, adept at simulating the formation processes of sedimentary rocks like limestone and shale.Initially, particle aggregates are generated through random distribution.Due to considerable initial overlaps and stress induced by corresponding commands in the PFC software, the volume of the particle aggregate undergoes modification using wall servo until the wall-transmitted pressure aligns with the preset servo value.This procedure facilitates the incorporation of pertinent contact parameters into the particles, resetting the pre-established inter-particle overlap to zero, and initiating incremental calculation mode, thereby mimicking the consolidation and cementation processes of rock particles in nature.
Leveraging the macro and micro parameters of rocks derived from laboratory experiments conducted by our research team, and drawing on experimental and numerical analyses of limestone and mudstone [17,18], three sets of rock micro parameters with distinct mechanical characteristics are calibrated.This calibration is undertaken to examine the influence of lithology on the failure process of rock slopes, as 5 depicted in Tables 1, 2, and 3.The first method, while relatively simple, necessitates a substantial particle count.In cases where the fracture size is significantly smaller than the mean particle size, the resulting pre-existing fracture exhibits reduced smoothness.The second approach yields smoother fracture surfaces but involves a more intricate modeling process.Additionally, this method faces the challenge of considerable unbalanced forces at the pre-existing fractures during servo pre-consolidation, leading to prolonged computation times relative to other techniques.The third method, more straightforward in nature, effectively aids in the formation of pre-existing joints with specific mechanical parameters.After evaluation, this study chose to employ the smooth joint model for generating pre-existing fractures.

2.2.2.
Rock slope modeling procedure.The precision of numerical simulation results enhances with the reduction in average particle size.However, a higher particle count markedly diminishes computational efficiency, necessitating an equilibrium between simulation precision and computational expediency.In simulating rock slopes, particle size in non-focal areas does not influence the outcome of the analysis.Thus, varying particle sizes may be employed in different model sections.The process for modeling rock slopes with intermittent fractures is delineated as follows (refer to Figure 2): (1) Construct a 2.7 m x 2.3 m rectangular boundary wall and within it, generate dividing walls to section the slope model into three zones.Particles of varying sizes (0.004 m to 0.0066 m, 0.0132 m to 0.02 m, and 0.025 m to 0.52 m) are then distributed in these zones.The smallest particles are utilized in the critical rock bridge area.uince negligible deformation occurs within the slope body during computation, the largest particles are used in this area, as they do not affect the results.
(2) The 'solve' command is employed to stabilize all particles.To enhance computational efficiency, the 'clean' command removes the dividing walls and rebalances the system.The boundary walls' normal movement is controlled using functions written in the FIuH language, ensuring that the force exerted by the particle aggregate on the walls reaches 1 MPa.
(3) Post pre-consolidation servo, all particles undergo gravity acceleration for initial self-weight balancing to establish the initial geostress.The top boundary wall is then removed, with normal constraints applied to the left, right, and bottom walls, simulating semi-infinite boundary conditions.
(4) Adhering to the slope's configuration, particles outside its range are eliminated using the DFN command, shaping a complete slope body.Concurrently, DFNs are created based on the centroid, trace length, and dip angle of pre-existing fractures, and the smooth joint model is applied to the intersecting particle contacts with the DFN, culminating in a slope model with intermittent fractures.
In this study, considering the impact of rock bridge quantity and lithology, slope models with 1, 3, and 5 coplanar rock bridges are established, utilizing the aforementioned three sets of rock micro parameters.The inclination of these coplanar rock bridges is uniformly set at 50°.

Mechanical properties and acoustic emission events of rock slopes
During rock failure, dislocation between internal grains under stress precipitates microcrack formation, releasing energy and generating elastic waves within the rock mass.Acoustic emission probes in laboratory tests monitor these waves, reflecting internal damage within the rock mass.These acoustic emissions correspond to bond break events in PFC's (Particle Flow Code) contact model.Figure 3 illustrates stress-strain curves and acoustic emission event variations in slope models with diverse variable combinations.
The nine groups of slope models demonstrate initial elastic deformation, pre-peak elasto-plastic yielding, and post-peak failure stages in their stress-strain profiles.In slope models with different contact parameters, the stress-strain curves show linear growth during the elastic deformation stage, though with varying slopes.The Rock III contact parameter slope model, indicative of a lower elastic modulus, has a gentler slope of linear elasticity compared to the Rock I parameter model.However, the Rock III parameter model exhibits more pronounced elasto-plastic yielding and post-peak failure stages, as depicted in Figure 3(c), (f), and (i).After reaching the initial peak, the stress-strain curve drops slightly before steadily rising and culminating in model failure after attaining the peak again.Figure 3(a), (d), and (g) suggest that the plastic yielding stage is inversely proportional to the rock's elastic modulus.Comparing models with varying rock bridges, more bridges result in a steeper post-peak curve drop, signaling more pronounced brittle failure in models with numerous bridges.Multiple bridges exacerbate stress and strain distribution unevenness within the rock mass during later loading stages, leading to instantaneous brittle failure at peak load.
In early loading stages, minor cracks appear within the rock mass, accompanied by sparse acoustic emissions.With increasing load, internal cracks expand and interconnect, culminating in rapid fracturing and a surge in acoustic emissions at peak load.These cracks include both shear and tensile types, with shear cracks predominant in early loading stages and tensile cracks in later failure stages.Notably, the peak in acoustic emissions aligns with the peak rock mass stress.However, some models exhibit a lag in acoustic emissions relative to the stress-strain curve, as seen in Figure 3(d) and (g).This indicates that even at maximum bearing capacity, some rock bridges retain the ability to support local loads, manifesting a progressive failure effect in the load-deformation process of models with multiple bridges.upecifically, while some bridges reach peak strength and fail completely, others remain in the elastic deformation or elasto-plastic yielding stages.Through deformation coordination and stress transfer, the remaining bridges begin to shoulder the majority of the load.

Influence of the number of bridges on the mechanical properties
In calibrating the micro parameters for the contact model, the mechanical properties of three distinct rock types are considered.Rock I demonstrates the highest strength, followed by Rock II, with Rock III exhibiting the lowest strength.Accordingly, slope models calibrated with these rocks' micro parameters However, the peak stress in models with identical parameter sets varies with the number of rock bridges.For models using Rock I's and II's parameters, an increase in rock bridges correlates with higher peak stress.In contrast, models with Rock III's parameters show a slight decrease in peak stress as the number of rock bridges increases.From the previous analysis, it is clear that the failure of slope models with multiple rock bridges is marked by progressive failure.This type of failure in rock masses is essentially a process where localized failure, due to uneven load distribution, evolves into total failure.In rock masses with multiple bridges, each bridge, under varying stress states and deformation stages, contributes differently to the mechanical strength of the mass, which is determined by the peak strength, residual strength, and the load borne by the elastic deformation of these bridges.
(a) Peak stress (b) Failure strain Figure 4. Peak stress and strain of slope models with different number of bridges

Failure modes of rock slope under different conditions
Figure 5 showcases the failure morphologies of slope models with various variable combinations.In the figure, the blue lines represent the failure of parallel bond contact models at specific locations, leading to particle cohesion disintegration.These findings align with the variations in acoustic emission events depicted in Figure 3, mirroring crack propagation within the slope model.As the number of rock bridges in the model increases, their interlocking effect intensifies, diminishing wing crack propagation at initial crack tips, as evident in Figure 5(a), (d), and (g).The pre-fabricated initial cracks modify the principal stress field of the intact rock mass, creating zones of compressive stress concentration between the inner endpoints of these cracks and confining the tensile stress field externally.This dynamic encourages shear failure in the rock bridges between initial cracks.Consequently, numerous microcracks surface in the rock bridge areas, as seen in Figure 5, indicating a significant reduction in rock mass strength where rock particles transition from a parallel bond contact model to a linear contact model.In models with multiple rock bridges, an increased prevalence of cracks near the upper loading end denotes more acute rock bridge failures, while fewer cracks at distant locations reflect the progressive failure behavior of the rock bridges.Furthermore, as the rock mass's elastic modulus diminishes, there is a notable increase in microcracks and a few wing cracks, suggesting that a lower elastic modulus diminishes the

Conclusion
The numerical research approach proposed for rock slopes with intermittent fractures utilizes the parallel bonded contact model, adept at simulating the mechanical behavior of rock particles effectively.The predefinition of fractures is executed using a smooth joint model, a process that is relatively straightforward and enables the creation of predefined joints with specific mechanical parameters.To optimize computational efficiency, the slope model is segmented, and particles of varying sizes are allocated.Employing this numerical approach, rock slope models are constructed, featuring different numbers of rock bridges and varying rock types.
This methodology facilitates the investigation of mechanical characteristics and failure modes of rock slopes under diverse variable combinations.The findings indicate that rock slopes with intermittent fractures display progressive failure effects and interlocking effects throughout the loading-deformation process.Additionally, the mechanical strength of the slope body in models with multiple rock bridges is determined by the combined influence of peak strength, residual strength, and elastic deformation of rock bridges at various locations.

Figure 1 .
Figure 1.Diagram of Parallel bond contact model [16] The model's primary control parameters include: normal stiffness ( n k ), tangential stiffness ( s k ), normal strength ( c  ), tangential strength ( c  ), and bond radius ( R ).The total contact force in the parallel bond is given as follows:

Figure 2 .
Figure 2. Rock slope with fractures modeling procedure

Figure 3 .
Figure 3. utress-strain curve and acoustic emission events variation of slope model characteristics, as depicted in Figure4.The peak stress of the slope model utilizing Rock I's micro contact parameters is around 40 MPa, while the models employing Rock II's and Rock III's parameters exhibit peak stresses approximately at 25 MPa and 18 MPa, respectively.

Table 1
Micro parameters of rock I

Table 2
Micro parameters of rock II

Table 3
Micro parameters of rock III Prefabricated crack modeling method.Presently, in Particle Flow Code (PFC), three primary methods are employed to create pre-existing fractures: (1) Formation of a complete particle aggregate using bonding parameters, followed by the deletion of particles at pre-determined fracture locations, based on center, dip angle, and trace length, after assembling a complete rock model.(2) Direct construction of an aggregate incorporating pre-existing fractures, specifically by establishing walls at the fracture sites, then executing particle placement, pre-consolidation, and bonding.(3) Utilization of the smooth joint model for generating pre-existing fractures.The built-in Discrete Fracture Network (DFN) commands in PFC facilitate the creation of DFNs anchored on the pre-existing fracture coordinates, applying the smooth joint model to the particle contacts intersecting within a defined proximity of the DFN.