Abstract
In order to study the influence of bedding structure on the mechanical properties of rock mass structure, this paper takes the layered carbonaceous slate in the Muzhailing tunnel as the research object, and studies the mechanical properties of layered slate under uniaxial compression and Brazilian splitting. In the numerical simulation, the discrete element software MUSEN (bonded particle model, BPM) is used to simulate the failure behavior of layered rock mass in the process of quasi-static compression. The main research contents of this paper are as follows: (1) the relationship between the microscopic Young’s modulus and the macroscopic Young’s modulus. (2) Influence of tangential shear strength and normal tensile strength on particle model. (3) The influence of different particle distributions on the simulation results. (4) The particle model with different bedding inclinations was simulated numerically. Through the comparative analysis of simulation and experiment, the results show that the mechanical response and failure mode have good consistency, and the splitting evolution process of layered slate with different bedding inclination angles is reproduced.
1. Introduction
Layered rock masses with structural planes are widely distributed in geotechnical engineering, among which there are usually well-developed bedding planes in rock masses such as shale, slate, and coal. These beddings of different scales cut the complete rock mass into transversely isotropic rock mass, and the resulting differences in rock mechanical properties have an important impact on the safety and stability of engineering rock mass [1, 2]. Because the layered rock mass has obvious transverse isotropy, the strength of the rock mass is not only related to its own strength, but the distribution and properties of the structural plane have a certain influence on the rock mass, and the form of damage and splitting is often related to the occurrence of the structural plane. Therefore, it is of great significance to carry out tests on the tensile and compressive properties of layered rocks and use efficient numerical software to improve the simulation efficiency for guiding large-scale engineering practice.
In order to study the properties and fracture modes of layered rock mass, the finite element method and the discrete element method (DEM) are gradually developed. In the finite element framework, numerical simulations can be used to predict the postpeak behavior of materials through elastic-plastic models [3–10]. The initiation, propagation, and penetration of rock defects under compressive loading can be simulated by near-field dynamics, conjugate bond near-dynamics, and meshless numerical methods. On the other hand, the numerical method can model and distinguish wing cracks, oblique secondary cracks, quasi-coplanar secondary cracks, and antiwing cracks [11–15]. However, there are not size effects and the difficulty of the parameter calibration in many simulation methods, such as nonlocal method [11–13], field-enriched finite element methods [15], and general particle dynamics [14]. On the other hand, the finite element numerical simulation is weak to study the multiscale and multicrack of materials. Therefore, the discrete element method is gradually favored by scholars. In this pioneering work, Cundall et al. treat rigid disks and spheres as discrete particles that interact to model the mechanical behavior of components [16, 17]. With the development of discrete element method and the gradual improvement of computer performance, computational efficiency is particularly important. The current discrete element software PFC has high computational efficiency in two-dimensional, but relatively low computational efficiency for large-scale three-dimensional models. Based on the discrete element software MUSEN, large-scale calculations can be achieved with GPU acceleration. In the physical model, particles are connected by solid bonds to simulate the failure behavior of layered rock mass during quasi-static compression, and the solid bonds act as interparticle cohesion.
For the development of BPM, scholars have made many contributions in theory and simulation. Based on the elastic assumption, many models can describe the stress field in a homogeneous spherical particle. Hertz calculated the pressure distribution in the contact area, the contact force, and the radius of the contact surface [18]. Huber describes the stress distribution in the half-space particle contact area and can be applied to calculate the contact volume of a sphere [19, 20]. Rumpf and Schönert calculate the spatial stress distribution for the entire sphere [20]. Antonyuk et al. outline various solutions for contact force-displacement behavior [21]. The mechanical failure behavior of particles is determined by the microscopic bonding mechanism, and the first systematic study of particle strength has confirmed the micro-macrocorrelation of particle bonds [22–24]. The tensile strength of particles is described by the well-known Rumpf [25] model, which takes into account the viscous force at the point of contact of a single spherical particle in a stochastic model. Bika et al. generalized the Rumpf [25] and Kendall [24] models and generalized them to porous particle structures [26]. The results show that the tensile strength of particles with a highly viscous liquid binder and high saturation depends on the capillary pressure [23]. Delenne et al. used the Mohr-Coulomb failure criterion to model and simulate the deformation and failure behavior of coupled cylindrical rods under different stress or loading modes (tension, compression, and shear)[27].
Numerically, the discrete element method is mainly used to study particle fracture. The motion of individual particles due to applied forces and moments in discrete time steps is calculated based on discrete element simulations. Originally, the discrete element method was developed for modeling ideal spherical particle populations. However, the method has been extended to model agglomerates composed of nonspherical objects or primary particles. To reproduce the shape and internal structure of the aggregates, BPM can be used effectively [28, 29]. BPM represents a particle as a set of initial particles connected by ideal elastic or viscoelastic solid bonds [28, 30–32, 32–34]. The results of the numerical simulations can provide detailed data on the external forces acting on the particles and the stresses inside the particles at the particle-scale, which cannot be obtained from experimental measurements. Similar to the BPM approach to study the mechanical properties of elastic particle composites, the approach of treating the material as a set of elastic springs can also be used. However, in contrast to BPM, this method cannot effectively describe the fracture of materials. BPM can also be used to analyze the effects of structural parameters that are difficult to separate experimentally. For example, the effect of solid bridge strength can be easily isolated in BPM, while the use of different binders in experiments may alter the internal structure of the particles, complicating the effect of solid bridge strength on agglomerate strength.
In experiments, many scholars have carried out relevant research on the tensile and compressive properties of layered rock mass. Liu et al.[35–37] summarized the effects of bedding and matrix on tensile and compressive strength by conducting indirect tensile tests, finite element, and discrete element numerical analysis of layered shale and schist. The results show that bedding has a significant effect on tensile strength and splitting failure. For layered slate, Liu et al.[38–41] revealed transverse isotropic strength characteristics and failure mechanism through uniaxial compression and Brazilian splitting tests, Ma et al. established a numerical model to simulate the failure process of shale specimens under triaxial compressive stress. The deformation and strength characteristics of the structure, the progressive failure process and the corresponding acoustic emission response are obtained by numerical simulation, and the rock failure process is reflected in detail [42–44]. Domestic scholars have carried out systematic experimental and theoretical research on gneiss, shale, bedded sandstone, etc.. Tavallalt and Vervoot conducted discrete element analysis of layered sandstone and found that the laboratory test results are consistent with the calculation results of rock failure patterns [45]. Youchang carried out three point bending fracture toughness tests under different loading rates on shale specimens with three typical bedding states, and systematically studied the fracture toughness development law of shale under different loading rates and bedding changes [46]. Gholami and Rasouli conducted Brazilian splitting and uniaxial compression tests on slate with different bedding inclinations, and obtained anisotropic characteristics of mechanical parameters such as failure mode, tensile strength, and compressive strength [47]. BPM has been widely used. In the discrete element method, each pair of particles can be connected by solid bonds. Dosta et al. define the constitutive relationship of chemical bonds, which can capture different macroscopic effects such as material hardening or softening [32, 33]. Zhao et al. have made an in-depth study of carbonaceous slate through tessellation polygon numerical simulation [48, 49]. The following points can be drawn from the analysis: the construction of different numerical models leads to obvious differences in numerical simulations. According to the numerical analysis of different bedded slate, the bedding structure has a great influence on the tensile and compressive strength of carbonaceous slate. In recent years, BPM has been increasingly used for 2D and 3D simulation and numerical analysis of rock materials, but BPM is rarely applied to layered carbonaceous slate.
The general framework of this paper is as follows: Firstly, the sample preparation, the principle of solid bonded particle model, the construction of bedding model, and the application of boundary conditions are described. Secondly, the relationship between the microscopic and macroscopic Young’s modulus was studied, the effects of normal tensile, and tangential shear strength on uniaxial compressive strength were investigated, and the effects of different models on the mechanical response of layered slate were analyzed. Finally, numerical simulations were carried out for samples with different bedding inclination angles and compared with the experiments (fracture mode and load-displacement response).
2. Carbonaceous Slate and Sample Preparation
The test samples were taken from the face of the construction tunnel of the Muzhailing Tunnel, and the rock mass at the sampling point had a thin-layered structure. The selected rock blocks have obvious bedding, and the bedding has different scales. The interlayer spacing is about 4.5-22.5 mm after rough measurement. The density of the sample is 2.688 g/cm3, and it is determined by X-ray diffraction that the sample is mainly composed of quartz and clay minerals. The uniaxial compressive strength of the horizontally bedded carbonaceous slate and the vertical bedded carbonaceous slate are about 49.2 MPa and 52.5 MPa. The horizontal elastic modulus and the vertical bedding elastic modulus are about 7.0 GPa and 7.8 GPa. The horizontal and vertical Poisson’s ratios were 0.20 and 0.23.
Figures 1 and 2 are the uniaxial compression and Brazilian splitting samples of layered carbonaceous slate, respectively. The uniaxial compression specimen has a bottom diameter of 50 mm and a height of 100 mm. The Brazilian split specimen is 50 mm in diameter and 25 mm in thickness. The sample processing mainly includes drilling cores, cutting cylindrical discs, and leveling samples. In order to reduce the influence of moisture on the rock mass properties, the prepared samples were allowed to stand for 60 d under ventilated conditions. Assuming that the included angle between the bedding plane and the horizontal is the bedding inclination angle, the layered carbonaceous slate samples under five inclination angles (0, 30, 45, 60, and 90) are used for the test. More specific mechanical properties and preparation process of the samples can be found in the literature [48, 49]. The test load-displacement response curve and the test results will be compared with the numerical simulation in Section 4 and Section 5. The next section mainly introduces the calculation principle of the numerical simulation software.
3. BPM Calculation Principle and Establishment of Physical Model
3.1. Particle-Particle Interaction
In the physical model, each pair of primary particles can be connected by one or more solid bonds. Each particle and each solid bond can possess unique geometric and material parameters to build complex structures of heterogeneous materials. In this study, soft spherical discrete elements are used, and the overlap between particles or between particles and boundaries can be explained as local material deformation. In general, three different types of interaction models are applied: (1) the interaction between the initial particles. (2) Interaction between particles and boundary geometry. (3) Solid bonds are used to describe the forces, moments, and torques acting on a single solid bond. Due to the nearly ideal elastic properties of carbonaceous slate, the linear elastic contact model is used for particle-grain contact, and the Hertz-Mindlin model is used for particle-boundary interface contact [34]. For particle-to-particle interactions, the forces are resolved in the normal and tangential parts, respectively. In the normal direction, is calculated as follows:
In the formula, is the normal tensile strength of the particles, is the coefficient of restitution between two particles, is the normalized contact vector between the centroids ( and ), is the overlapping volume, is the relative velocity between two particles, and is the equivalent mass of the two particles with mass and mass . In the tangent direction, is calculated as follows:
In the formula, is the tangential shear strength of the particle, is the tangential force of the previous iteration, is the increase of the tangential overlap at the current time step, and is the relative velocity of interacting particles. Furthermore, is limited by the interparticle friction coefficient , which is related to sliding friction
If equation (5) is not satisfied, then equation (3) is modified according to
In addition, the resultant force between particles is the sum of the normal force vector and the tangential force vector:
3.2. Solid Bond Model
The solid bond between the initial particles is an ideal cylindrical solid with initial length and radius. The radius cannot exceed the minimum radius of the contact pair, and the elastic bond has additional damping forces [28, 32] in the normal and tangential directions. One of the challenges of using the elastic bond model is the almost undamped oscillation during the simulation. To reduce this effect, artificial dampers were added. Generally speaking, the stress acting on the bonds is the result of the interaction between the corresponding particles. In order to simulate the fracture process of the material and the stress in the binder, a comparison is made with the tensile/compressive strength and the tangential strength of the material.
In the formula, and are the resultant normal and tangential forces on the solid bond, namely, the axial force and shear force at both ends of the solid bond, which can be expressed as . The sum of the particle-particle vector resultant force and the solid bond vector resultant force is the total force of the particle, namely, . In addition, the acceleration of particles can be obtained according to Newton’s second law of motion. and are the torsional and bending moments acting on the bond, and are cross-sectional area and radius, and and are torsional inertia and rotational inertia, respectively.
3.3. Construct Physical Model
The physical models at different bedding inclinations are shown in Figure 3, and the application of model boundary conditions is shown in Figure 4. Figure 4(a) represents the particle-solid bond model, Figure 4(b) represents the solid bond (green represents the matrix solid bond, red represents the bedding solid bond), and Figure 4(c) represents the particle model. In the model, the particle diameter is 1 mm, the porosity is 0.38, the number of particles is 300975, and the number of solid bond elements is 3141792. The model building method is divided into the following steps: (1) firstly, a geometric shape with a bottom diameter of 50 mm and a height of 100 mm is formed in MUSEN. (2) According to the diameter and porosity, corresponding particles and solid bonds are formed inside the geometry, and GPU acceleration can be used in the process of particle formation. (3) Extract the formed particles and solid bonds in text format, and select the equally spaced bedding solid bond elements through programming, and finally form the model shown in Figure 3.
(a)
(b)
(c)
There are two advantages to constructing the bedding model in this way: Firstly, it is guaranteed that the same model is used in the five bedding simulation processes, and its internal structure is the same. Secondly, it can be seen from Figure 4(b) that the bedding solid bonds are serrated in the model, which is closer to the rock mass material structure. Figure 4 shows the application of boundary conditions when the bedding is 30. The upper and lower plates are loading plates. In order to be consistent with the test conditions, the lower plate is fixed, and the upper plate is 0.1 mm/min. The speed moves down. In the same way, the Brazilian splitting model with different bedding inclination angles is constructed as shown in Figure 5. Among them, the particle diameter is 1 mm, the number of particles is 75264, the number of solid bond elements is 925667, and the total number of elements is 1000931.
4. Uniaxial Compression
4.1. Parameter Calibration
4.1.1. Calibration of Micro and Macro Young’s Modulus
Taking the homogeneous slate model as the specimen, the relationship between the micromodulus and the macromodulus is shown in Figure 6. Among them, the microscopic Young’s modulus represents the elastic modulus of the particles during numerical simulation, and the ordinate Young’s modulus represents the test modulus obtained by calculating the slope according to the load-displacement curve and converting, that is, the test Young’s modulus. In numerical simulation, the selection of microscopic Young’s modulus is very important. In order to correspond the Young’s modulus of the specimen to the microscopic Young’s modulus, we performed 24 sets of simulations for correction. The specific method is as follows: under the condition of the same geometric size, six models are constructed with the same particle diameter and porosity, and different microscopic Young’s moduli (2 GPa, 8 GPa, 12 GPa, and 16 GPa) are used for each group of models. After the model (corresponding to the four moduli) is calculated, the macroscopic Young’s modulus can be obtained from the stress-strain curve. The simulated boundary conditions are the same as in Figure 4, the lower steel plate remains unchanged, the top plate moves downward, and the speed is 0.1 mm/min. From the curve shown in Figure 6, it can be seen that the relationship between the micromodulus and the macromodulus is roughly linear. Therefore, when the Young’s modulus of the specimen is measured, the microscopic Young’s modulus can be obtained through interpolation calculation. For this uniaxial compression test, the Young’s modulus is about 7.0 GPa, and through parameter correction, the microscopic elastic modulus is 6.35 GPa.
4.1.2. Influence of Solid Bond Tensile and Shear Strength on Mechanical Properties of Rock Mass
As shown in Table 1, the ratio of shear strength to tensile strength is increased or decreased proportionally to study the effect of strength ratio on the mechanical properties of rock mass. The simulated boundary conditions are the same as in Figure 4, the lower steel plate remains unchanged, and the top plate moves downward at a speed of 0.1 mm/min. The final load-displacement curve is shown in Figure 7. According to the curve, it can be seen that with the gradual increase of the tensile shear strength of the solid bond, the corresponding load and displacement at the peak value increase. In addition, as the tensile shear strength gradually increases, the cohesion between the particles increases, resulting in an increase in the peak load. The speed should be selected as small as possible to keep the model in quasi-static loading to eliminate the influence of kinetic energy on the results. When the speed is too large, the load-displacement curve will rise intermittently in waves, and as the speed gradually decreases, it will gradually show a smooth curve. After testing, when the speed is 0.1 mm/min, the model is roughly in the quasi-static loading process. The tensile shear strength of the solid bond has a great influence on the load peak value. Therefore, in the selection of normal tensile strength and tangential shear strength, based on the above principles and through the trial and error method, the agreement between the numerical simulation and the test is enhanced.
4.2. The Effect of Particle Models
According to the principle of randomness, five different particle distribution models are formed to analyze the influence of different particle models on the mechanical response of rock mass. In the numerical simulation, the boundary conditions are consistent with Figure 4, the lower steel plate remains unchanged, and the top plate moves downward at a speed of 0.1 mm/min. Figure 8 is a schematic diagram of a model failure process under uniaxial compression. The corresponding load-displacement curves under five different particle model conditions are shown in Figure 9. It can be seen from the curves that although the particle models are different, the final mechanical response curves are roughly similar, but there are also certain influences. Therefore, in the comparison between this experiment and numerical simulation, in order to reduce the influence of the model, the same particle model is used for numerical calculation and analysis.
(a)
(b)
(c)
4.3. Comparative Analysis of Uniaxial Compression Test and Simulation
In the selection of the tangential strength and the normal strength of the bedding model, the shear strength of the rock material is greater than the tensile strength, so the value of the tangential strength is greater than the normal strength. The cohesiveness of the matrix is higher than that of the bedding, so the solid bond strength of the bedding is greater than that of the matrix solid bond. In the selection of normal and tangential strength, based on the above principles and through the trial and error method, the degree of agreement between the numerical simulation and the test is enhanced. The Table 1 is the solid bond normal and tangential strength parameters of the bedding model.
Figure 10 shows the comparison of the load-displacement curves between the numerical calculation and the test under different bedding inclination angles. The following points can be drawn: (1) In the test, the peak loads of 0°-90° are 98.374kN, 81.457kN, 78.673kN, 80.751kN, and 102.595kN in turn. As the inclination angle gradually increases, it roughly presents a “U”-shaped distribution. (2) Both the simulated curve and the test show the development trend of linear elasticity and instantaneous splitting and falling, and the sample exhibits obvious brittle splitting failure characteristics. (3) The simulation and test load-displacement curves are roughly consistent.
(a) 0° bedding
(b) 30° bedding
(c) 45° bedding
(d) 60° bedding
(e) 90° bedding
4.4. Fracture Mode
Figure 11 shows the failure process of samples and models under different bedding inclinations. According to the failure evolution, the following points can be drawn: (1)Figure 11(a) is the 0° bedding failure process. In the initial stage, the bedding solid bonds are broken and the matrix solid bonds perpendicular to the bedding are also broken gradually. With the accumulation of roof pressure, the model occurs along the matrix solid bonds, accompanied by a small amount of bedding solid bonds. The model eventually suffered split tensile failure(2)Figure 11(b) is a diagram of the 30° bedding failure process. In the initial stage, shear slip occurs on the bedding structure surfaces at both ends. As the pressure of the upper roof increases gradually, the middle structural plane is gradually damaged, and the splitting damage occurs along the normal direction of the bedding structure plane. The final model shows the staggered splitting failure of the matrix structure plane and the bedding structure plane(3)Figure 11(c) and Figure 11(d) are 45° and 60° bedding failure process. In the initial stage, the crack occurs with slip shear failure along the upper or lower end of the bedding structure plane. With the gradual increase of the load, the shear force gradually increases, and the intermediate bedding structure surface is gradually damaged. In the final model, a through crack is formed along the bedding structure plane, and the failure mode is shear slip failure(4)Figure 11(e) is a diagram of the 90° bedding failure process. In the initial stage, the matrix structure plane and the bedding structure plane at both ends of the rock mass fracture almost simultaneously. With the accumulation of time, the cracks undergo splitting and tensile failure along the bedding structure, eventually forming one or more through cracks(5)Through numerical analysis, the simulation results are close to the fracture mode of the numerical model. Therefore, the numerical model can better reproduce the splitting failure process
(a) 0° bedding
(b) 30° bedding
(c) 45° bedding
(d) 60° bedding
(e) 90° bedding
4.5. Effect of Solid Bond Force on Model Failure
Figure 12(a)–12(f) is a diagram of the failure process of the solid bond element inside the 0° bedding model, and the color bar represents the magnitude of the solid bond force. As can be seen that when the solid bond color changes to red, the bond at this position will be removed from the model according to the equation (8) and the equation (9). From the process of Figures 12(a)–12(c), it can be found that the bedding solid bonds are broken first. With the accumulation of time, as shown in Figure 12(d)–12(f), the solid bond force along the normal direction of the bedding gradually increases, and finally a through crack is formed. Figures 13(a)–13(f) is a diagram of the failure process of the solid bond element inside the 30° bedding model, and the color bar represents the magnitude of the solid bond force. As can be seen that the model is fractured along the edge bedding in the initial stage. With the accumulation of time steps, slip fracture occurs in the inward bedding, and the solid bonds are partially fractured along the normal direction of the bedding. Finally, the model is unstable in the slip-splitting mode.
5. Brazilian Split
5.1. Schematic Diagram of the Comparison of Load Displacement and Fracture Mode
Under different bedding inclinations, the test and simulation load-displacement curves and the comparison of fracture failure modes are shown in Figure 14. The following points can be drawn from the result analysis: (1) in the test, the peak loads from 0°-90° are 2603.92kN, 2242.22kN, 1767.58kN, 1562.03kN, and 1398.08kN. When the inclination angle changes from 0°-90, with the gradual increase of the inclination angle, the peak load decreases sequentially. After the peak, the load-displacement curves at different bedding inclinations drop instantaneously, and damage or softening can better explain its physical meaning. In the numerical simulation, when the load does not reach the peak value, the curve grows linearly and elastically, and then drops rapidly after the peak value. The comparison shows that the simulation is in good agreement with the experiment. (2) When the inclination of the bedding is 0°, the model undergoes tensile splitting failure along the line connecting the two loading ends through the matrix and the bedding. When the inclination of the bedding is 30° and 45°, the failure of the specimen roughly presents a broken line type failure. When the bedding inclination angle is 60°, the slip shear failure of the specimen occurs almost along the bedding, the end breaks along the matrix, and the cracks roughly show a broken line distribution. When the bedding inclination angle is 90° splitting failure occurs along the bedding, and the cracks are distributed in a straight line. Through the comparison of fracture modes, it can be seen that the numerical simulation failure form of the bonded solid model is roughly similar to the test.
(a) 0° bedding
(b) 30° bedding
(c) 45° bedding
(d) 60° bedding
(e) 90° bedding
5.2. The Process of Split Evolution
Figure 15(a) is a diagram of the 60° type bedding failure process. In the initial stage, shear slip failure occurred in the upper part along the bedding plane, and then the bedding in the lower part also began to fracture. With the increase of external force, the fracture mode of the model is mainly slip failure and splitting failure is supplemented. The final fracture mode is roughly fold-line or arc-shaped. Figure 15(b) is a diagram of the 90° type bedding failure process. In the initial stage, the crack gradually propagates through the center of the circle and the upper and lower ends. With the gradual increase of the external force, a penetrating crack appears in the middle, and finally the composite splitting tensile failure occurs and the splitting is distributed in a straight line.
(a) 60° bedding
(b) 90° bedding
In order to analyze the fracture process of the model from the magnitude of the solid bond force, a brief analysis of the bedding inclination angle 60° sample is carried out here. Figures 16(a)–16(h) is a diagram of the failure process of solid bond elements inside the 60° bedding model. It can be seen that in the initial stage, the stress in the upper part along the bedding plane first appeared in red. As the external force gradually increases, the disappearance of the upper red line means that the bedding solid bond is removed from the model. Over time, red solid bonds also appeared at the bottom layer. The fracture modes of the final model are roughly symmetrically distributed on a broken line or arc.
5.3. Advantages of 3D Numerical Software MUSEN
Through tension and compression test and numerical study, we can summarize the following characteristics: (a)Model construction. Figure 17 shows the sample and particle model under uniaxial compression. Figure 17(a) represents the slate sample, Figure 17(b) represents the three-dimensional particle model, and Figure 17(c) represents the two-dimensional particle model [48, 49]. Through the comparative analysis of the above two model figures, we can draw the following two points: (i) bedding structure. Compared with the two-dimensional model, the bedding structure of the three-dimensional model is more consistent with the actual sample. (ii) Through three-dimensional numerical analysis, the propagation state of the internal crack of the material can be reproduced, as shown in Figure 12, which is also the most important place and characteristic of this paper(b)Computational efficiency. The software MUSEN is used to study the mechanical behavior of layered carbonaceous slate. Although the sample is a small-scale specimen, as shown in Figure 17(b), the total number of elements in the numerical model is 3141792 (the number of particles is 300975). In terms of the number of particles, the size of the particle model is relatively large. When PFC3D is used to form a model with a number of 150,000 particles, the software will frequently exit in most cases. However, particle models can be stably and rapidly formed by CPU or GPU in MUSEN. Therefore, the numerical study of rock mass by software MUSEN provides a certain reference for engineering practice(c)Parameter correction. In this paper, a series of numerical simulations have been done for different particle models. On the one hand, the influence of particle model on numerical simulation is explored. On the other hand, the experimental Young’s modulus is transformed into microscopic Young’s modulus by parameter correction and interpolation method, which lays a foundation for layered numerical simulation, as show in Figure 6
6. Conclusion
The peak strength and fracture mode of slate are explored through tensile and compression tests and systematic numerical simulation of layered carbonaceous slate with different bedding inclination angles, and the following conclusions can be drawn: (1)Through parameter corrections to the macroscopic Young’s modulus and the microscopic Young’s modulus, a linear relationship is obtained(2)The construction of a sawtooth particle model in line with the actual material is of great significance to the numerical simulation results. On the one hand, the consistency between the numerical simulation and the fracture mode of the specimen is enhanced. On the other hand, the interaction between the particle model and the simulation parameters promotes the consistency on the load displacement between the numerical simulation and the specimen(3)The load-displacement curves of each group of layered slate have the development trend of linear elasticity-instantaneous splitting and falling, and the samples all show significant brittle splitting failure characteristics. In uniaxial compression samples, as the inclination angle gradually increases, the compressive strength decreases first and then increases, which roughly presents a “U”-shaped distribution. In the Brazilian split specimen, when the bedding inclination angle changes from 0° to 90°, the peak load decreases gradually(4)It can be seen from the fracture modes of the uniaxial compression specimens that with the change of bedding inclination, the specimens mainly exhibit sliding shear failure along the bedding plane and a combined tensile-shear failure mode affected by the matrix and bedding. For Brazilian split samples, the splitting failure modes of layered carbonaceous slate can be divided into split tensile failure, shear slip failure, and composite failure with the change of bedding inclination angle
The mechanical response and splitting failure form of layered slate based on BPM are in good agreement with the experimental comparison, and the splitting evolution process of layered slate under different bedding inclination angles is reproduced
Data Availability
The focus of this paper is to use the three-dimensional numerical software MUSEN to conduct in-depth research on the internal fracture characteristics of layered slate. On the copyright of the data, (a) the experimental data have been included in the article and are freely available. (b) For more specific experimental data, it is freely available in the following two articles, their doi are [10.1007/s12517-022-10630-4] and [10.1007/s12665-021-10106-w], respectively. On the other hand, these two articles have also been cited in this paper.
Consent
Informed consent was obtained from all individual participants included in the study.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
The present work is supported by the National Key Research and Development Program of China (No. 2016YFC0600901) and the National Science Foundation of China (Nos. 41172116 and U1261212).