Particle Flow Code Simulation of the Characteristics of Crack Evolution in Rock-Like Materials with Bent Cracks

The distribution and propagation of rock cracks have a signi ﬁ cant impact on geotechnical engineering. Taking rock-like materials with bent cracks as the research object, the particle ﬂ ow code in two dimensions numerical simulation method was used to study the impact of the bend number on rock-like materials strength and crack evolution. According to the results, when the bend number was 1, 3, and 7, the strength of the specimens gradually increased; the elasticity modulus did not change signi ﬁ cantly with the crack bend number. Uniaxial compression generated tensile cracks in all the specimens with bent cracks, but in terms of failure mode, the specimens with 0 bend tended to su ﬀ er penetrating failure along the fracture strike, while the specimens with 1, 3, and 7 bend tended to su ﬀ er penetrating failure along the diagonal direction. Both the fractal dimension and bend number were positively correlated with strain; with the gradual increase of the stress percentage, the damage variable of the specimens gradually increased at a growing rate. The research results provide a reference for predicting the stability of the underground engineering surrounding rocks containing bent cracks.


Introduction
Crack defect propagation has long been a focus of research for rock mechanics researchers [1][2][3]. For in-depth studies, scholars have performed a lot of advanced tests and pioneering work from different aspects by laboratory tests [4,5], numerical simulation [6][7][8], and theoretical derivation [9,10]. These experimental and numerical results enhance the understanding of the influence of joints on the mechanical behavior of rock masses.
The earliest research on crack propagation started with a single crack model. In 1920, Griffith pointed out that material damage was caused by crack propagation when external load led to stress concentration at the tip of these defects. Nolen-Hoeksema et al. [11] observed a failure mode in marble caused by crack propagation at the oblique crack tip, discovering that crack propagation was asymmetric, and pointed out that the propagation of penetrating cracks on specimen surface could well reflect the inside situation. Zhou et al. [12] conducted uniaxial compression tests on prismatic specimens with preexisting flaws and used the digital image cor-relation method to detect the process zone nucleation characteristics in granite; the complete cracking processes of granite were classified into six levels, which were distinguishable by five characteristic stresses. Saadat et al. [13] used a grain-based distinct element model to investigate the influence of grain size on the fracturing response of precracked granite, developed, and implemented a cohesive model in distinct element codes to mimic the elastic and softening response of the intragrain contacts in GBM. Wong et al. [14,15] conducted a uniaxial compression test on a rocklike material with preexisting cracks, finding that the inclination of preexisting cracks, the inclination of rock bridges, and the friction coefficient of cracks could affect the final failure mode of the specimens. Huang et al. [16] used laboratory tests and discrete element method simulations to investigate the influence of angle, spacing, joint length, and rock bridge length on uniaxial compressive strength, finding that the joint inclination angle has the most significant influences on uniaxial compressive strength and young's modulus. Yi et al. [17] proposed a new electrical method of conductive carbon-film to continuously measure crack propagation rate of brittle rock under coupling condition. A self-designed coupling testing system, finding that the higher the temperature is, the earlier the crack is initiated and the larger the crack propagation rate and accelerated velocity are. Zhao et al. [18] carried out numerical simulations of a rock specimen containing a single preexisting crack by the expanded distinct element method; the simulation results indicated that the crack inclination angle and confining pressure had a great influence on the tensile and shear properties, peak strength, and failure behaviors. Morgan et al. [19] conducted a series of unconfined compression tests on prismatic specimens with two preexisting flaws and various bedding plane orientations and used high-speed and high-resolution imagery to capture crack initiation; finding that in specimens with bedding planes aligned perpendicularly to the maximum applied compressive stress, the effects of the flaw angle and flaw pair bridging angle on the cracking processes were similar to those observed in previously tested rocks without bedding planes.
The previous studies focus only on neat cracks distributed in parallel or crosswise. However, cracks are usually distributed in an irregular sequence in actual engineering. Therefore, it is of great significance to study the influence of bent cracks on rock mass failure. In this paper, we study the influence of bent cracks on the rock strength and fracture evolution by Particle Flow Code in two dimensions (PFC 2D ) numerical simulation. First, we used the trial-and-error method to reveal the microparameters that could reflect macromechanical properties. Then, we simulated specimens with bent cracks to study the impact of different bend numbers on crack propagation. This study can provide a reference for the crack propagation of underground engineering with bending joints and cracks.

Numerical Modeling and
Parameter Verification PFC 2D is a piece of 2D simulation software based on discrete elements, commonly used to simulate the mechanical behavior of rocks and rock-like materials [20]. It uses rigid discs to simulate the particles in rock specimens and builds different connection models to simulate the connection among particles [21]. Interparticle friction is activated by judging interparticle bonded failure. The program can track and record the process of crack initiation and propagation and judge the type of cracks [22].

Numerical
Modeling. The uniaxial compressive strength of cement-gypsum specimens referred to in this paper is 3.91 MPa. Figure 1 shows the schematic diagram of the numerical simulation model. Its size is the same as that of laboratory test specimens, the cross-section size is 150 mm, and the height is 300 mm. As for the crack size, its projected length is 40 mm and width is 1 mm. To ensure sufficient crack smoothness and calculation speed, several attempts were made, finally realizing uniform distribution by controlling the particle radius at 0.5~0.9 mm. Parallel bonding model was chosen as the contact mode, and each specimen for numerical simulation contained 25,569 particles and 62,029 contacts. To systematically characterize cracks, open vertical cracks that are the same as one another in terms of projected length are made in this paper by removing particles to keep the linear length between both ends of a crack unchanged. n = 2 i − 1, (i = 1, 2, 3, 4), according to the calculation results, the crack bend number is 0, 1, 3, and 7, respectively. The geometric size and morphology of cracks are shown in Figure 1.
The walls on the upper and lower ends were given a constant speed by simulating the loading mode adopted in laboratory testing, and a load was applied to the specimen by controlling wall displacement. To study the effect of the loading rate, four loading rates, including 0.05 m/s, 0.1 m/s, 0.15 m/s and 0.2 m/s, were set.

Microparameter
Verification. The macroelasticity modulus is mainly linearly related to the contact modulus and bonding modulus of the particles, while Poisson's ratio is logarithmically related to the particle stiffness. The compressive strength is determined by the bonding stress ratio of the particles. Moreover, the ratio is mainly affected by the tangential bonding force of the particles when it is less than 2. The ratio is mainly affected by the normal bonding force when it is greater than 2, and it is hardly affected by the friction coefficient of the particles. In addition, the postpeak curve shape is affected by the friction coefficient.
In the parallel bonding model, parallel bonding between particles can not only transmit forces and vectors between particles but also transmit forces at contact points. Therefore, the parallel bond model is suitable for simulating the mechanical strength and fracture behavior of rock-like materials with bent cracks.
Based on laboratory test data of rock-like materials, granular flow simulation microparameters were obtained after repeated debugging, as shown in Table 1. Figure 2 shows the uniaxial compression stress-strain curve obtained by numerical simulation and laboratory test, and the overall variation trend of the two curves is similar. Table 2 presents the uniaxial compressive strength and the elastic modulus of calibrated specimen corresponding to Table 1. The deviation obtained for the uniaxial compressive strength and elastic modulus is 0.26% and 0.76%, respectively. The simulation results can better fit the laboratory test results.

Macromechanical Properties of
Specimens with Bent Cracks 3.1. Stress-Strain Curve Analysis. Figure 3 is the uniaxial compressive stress-strain curve of specimen with bent cracks at different bend numbers. When the bend number was 0, 1, 3, and 7, the peak stress was 3.24 MPa, 3.01 MPa, 3.45 MPa, and 3.75 MPa, respectively, while the elastic modulus was 0.309 GPa, 0.316 GPa, 0.319 GPa, and 0.325 GPa, respectively. As can be seen from Figure 3, the peak stress of the specimen with preexisting vertical cracks was 3.24 MPa, while the 2 Geofluids corresponding axial strain was 0.0129. In the OA segment, the stress increased from 0 to 64% of the peak strength, leading to linear elastic deformation of the specimen, with the elastic modulus of the specimen basically unchanged. In the AB segment, the stress was equal to 64%-100% of the peak strength. During this period, microcracks were initiated on the specimen surface and then propagated freely with inner cracks, evolving into penetrating cracks; when the stress increased to the peak strength value, the elastic modulus gradually decreased, with the specimen at the elastic-plastic stage. In the BC segment, cracks continued to propagate, and with the fracturing of the crack tip, the specimen ushered in an obvious yield stage, with the volume strain turning from compression to expansion. At this stage, the crack tip continued to expand, but unsteadily, and the rock specimen suffered increasing macrodamage along with crack propagation while a change took place in the stress state of the crack tip; the stress intensity rebounded but, then, fell again and changed in coordination with postpeak rock specimen damage. When the crack bend number was 7, it already   3 Geofluids rebounded before the peak point and, then, fell again after reaching the peak point. At this stage, crack surface friction hindered rock damage, manifested as periodic oscillation of postpeak strength. Meanwhile, the postpeak stress-strain curve of cracked rock masses fluctuated violently. In the CD segment, secondary cracks propagated quickly, so that cracks completely penetrated from top to bottom, with the specimen damaged and destabilized.
As can be seen from Figure 4, the peak stress gradually increased with the bend number from 1 to 7, indicating that as the bend number increased, the specimens were more inclined to generate shear damage on preexisting crack surface. The vertical cracks are parallel to the principal stress direction, which have little impact on crack propagation at the early stage. Moreover, in the loading process, the interparticle force is hindered from spreading between the two sides of vertical cracks.
The elastic modulus changed slightly with the bend number and the difference less than 4.9%, suggesting that the change in the bend number did not lead to an obvious change in the elastic modulus. This is primarily determined by the contact modulus and parallel bonding modulus of particles set during modeling.

Loading Rate and Mechanical
Properties. The loading rate has an important effect on the mechanical properties and deformation of cracked rock masses, so it can be changed to study the mechanical properties of specimens with bent cracks at different strain rates. The strain rate is calculated as follows: where ε * is the strain rate, s -1 ; ε is strain; Δl is the length variation of specimens, m; l is specimen length, m; t is loading time, s; and v is the loading rate, m/s. Since the specimen length was 0.3 m, the strain rate corresponding to the loading rate 0.05 m/s, 0.1 m/s, 0.15 m/s, and 0.2 m/s was set to 0.167 s -1 , 0.333 s -1 , 0.5 s -1 , and 0.667 s -1 , respectively. Figures 5 and 6 show curves of peak strength and elastic modulus changing with the strain rate at different bend numbers. As shown in Figure 5, with the rise of the strain rate, the peak strength of the specimen increased, and when the bend number was 0, 1, 3, and 7, it increased by 4.98%, 6.94%, 7.32%, and 6.71%, respectively. As shown in Figure 6, the elastic modulus was insensitive to the change of the strain rate, with errors of 0.65%, 1.28%, 1.27%, and 0.93%, respectively, indicating that the cracked rock-like specimens mainly showed rock-like mechanical properties at the elastic stage.
Linear function σ max = aE + b and power function σ max = a + bE c were used separately to fit the above results, as shown in Table 2. Figure 5 shows the strain-rate-dependent curve of peak stress fitted in a power function. As can be seen from Table 3, fitting by the power function is more in line with the relationship between the peak strength and strain rate than the linear fitting.

Law of Crack Propagation in Specimens with Bent Cracks.
To better reveal damage and failure mechanism of the specimens, we made a detailed analysis of the process of crack initiation and propagation at different bend numbers. For an easy comparison, the strains at the periodic transition points of the four were taken as special points to record the process of crack initiation and propagation Figures 7-10 show the process of crack initiation and expansion at different bend numbers. As can be seen from the figure, the final crack mode of the specimens is affected by the bend number. As shown in Figure 7, when the bend number was 0, the cracks had little impact on the fracture mode, but the crack width and distribution increased. The initial cracks were per-pendicular to the crack loading surface and parallel to the force loading direction. In this layout form, the tensile stress was obviously concentrated at the crack tip, tearing the cracks along the side direction. Under the condition of uniaxial compression, new cracks were formed on the specimen surface owing to the internal dilatancy effect of the specimens. As loading continued, the plastic zone expanded, with the open cracks propagated continuously, generating shear cracks near the tensile cracks.
When there was only 1 crack bend, as shown in Figure 8(a), ε = 0:0067, a damage zone appeared at the Vshaped crack tip first. When ε = 0:0122, the tensile cracks propagated quickly at the V-shaped crack tip, as shown in Figure 8(b). The specimen was under the peak stress, with two obvious tensile cracks formed in its surface. With continuous crack propagation, new tensile cracks were initiated near the V-shaped cracks, and as time went by, the upper and lower V-shaped cracks got interconnected, forming a triangular broken area. As shown in Figure 8(c), when ε = 0:0183, the specimen was in a postpeak stress state. At this stage, the specimen stress decreased slowly, the surface damage gradually worsened. As shown in Figure 8(d), when ε = 0:0235, the residual stress was equal to 51.2% of the peak stress. With the propagation and penetration of macrocracks, the cracks were compacted by the load, with almost all of the preexisting cracks compacted and the macrocracks penetrating through the specimen in the direction of thickness, causing the splitting failure.
When the bend number was 3, as shown in Figure 9(a), when ε = 0:0069, there were microcracks formed at the crack tip. As shown in Figure 9(b), when ε = 0:0122, with a stress peak appearing. Tensile cracks were formed near the preexisting cracks, and the shear slip band continued to develop, with new tensile cracks formed near the shear cracks. As shown in Figure 9(c), when ε = 0:0183, the specimen was in a postpeak stress state, while the macrocracks propagated and penetrated, with tensile cracks and penetrated through the shear zone around the preexisting cracks. As shown in Figure 9(d), when ε = 0:0235, tensile failure occurred at the crack tip, and the cracks penetrated through the specimen surface, depriving the specimen of its bearing capacity.
When the bend number was 7, as shown in Figure 10(a), ε = 0:0067, with microcracks initiated throughout the specimen. As shown in Figure 10(b), when ε = 0:0122, the specimen was in a peak stress state, with tensile cracks formed near the preexisting cracks, the shear slip zone continuing to develop, and the new tensile cracks formed around the shear cracks. As shown in Figure 10(c), when ε = 0:0183, the specimen was in a postpeak stress state, while the macrocracks propagated and penetrated, with tensile cracks and penetrated through the shear zone around the preexisting cracks. As shown in Figure 10(d), when ε = 0:0235, cracks propagated throughout the specimen, and the stress dropped to 1.74 MPa. The shear slip zone near the WW-shaped cracks improved the bearing capacity of the specimen, giving it a high residual stress.
The local damage of the specimen became more obvious in the loading process. Due to the impact of preexisting cracks, the microcracks on the specimen surface almost 5 Geofluids occurred near the preexisting cracks. After the penetration of macrocracks, the specimen suffered a postpeak stress. At this time, there remained some residual strength, and new cracks were formed near the initial macrocracks. Tensile cracks were formed at the crack tip first, or tensile cracks were formed owing to the impact of stress concentration near the loading end. Most of the macrocracks were tensile cracks, accompanied by shear slip. In the process of crack propagation, secondary cracks were formed and interpenetrated through one another, causing the generation of macrocracks, finally getting the specimen damaged.
According to the simulation results, the failure mode of specimens containing different members of crack bends can be divided into 4 stages, including crack initiation, crack propagation, overall rupturing, and instability due to failure. As for I-shaped cracks, because they are parallel to the principal stress direction, their final impact fell only on the specimen strength, causing penetrating failure along the crack trend. For V-shaped, W-shaped, and WW-shaped cracks, the shear-slip band that appears near the preexisting cracks caused diagonal damage on the specimen.

Fractal Characteristics of a Specimen with Bent Cracks.
The fractal dimension of the surface topography of the specimen could be calculated based on particle images. Then, a corresponding MATLAB program was developed and used for batch processing of images with characteristic significance. After gray-scale processing, binaryzation, and box  7 Geofluids dimension calculation, the fractal dimension was figured out. Fractal can be used to characterize the crack propagation path. The change of fractal dimension indicates the change of damage degree. The fracture crack propagation of rocklike materials can be characterized by the fractal dimension of the surface cracks. The larger the fractal dimension is, the more serious the damage is.
As can be seen from Figure 11, on the whole, the fractal dimension was positively correlated with the bend number and strain. As the strain increased, the fractal dimension corresponding to each bend number increased from 1.671, 1.683, 1.705, and 1.712 to 1.865, 1.823, 1.834, and 1.856, respectively. The fractal dimension was divided into three sections along with the growth of strain, and the slope of the three sections increased first and then decreased.
In the simulation process, we recorded the bend number within a certain range of specimen changing with the step of time and then used the cumulative number of microcracks to characterize the damage variable. As shown in Figure 12, with the gradual rise of the stress percentage, the fractal dimension of the specimen under uniaxial stress gradually increased. In the meantime, the damage variable of the specimen gradually increased at a growing rate. It increased obviously when the fractal dimension rose to about 1.75 and later increased faster. According to the analysis results, as the stress percentage rose, damage gathered at the damaged parts of the specimen, forming macrocracks. Before the load reached the peak stress, the microcracks inside the specimen 8 Geofluids gradually evolved into macrocracks. In this process, the fractal dimension reached a relatively fixed value, while the damage variable increased rapidly. As the bend number increased, the change of the damage variable gradually became smooth, and with the increase of the fractal dimension, damage gathered inside the specimen, i.e., the microcracks gradually evolved into macrocracks. As the load increased, the damage variable increased at a decreasing rate, and its distribution became increasingly sparse, but when it was about to reach the peak, the damage variable increased faster. This shows that at the later stage of specimen deformation and failure, the increase of the same axial stress would cause greater damage to the specimen.

Conclusions
Taking rock-like materials as the research object, the PFC 2D numerical simulation method was used to study the impact of the crack bend number on rock strength and crack evolution. The main conclusions are as follows: (1) As the crack bend number increases from 1 to 7, the specimen strength increases linearly in a gradual way, but the vertical crack strength is lower than that of W-shaped cracks (2) The peak strength of the rock masses with bent cracks is positively correlated with the strain rate while there is an insignificant correlation between the elastic modulus and the strain rate (3) In the loading process, tensile cracks and shear cracks are formed in all specimens with bent cracks, of which the specimen with 0 bend cracks tend to suffer penetrating failure along the fracture strike while the specimens with 1, 3, and 7 bends tend to suffer penetrating failure along the diagonal direction (4) The fractal dimension is positively correlated with the bend number; with the gradual rise of the stress percentage, the damage variable of the specimens gradually increases at a growing rate

Data Availability
The 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.