Investigation of Mechanical Characterization of Transversely Isotropic Rock Mass Based on Rigid Block Discrete Element Method

The laminated joints of transversely isotropic rock mass show strong discontinuity, inhomogeneity, anisotropy, and nonlinearity. This paper proposes a novel rigid block discrete element method (RB-DEM) to simulate the transversely isotropic rock mass. The balls in traditional DEM simulation are replaced by rigid blocks generated using FEM mesh. The contact on the laminar joints can be described accurately. The rigid block discrete element method (RB-DEM) results are consistent with the laboratory test. The e ﬀ ects of joint parameters, joint spacing, and sti ﬀ ness ratio on the uniaxial compression properties of transversely isotropic rock mass are discussed. This study provides a novel and e ﬀ ective tool for the analysis of transversely isotropic rock mass.


Introduction
The transversely isotropic rock mass is widespread in geological materials (as shown in Figure 1). Isotropy refers to the property wherein an object's physical and chemical properties do not change with its orientation [1,2]. A transversely isotropic rock mass is a particular case of anisotropic rock masses. A rock mass with the same elastic properties in all directions in a plane is called an isotropic surface. In contrast, the mechanical properties in the direction perpendicular to this surface are different, and the rock mass with such properties is called a transversely isotropic rock body. Researchers have studied the lamellar structure of phyllite in-depth and pointed out that rocks containing phyllite structures can be reduced to transversely isotropic mass. The transversely isotropic elastomer model describes the deformation and damage analysis of natural geological bodies formed by sedimentation like this. Due to the specificity of its engineering performance, it is often encountered in the stability analysis of slopes and underground works. Therefore, researchers must study and clarify the damaging evolution of transversely isotropic rock mass under compressive conditions [3][4][5][6].
Transversely isotropic phenomena are relatively common in geological materials, and anisotropic traits have a crucial influence on the stress-strain analysis and damage mechanical behavior of rock mass. Scholars have done much research work at home and abroad. The current study mainly includes theoretical, experimental, and numerical simulations. For academic studies, Sakurai et al. proposed the displacement-strain feedback analysis method based on the retrospective analysis method to derive the complete initial state of stress and Young's modulus from a set of relative displacements measured between adjacent monitoring points [7]. Fan and Ye established the three-dimensional equation of state in the local coordinate system, obtained the accurate static and dynamic solutions of the three-layer orthogonal thick plate, and proposed the orthotropic elastic theory [8]. Valanis and Haupt proposed an internal time plasticity model based on the theorem of irreversible thermodynamics [9]. These theoretical studies provide a crucial theoretical basis for clarifying the mechanical behavior of transversely isotropic rock mass.
In terms of experimental research, Amadei studied and verified the relationship between five elastic parameters and stress-strain of transversely isotropic rock mass under different bedding angles [10]. Nguyen et al. recorded the process of inclined cracks in Naples tuff under uniaxial compression based on the digital image correlation method [11]. Based on fracture mechanics, Zhou et al. established a damage model considering the random distribution of microcracks in the rock mass. They studied the influence of microdefects on the strength and stiffness of rock mass [12]. Tien and Tsao implemented a split Hopkinson pressure bar (SHPB) test on artificial layered rock specimens with different dip angles [13]. The relationship between Young's modulus, overall strength, rock deformation, and dip angle was studied using the incomplete bonding interface composition model (IBICM) [14].
For numerical simulation, Chen et al. proposed the random field simulation using the Karhunen-Loève (K-L) expansion [15]. Faizi et al. conducted triaxial numerical simulations based on the discrete element method for different dip angle models and investigated the mechanical behavior of transversely isotropic rock mass [16]. Zhao et al. calculated the damage zone at the laminate tip based on the Kachanov equation and discussed the crack connections between the laminate surfaces. The transmission coefficients of the stress waves at the level are derived [17].
Because the cross-sectional isotropic rock is more fragmented, it is difficult to collect the experimental specimens in the field, and it is equally difficult to produce multiple rock specimens in the laboratory. The development of numerical calculations has made it possible to study the compressional damage processes and damage mechanisms in the transversely isotropic rock mass. Traditional computational models do not simulate the crack formation process, and now, we usually use discrete element models or granular flow models to simulate transversely isotropic rock mass problems. The limited number of computational units in the granular flow model has limitations in finely modeling the mechanical processes of progressive failure of the intermittent transversely isotropic rock mass. In contrast, the discrete element model can consider the cracking of through joints and the movement of blocks, using efficient solving methods to find and classify contacts, enabling further improvements in simulation efficiency. This paper proposes a rigid block DEM modeling method (RB-DEM). Compared with conventional DEM, this method has the advantages of simple model generation, high computational efficiency, and full consideration of the influence of the joint layers. In addition, using the DEM model allows for a more accurate simulation of crack extension in transversely isotropic rock masses. Adjacent boundaries no longer use the SJ smooth joint model but instead use particle-to-particle contact to transfer forces [16]. By setting different strength parameters for the blocks and joints, we can simulate the columnar transverse isotropic rock properties more realistically. By comparing and verifying that the numerical results are close to the field test results, the paper discusses the influence of various factors on the uniaxial compression properties of the rock mass.   2 Geofluids interaction between particles and particles or between particles and boundaries, as well as the macroscopic motion of the whole discrete system. Rigid blocks, which can simulate closed polygons, flow polygons, and irregular polygons, are introduced into the DEM calculation. In this paper, the balls in traditional DEM simulation are replaced by rigid blocks generated using an FEM mesh [18]. In a DEM, rigid blocks can be created by countless methods, such as creating rigid blocks from vertices or importing rigid blocks from geometry. As shown in Figure 2, this paper uses a closed polygon rigid block with a shape similar to the rigid block model proposed by Francesco et al. Compared with the traditional method, the contact detection algorithm of this rigid block has better robustness and can accelerate the computation [19,20]. The discrete element solves Newton's laws of motion using explicit differential dynamics, and its program uses an efficient solution method to detect and classify contacts. With contact detection, contact will be generated if a pair of rigid blocks are likely to collide. Once contact is created, the algorithm will apply different contact detection algorithms depending on the type of object in contact at both ends. The program can maintain data structures and memory allocations, so tens of thousands of contacts can be processed quickly.

Adhesive Properties of Transversely Isotropic Rock
Mass. As shown in Figure 3, the linear parallel bonding model is used as the linear model between the contacts in this paper [21][22][23]. It is a linear-based model that can be installed at both ball-ball and ball-facet contacts. The linear parallel bond model provides the behavior of two interfaces: an infinitesimal, linear elastic (no-tension), and frictional interface that carries a force and a finite-size, linear elastic, and bonded interface that carries a force and moment (see Figure 3). The first interface is equivalent to the linear model: it does not resist relative rotation, and slip is accommodated by imposing a Coulomb limit on the shear force. The second interface is called a parallel bond because, when bonded, it acts in parallel with the first interface. When the second interface is bonded, it resists relative rotation, and its behavior is linear elastic until the strength limit is exceeded and the bond breaks, making it unbonded. When the second interface is unbonded, it carries no load. The unbonded linear parallel bond model is equivalent to the linear model. The model can reasonably simulate the mechanical properties of selected objects with parallel bonded members that transfer forces and moments between rigid blocks and do not prevent sliding from occurring [24].
The contact interaction law that defines the relationship between contact forces and contact overlap properties is essential for DEM. As shown in Figure 4, in this paper, when the contact members are all rigid blocks located in the same group, the contact type is defined as a rock contact. When the contact members are rigid blocks located in two different groups, the contact type is defined as a joint contact. The effect of various joint parameters on the strength of the nodular rock mass is discussed in the next section.

Rigid Block Generation.
The size of the geometric model in this paper is 100 mm × 50 mm, as shown in Figure 5  3 Geofluids mesh, the model is preprocessed, as shown in Figure 5(b). In order to avoid abnormal meshes from interfering with DEM calculations, the corresponding software is used to generate well meshes in the model. All mesh elements are triangles, creating a structural mesh of well-shaped cells in Figure 5(c). The next step is to import the model into FLAC3D, list the blocks for each group of rock masses, and convert them to the corresponding block collection. In Figure 5(d), we can see the final rock mass model, where each rigid block consists of a set of blocks transformed from many blocks.

Parameter Calibration.
Due to the influence of various complex factors like the mineral composition of the rock mass, the development level of joints, and the orientation of joints, obtaining a suitable set of rock parameters becomes extremely difficult. The microscopic parameters selected in the DEM cannot be determined directly from the physically tested rock properties and need to be calibrated by an iterative trial-and-error approach. This paper utilizes the pairwise algorithm proposed by Thurstone for solving psychometrical problems to find the most suitable parameters [25]. Assume that each dimension is orthogonal. When β = 90°, since the rock matrix mainly bears the applied load in this direction, the weak joint layer has the least influence on the effective modulus and rock strength. Therefore, it is used to calibrate the effective modulus of the experimental rock sample close to Zhang, the cohesive force, and tensile strength [26]. The parameters of the experimental rock samples are approximated by changing the effective modulus, cohesion, and tensile strength. Several iterations of the above steps are performed to determine the best combination of effective modulus, cohesion, and tensile strength. When β =15°, 30°, 45°, 60°, and 75°, the damage mainly occurs on 4 Geofluids the weak joint layer of the body. Therefore, the effective modulus, cohesion, and tensile strength of the rock matrix are kept unchanged, and the parameters of the experimental rock samples are approached by changing the effective modulus, cohesion, and tensile strength of weak joints. Several iterations of the above steps are performed to determine the optimal combination of the effective modulus, cohesion, and tensile strength for weak joints. When β = 0°, the normal stiffness ratio becomes the fundamental factor that determines the difference between the stress-strain curves at β = 0°and β = 90°since the damage occurs mainly in the rock matrix of the rock mass. The effective modulus, cohesion, and tensile strength of the rock matrix and weak joints can be kept constant to approximate the stress-strain curve of the experimental rock sample by changing the normal stiffness ratio. Several iterations of the above steps are performed to determine the normal stiffness ratio of the weak joints and the rock matrix. Compared with the traditional orthogonal design, pairwise optimization can reduce the input cost and increase efficiency. In addition, the larger the number of dimensions, the more significant the results, which applies to this paper. Table 1 lists the detailed parameters of the rock mass model in pairs [27].

Mechanical Characterization of Transversely Isotropic
Rock Mass. In this paper, a total of seven uniaxial compression test simulations were performed to investigate the effect of seven anisotropic angles (β = 0°, 15°, 30°, 45°, 60°, 75°, and 90°) on the properties of the transversely isotropic rock mass. During each simulation, the change in stress value is determined by the synthetic stress acting on the wall, while the change in strain is calculated from the wall displacement [28,29].
After the rock is loaded, the strain increases with the rise of stress. When the pressure increases to the rock strength value, the rock fails. The stress-strain curve of solid magmatic rocks, limestone, and other rocks are nearly linear, showing elastic deformation and sudden brittle damage, which is in good agreement with the simulation results of this paper [26].    Figure 6. Comparing the experimental and simulation results, we can see that the two have a high degree of fit, generally showing a V shape.
With the change of inclination angle, the general trend of the compressive strength curve is low in the middle and high on both sides. The curve decreases in the dip angle of the rock body from 0°to 60°, and the curve is relatively flat.
The inclination angle of the rock mass decreases from 0°to 60°, and the curve is relatively gentle. The strength of the rock sample is the largest at 0°and the smallest at 60°. The inclination angle of the rock mass rises within 60°~90°, and the curve is relatively steep. From Figure 6, it can be seen that when the dip angle of rock joints β = 90°, the corresponding peak strength is 50 MPa, while when out when the dip angle of rock joints β = 60°, the peak strength of the rock is 20 MPa. The ratio of the corresponding peak strength is 2.5. The gap between the two is too large. It reflects that the dip angle of the rock mass will have a great impact on the strength of the rock mass. In engineering practice, we should pay more attention to the influence of rock mass inclination on engineering safety [30].

Failure
Process. The data obtained during the simulation test were organized and plotted into a stress-strain curve. Figures 7 and 8 show the stress-strain curves of the uniaxial compression test of the columnar transversely isotropic rock mass with a dip angle of 0°to 90°. According    [31]. It can roughly be divided into five stages: microcrack compaction and closure stage, linear elastic deformation stage, microcrack initiation and propagation stage, volume expansion and crack unsteady propagation stage, and failure and postpeak deformation stage. In the early stage, the original cracks and pores inside the sample gradually closed with the increase in load. At this stage, the stress-strain relationship of the rock appears as a concave upward curve, but no internal microcracks are established in the simulation model, so the crack closure stage does not appear in Figure 7 [32 -34]. Once most of the primary cracks are compacted and closed, the deformation of the rock enters the linear elastic deformation stage, and the stress-strain curve shows a straight line segment. The linear elastic deformation stage of rock is a stage of microcrack initiation. When the stress level applied to the rock exceeds the sprouting stress threshold σ ci for new cracks, the microcracks inside the rock start to sprout and show stable expansion [35][36][37]. At this stage, the microcracks are in a steady state of expansion and the crack damage threshold. When the stress level exceeds the crack damage threshold σ cd , the microcracks inside the rock sample begin to penetrate gradually, form a crack network, and finally form a macrofracture surface. The crack damage threshold corresponds to the inflection point on the volume strain curve where the rock volume shifts from compression  Geofluids to expansion. The specimen reaches peak strength after σ p, and the rock sample is damaged. Many new microcracks in the rock are generated, expanded, and merged, failing the specimen. For uniaxial compression tests, the peak rock strength is σ p . It is an important and commonly used mechanical index in rock engineering. When the rock sample exceeds the peak strength σ p , the stress-strain curve of the rock enters the postpeak softening stage, and the stress decreases with the increase of the strain [38,39].

Failure Mode.
The damage modes of transverse view isotropic rock masses in uniaxial compression tests are shown in Figure 9. As the column changes from rotation too, the loading direction leads to anisotropic damage modes of transverse view isotropic rock masses. Overall, the damage modes of uniaxial tests can be divided into three categories [40][41][42]: splitting damage perpendicular to the column axis, shear slip damage, and splitting damage along the column axis. When the column inclination is perpendicular to the loading direction, for example, β = 0° (Figures 9(a) and 9(b)), the columnar transversely isotropic rock mass mainly suffers from splitting failure along the axial direction of the specimen. Most anisotropic rock samples show split failure under uniaxial compression, which is the characteristic of uniaxial compression failure. When the inclination of rock mass is medium or steep, for example, β = 30°, 45°a nd 60° (Figures 9(c)-9(e)), shear slip along the inclined rock mass face is the predominant structural damage mode. For steeper dips, shear stress plays a key role in driving the shear slip damage along the middle of the rock mass rock specimen. When the axis of the rock mass is approximately parallel to the loading direction, for example, β = 75°and 90° (Figures 9(f) and 9(g)), the transverse-view isotropic rock mass exhibits another type of cleavage damage. In this dip range, the axial stress causes the transverse-view isotropic rock mass to undergo cleavage damage along the nearvertical or vertical polygonal rock face, resulting in the collapse of the rock mass.

Effect of Joint Parameters.
A rock mass is a combination of rock blocks and structural planes. The strength of the rock mass is the comprehensive strength considering the effect of the structural plane. Factors such as the occurrence, joint length, and connectivity of the joint plane have an important influence on the strength of the rock mass. Since the essence of the joint plane is the weak plane in a specific direction developed in the rock mass, the strength of the structural plane is weaker than that of the rock block. Joint planes are usually divided into general joint planes and weak joint planes. The strength of the weak structural face is usually in the range of 0~0.02 MPa, while the strength of the general structural face is usually greater than 0.03 MPa. In summary (as shown in Table 2), the strength ratio of the joint plane and rock block is selected as 0.25, 0.5, 0.75, and 1 to study the influence of the joint surface on the strength of rock mass with different dip angles.
In Figure 10, we can see the stress-strain curves of rock masses with different joint strengths under uniaxial compression. When the joint strength is constant, with the change of the dip angle, the overall trend of the rock mass compressive strength curve is that the middle is low, and the two sides are high. When the dip angle is in the range   9 Geofluids of 0°~60°, the curve decreases, and the strength of the rock mass sample is the smallest at 60°. When the dip angle is in the range of 60°~90°, the curve rises, and the rock mass sample has the highest strength at 90°. When the dip angle is constant, as the joint strength increases, the peak strength of the rock mass generally increases. When the dip angle is 0°and 90°, the change in the joint strength leads to a tinier change in the rock mass strength. When the dip angle is 15°, 30°, 45°, 60°, and 75°, the alteration in the joint strength leads to a larger change in the rock mass strength. The reason for the above results is that when β is 0°and 90°, the main failure mode of the rock mass is splitting failure, and the strength of the rock mass depends on the strength of the inner rock block. When β is 15°, 30°, 45°, 60°, and 75°, the main failure mode of the rock mass is shear slip failure along the joint plane, and the rock mass strength is more affected by the joint strength. Figures 11(a)-11(d) show the stress-strain curves obtained by simulating the rock mass with different joint strengths under uniaxial compression conditions. As can be seen from the figure, as the joint strength changes, the stress-strain curve also varies. It can be roughly divided into three stages. In the early stage, the rock enters the linear elastic deformation stage, and the stress-strain curve shows a straight line segment. With the increase of the joint strength, the deformation gradually increases, and the two show a good linear relationship. When the stress exceeds the elastic deformation stress limit, the stress-strain curve deviates from the elastic deformation curve. During this stage, the rock gradually enters the yield stage. Since the joint spacing is different, the yield strength is different. With the increase of the joint spacing, the yield strength also increases accordingly. When the stress exceeds the ultimate stress, the stress decreases, the strain softens, and the stress-strain curve  10 Geofluids shows the residual strength. Even at this stage, joint strength has a great influence on residual strength. The residual strength increases as the joint spacing increases. Therefore, it has influences on the three stages of the stress-strain curve. The elastic deformation, yield strength, and residual strength all increase with the increase of joint strength.

4.2.
Effect of Joint Spacing. The rock mass is a geological material full of various discontinuous surfaces, and joint surfaces are typical representatives of such discontinuous surfaces. The strength of joints is influenced by the mechanical properties of rocks and joints on the one hand and by the geometric characteristics of joints (such as arrangement, yield, and spacing) on the other hand.
The researchers conducted many physical simulation experiments to study the correlation between the performance of transversely isotropic rock mass and the geometric parameters of the joint network.
According to the research results of Einstein and Hirschfeld and Prudencio and Jan on the strength and deformation characteristics of rock masses with different joint spacings [43,44], this paper selects three joint spacings of 5 mm, 20 mm, and 50 mm (as shown in Figures 12(a)-12(c)) to study the influence of the joint spacing on the strength of rock mass with different dip angles.
In Figure 13, we can see the strength curves of rock masses with different joint spacings under uniaxial compression. When the joint spacing is constant, with the change of the inclination angle, the general trend of the compressive strength curve is that the middle is low, and the two sides are high. When the dip angle is in the range of 0°~45°, the curve decreases, and the strength of the rock mass sample is the smallest at 45°. When the dip angle is in the range of 45°~90°, the curve rises, and the rock mass sample has the highest strength at 90°. When the inclination angle remains unchanged, as the joint spacing increases, the peak strength of the rock mass generally increases. The reason for the above results is that when the joint spacing increases to a certain extent, the strength of the rock mass is equal to the strength of the rock block. When the joint spacing reduces to a certain extent, the rock mass strength is equal to the joint strength throughout. Figures 14(a)-14(c) show the stress-strain curves obtained by simulating the rock mass with different joint spacing under uniaxial compression. As can be seen from the figure, as the joint spacing changes, the stress-strain curve also varies. It can be roughly divided into three stages. In the early stage, the rock enters the linear elastic deformation stage, and the stress-strain curve shows a straight line segment. With the increase of the joint spacing, the deformation gradually increases, and the two show a good linear relationship. When the stress exceeds the elastic deformation stress limit, the stress-strain curve deviates from the elastic deformation curve. During this stage, the rock gradually enters the yield stage. Since the joint spacing is different, the yield strength is different. With the increase of the joint spacing, the yield strength also increases accordingly. When the stress exceeds the ultimate stress, the stress decreases, the strain softens, and the stress-strain curve shows the residual strength. Even at this stage, the joint spacing has a great influence on the residual strength.
The residual strength increases as the joint spacing increases. Therefore, the joint spacing influences the three stages of the stress-strain curve. The elastic deformation, yield strength, and residual strength all increase with the increase of joint spacing.

4.3.
Effect of Normal-Tangential Stiffness Ratio (k n /k s ). The discrete element method (DEM) replaces the objects of a discrete system with discrete units of a specific shape, and it focuses on the contact interaction between particles and particles or between particles and boundaries, as well as on the macroscopic motion of the whole discrete system [45,46]. The determination of the fine apparent parameters is very important when using discrete element models for calculations and plays a decisive role in the accuracy of discrete element numerical simulations. For the apparent parameters, we can obtain them through indoor and outdoor tests. We determined the mesoparameters by observing the macromechanical features characterized by the interaction of the particle set and the interaction of the microscale-modeled components. The complexity of particle interactions in discrete elements has led to macroscopic parameter relationships that have not yet been established. Therefore, it is very important to study the relationship between macroand fine parameters of rock mass for us to use a discrete element model to simulate the mechanical properties of the rock mass.
The researchers conducted a large number of physical experiments to study the mesoscopic parameters k n /k s of the rock mass. According to Rojek et al. research on rock mass strength and deformation characteristics under different k n /k s values. In this paper, four different ratios of 0.3, 0.5, 1.0, and 1.5 are selected to study the influence of different values of k n /k s on the strength of the columnar transversely isotropic rock mass [47].

Geofluids
In Figure 15, we can see the uniaxial compression strength curves of the rock mass under different stiffness ratios. When the stiffness ratio is constant, with the change of the inclination angle, the general trend of the compressive strength curve is that the middle is low, and the two sides are high. When the dip angle is in the range of 0°~45°, the curve decreases, and the strength of the rock mass sample is the smallest at 45°. When the dip angle is in the range of 45°~90°, the curve rises, and the rock mass sample has the highest strength at 90°. When the column inclination angle is constant, with the increase of the stiffness ratio, the peak strength of the jointed rock body generally shows a trend of first increasing and then decreasing. Figures 16(a)-16(d) show the stress-strain curves obtained by simulating the rock body subjected to uniaxial compression under different stiffness ratios. As can be seen from the figure, as the stiffness ratio changes, the stressstrain curve also varies. It can be roughly divided into three stages. In the early stage, the rock enters the linear elastic deformation stage, and the stress-strain curve shows a straight line segment. With the increase of the stiffness ratios, the deformation gradually increases, and the two show a good linear relationship. When the stress exceeds the elastic deformation stress limit, the stress-strain curve deviates from the elastic deformation curve. During this stage, the rock gradually enters the yield stage. Since the stiffness ratios are different, the yield strength is different. With the increase of the stiffness ratios, the yield strength also increases accordingly. When the stress exceeds the ultimate stress, the stress decreases, the strain softens, and the 12 Geofluids stress-strain curve shows the residual strength. Even at this stage, the stiffness ratio has a great influence on the residual strength.
The residual strength increases as the stiffness ratio increases. Therefore, it influences the three stages of the stress-strain curve. The elastic deformation and the residual strength increase with the increase of the stiffness ratio, while the yield strength decreases with the increase of the stiffness ratio.

Conclusion
In this paper, rigid block DEM simulations are performed on transversely isotropic rock mass models with different joint parameters, joint internal friction angles, and joint spacing under uniaxial compression. The study focuses on crack propagation and stress-strain curves for a transversely isotropic rock model and draws the following conclusions: (1) The joint strength has an influence on the three stages of the stress-strain curve. The elastic deformation, yield strength, and residual strength all increase with the increase of joint strength (2) When the joint spacing increases to a certain extent, the strength of the rock mass is equal to the strength of the rock block. When the joint spacing reduces to a certain extent, the rock mass strength is equal to the joint strength throughout (3) The stiffness ratio has an influence on all three phases of the stress-strain curve, and the elastic deformation and the residual strength increase with the increase of the stiffness ratio, while the yield strength decreases with the growth of the stiffness ratio

Data Availability
The data used to support the findings of this study are included within the article, and further data or information required is available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.