Strength and Failure Properties of Preflawed Granite under Coupled Biaxial Loading and Unloading Conditions

In this study, numerical simulations of uniaxial compression, biaxial compression, and biaxial unloading were performed on granite specimens that contained di ﬀ erent prefabricated defects. The microscopic parameters in numerical models were veri ﬁ ed by the uniaxial compression experiments on the intact standard cylindrical granite specimen and the square granite specimens with prefabricated defects. The in ﬂ uences from di ﬀ erent stress paths, di ﬀ erent shapes of prefabricated defects, di ﬀ erent numbers of defects, and di ﬀ erent distribution of defects on the strength, deformation, and crack initiation stress characteristics of the rock specimens were investigated. Furthermore, the initial cracking and cracking stage distributions, cumulative crack amounts, ultimate failure modes, and crack propagation fractal dimensions of specimens with di ﬀ erent prefabricated defects under biaxial unloading conditions were analyzed and compared. The experiment was divided into three stages to analyze crack evolution mechanisms. The results show that most cracks appeared after peak strength, and di ﬀ erent shapes, the number of defects, and the relative defect positions signi ﬁ cantly a ﬀ ected crack initiation, crack propagation, and crack coalescence.


Introduction
As rocks form, various primary fractures and joints in natural rock are present due to influences from various complex geological environments, and existing primary fractures and joints affect rock strength and deformation characteristics [1]. Disturbances from underground engineering openings break the original stress distributions and intact rock mass causing defects already existing in the rock to expand and intersect each other and resulting in the stress redistribution, which eventually leads to rock failure [2][3][4]. To avoid possible damage from destabilization in underground engineering construction and ensure the safety of human life, property, and production, it is of the utmost importance to study various mechanical characteristics and crack propagation laws of rocks containing defects under different stress paths.
Many scholars have conducted a large number of indoor experiments on rock specimens containing defects and investigated the effects from different angles, different num-bers of fissures, different combinations of holes and fissures, and different stress states on rock crack expansion patterns and mechanical properties [5][6][7][8][9][10][11]. Haeri [12] investigated the effects of crack inclination angle and crack length on the fracturing processes of the Brazilian disk specimen with a precrack through central straight. Yang et al. [13,14] performed uniaxial compression on sandstone specimens containing single elliptical holes with different angles and double elliptical holes at different relative positions and analyzed the associated effects on deformation characteristics and crack initiation, propagation, and coalescence. Zeng et al. [15] produced sandstone specimens with differently shaped holes, such as hexagonal, square, and trapezoidal, to study the effects of differently shaped holes on the strength and crack characteristics of the sandstone under uniaxial loading conditions. Li et al. [16] combined double fissures with circular holes and conducted uniaxial compression tests to study expansion patterns for wing and secondary cracks. Yao et al. [17] conducted the uniaxial and conventional triaxial compression experiments on intact, single, and double fractured granites and compared their mechanical properties.
Other scholars studied crack coalescence mechanisms in rock-like materials containing double fissures [18][19][20][21]. Li et al. [22] studied the unloading responses for preflawed rock specimens under different unloading rates. Prudencio and Jan [23] conducted conventional triaxial tests on nonpenetrating defected rock materials and found that the failure mode depended on several factors, including the defect shapes, principal stress directions, and ratios of the intermediate stress to the compressive strength of intact material.
The discrete element method (DEM) is a numerical method based on the assumption of discontinuity, which is particularly suitable for solving discontinuity problems in fractured rock mass and also plays an important role in studying crack propagation [24][25][26][27][28][29][30]. Haeri et al. [31] investigated the effect of normal load on the crack propagation from pre-existing joints using particle flow code (PFC) 2D software. Zhang et al. [32] used a numerical model to simu-late the acoustic emission properties of rock-like materials containing single fissures at different loading rates. Jin et al. [33] simulated the effects of single fissures at different angles on crack initiation, deformation fields, and energy fields of rock-like materials under uniaxial compression. Yang and Huang [34,35] used PFC 2D software to simulate sandstone and rock-like specimens containing two unparallel fissures and analyzed the stress and displacement fields and failure modes under uniaxial loading. Zhao and Zhou [30] simulated the mechanical properties and failure modes of rock materials containing open and filled fissures under uniaxial compression. Wang et al. [36] studied the effects of normal loads and dip angles of prefabricated fissures on crack propagation under uniaxial compression. Liu et al. [37] used a density-based clustering method to detect microcrack clusters, and the well-known statistical tool of Principal Component Analysis (PCA) was applied in an innovative method in this experiment to define the macrocracks. Specimens containing two cavities were investigated for their normal load effects one crack propagation. Xu  Figure 1: Excavation-induced loading and unloading stress conditions, including the situations (a) before excavation, (b) after circular roadway excavation, (c) variations in stress near the circular roadway as the distance from the center of the roadway increases, and (d) variations in stress near the circular roadway with time (σ r is the radial stress, σ θ is the tangential stress, r is the distance between the center of the roadway and the rock element, a is the radius of the circular roadway, θ is the angle of the line between the rock mass and center of the circular roadway and horizontal line, p 0 is the hydrostatic in situ stress, and t is the working face advanced time). et al. [38] established a DEM numerical model for coal samples with two prefabricated fissures, and a numerical test of biaxial compression was conducted on fractured coal samples. When the effects from the Rock Bridge Dip Angle (RBDA) >30°, the number of tensile and shear cracks decreases as RBDA increasing. In addition, high confining pressure can promote increase in peak deviatoric stress. Yang et al. [39] investigated the failure behavior in both inner and outer zones around a circular opening in a non-persistently jointed rock mass with biaxial compression applied through numerical simulations. Sagong et al. [40] studied the rock fracture and joint sliding behaviors of jointed rock masses with an opening under biaxial compression through numerical analyses.
Voyiadjis et al. [41] developed an anisotropic elastoplastic-damage model for quasi-brittle materials within the thermodynamics framework. The damage criterion can capture larger axial damage in uniaxial tension and smaller axial   damage in uniaxial compression as compared with the lateral damage. Latha and Garaga [42] presented the results from elasto-plastic numerical simulations of jointed rocks using both the equivalent continuum and discrete continuum approaches. Nazerigivi et al. [43] checked the pure mode I (tension) and pure mode II (shear) by compressing the Central Straight-through Crack Brazilian Disk (CSCBD) specimens with central cracks at an angle of 0°and 27°with respect to the loading axis, and Sabri et al. [44] prepared some Single-Edge Cracked Brazilian Disks (SECBD) Granite specimen and CSCBD rock-like specimens and studied specimen fracturing processes.    Figure 3: The TRW-3000 true-triaxial testing apparatus, including (a) uniaxial testing and (b) an overview of the testing system.

Lithosphere
Before excavations for underground openings, the rock mass has an in situ stress state applied to it. Then, stress redistribution will occur when conducting the opening excavation. As the surrounding rock gets closer to the opening surface, the radial stress decreases, and tangential stress increases. Finally, the rock mass may suffer from unloading failure. As shown in Figure 1, for the rock in the excavation damage zone around the opening, the stress paths include radial stress decreasing, and the tangential stress first increases before decreasing when the stress reaches its peak, causing rock failure to occur. However, previous studies pri-marily focused on uniaxial compression applied to preflawed rock, while not focusing on coupled biaxial loading and unloading conditions (lateral stress decreasing and axial stress increasing). Therefore, it is necessary to investigate the strength and failure characteristics of rock specimens that contain different defect combinations under coupled biaxial loading and unloading conditions. In this study, a set of micromechanical parameters in the PFC was calibrated, and the systematic numerical simulations that isolated coupled biaxial loading and unloading were conducted for specimens containing single and combined  5 Lithosphere defects by simulating results for uniaxial compression as well as biaxial compression. Finally, the mechanical properties and crack coalescence processes under coupled biaxial loading and unloading conditions were analyzed in detail. The main novelty of this paper is that the biaxial unloading stress paths explored in this paper have never been investigated before, and the locations of the unloading point are derived jointly from the results of uniaxial compression strength and biaxial compression strength. Double-defect specimens with different relative positions studied in this paper have rarely been mentioned in previous studies, so it is necessary to investigate crack propagation and mechanical properties of specimens containing said defects.    The cylindrical specimen was tested in uniaxial compression with a displacement-controlled loading speed of 0.01 m/s. After a series of repeated tests, a set of suitable microscopic parameters (Table 1) was obtained to establish numerical granite specimens. The particle contact modulus and parallel-bond modulus are both 20 GPa, and the ratios of normal to shear stiffness of the particle and the parallel bond are both 1.7. The particle friction coefficient is 0.5, and the mean value and standard deviation of the parallel-bond   Table 2 shows the comparison of macroscopic mechanical properties for standard cylindrical granite specimens between indoor tests and numerical tests. From Table 2, the uniaxial compressive strength, Young's modulus, and Poisson's ratio of the standard cylindrical specimens are approximately equal under the numerical simulated and experimental values. The error values of experimental uniaxial compressive strength and Young's modulus from the simulated values are 0.07% and 0.66%, respectively. Figure 2 shows the comparisons of experimental and numerical stress-strain curves for standard cylindrical granite specimens, as well as comparisons of failure modes. It can be seen from Figure 2 that the two stressstrain curves are well fitted. However, the numerical model cannot simulate the initial compression-density phase of the stress-strain curve, resulting in some discrepancies between the final peak strain and experimental value. Both the experimental and numerical models exhibit the axial cleavage failure mode. Furthermore, two square granite specimens with dimensions of 100 × 30 × 100 mm were made, and the specimens were carefully polished in sections. Two square granite specimens and the standard cylindrical specimen were removed from the same intact granite blocks with a homogeneous matrix. One fissure and a single circular hole were formed by cutting the square specimens with a waterjet. The TRW-3000 true-triaxial testing apparatus, which was designed and manufactured by Central South University, performed the uniaxial compression tests on two types of square granite specimens, as shown in Figure 3. The PCI-II acoustic emission system was employed with a GTMicro-300 sensor, and the displacement was recorded using the TRW-3000 displacement sensor. The axial force was transmitted through the loading lead block for uniaxial loading.

Validation of Microparameters.
The principal stress σ 1 was continuously increased at the loading speed (0.5 kN/s) until the rock specimens ultimately failed. As shown in Figure 4, the crack propagation and failure mode of the granite specimens were similar to those that resulted in the numerical specimens. The mechanical properties of the numerical specimens are quite close to those of the real granite specimens, and the curves are in good agreement, so the numerical granite specimens can reproduce the mechanical behaviors observed in the physical granite specimens.

Specimen Model.
After the microscopic parameters were calibrated, the square numerical specimens of 100 mm × 100 mm were shaped, and the defects of different shapes were formed by deleting spherical particles. The complete square numerical specimen contained 25,920 round particles and 62,275 contacts. The complete square numerical specimen and the prefabricated defects specimen are shown in Figure 5 once the particles were deleted. The defect dimensions were a fissure width of 2 mm, fissure length of 16 mm, fissure inclination of 45°, and circular hole diameter of 10 mm. The natural rock in mining, tunneling, and petroleum engineering commonly have many defects, such as openings, fissures, and joints. These defects will develop and connect with each other under excavation-induced stress redistribution conditions. To investigate development characteristics of single defects and the connection mechanism among defects, different combinations of holes and fissures were considered in specimen models. The specimen models with different prefabricated defects and relative positions of different prefabricated defects are described in detail in Figure 6. Figure 1, the stress in the surrounding rock around the opening will change with the radial direction decreasing and tangential direction increasing. Thus, the coupled biaxial loading and unloading process is a common stress path that the surrounding rock around underground opening can be 9 Lithosphere subjected to. Meanwhile, the uniaxial compression and biaxial compression should be conducted to obtain the uniaxial compressive strength and biaxial compressive strength of rock specimens to determine the stress levels before unloading. The three stress paths were used in this test: uniaxial compression, biaxial compression, and coupled biaxial loading and unloading, as shown in Figure 7. The loading and unloading procedures were written to apply uniaxial compression, biaxial compression, and coupled biaxial loading and unloading to the numerical specimens in three stress paths in sequence. In uniaxial compression, σ 1 was loaded at a loading speed of 0.01 m/s until the rock failed. In biaxial compression, σ 1 and σ 2 were first increased to 20 MPa simultaneously, while σ 2 was kept constant at 20 MPa, and σ 1 was continuously increased at a loading speed of 0.01 m/s until the rock failed. The preloading process of coupled biaxial loading and unloading is the same as that of biaxial compression, but when σ 1 reached half of the uniaxial compression strength and biaxial compression strength sum combined, the servo program was rewritten to unload σ 2 at a rate of 1 MPa/ms, and σ 1 continued increasing at a loading speed of 0.01 m/s until rock failure to simulate the excavation-induced stress change path in the surrounding rock of the underground opening. The numerical software PFC2D was calculated in the simulation by time step, and the single test time for this simulation was in the tens of milliseconds for a given loading speed of 0.01 m/s. In this numerical simulation, the time step is equal to 0:65 × 10 −8 s/step, and 0.01 m/s can be converted to 6:5 × 10 −8 mm/step. The study shows that the loading speed of 0.2 m/s in PFC2D is capable for quasi-static loading of the rock specimen, and the loading speed of 0.01 m/s in the test fully satisfies the quasi-static loading [45][46][47]. Table 3 lists the strength and deformation parameters of numerical specimens containing different prefabricated defects under three types of stress conditions. The peak value of axial stress σ 1 during uniaxial compression is defined as the uniaxial compression strength (σ ucs ). The strain value corresponding to σ ucs is the peak uniaxial compression strain (ε u ), and the axial stress σ 1 when the number of cracks reaches 1% of the total is defined as the initial crack stress for uniaxial compression (σ uci ). The peak axial stress σ 1 in biaxial compression is defined as the biaxial compressive strength (σ bcs ). The strain value corresponding to σ bcs is the peak uniaxial compression strain (ε b ), and the axial stress σ 1 when the number of cracks reaches 1% of the total number of cracks is defined as the initial crack stress for biaxial compression (σ bci ). Since the procedure of coupled biaxial loading and unloading before the unloading process is the same as that of biaxial compression, and the initial cracks appear before the unloading procedure is triggered, the initial crack stress of biaxial compression and biaxial unloading are the  10 Lithosphere same. Half of the σ ucs and σ bcs sum combined is defined as σ ðucs+bcsÞ/2 , and the peak axial stress σ 1 under coupled biaxial loading and unloading is defined as the peak biaxial unloading strength (σ bu ). When σ 1 reaches σ ðucs+bcsÞ/2 during biaxial unloading, the unloading procedure is triggered, and the lateral stress σ 2 is unloaded at a rate of 1 MPa/ms from the previously maintained level of 20 MPa. When σ 1 reaches the biaxial unloading strength, the corresponding peak strain is defined as the peak biaxial unloading strain (ε bu ), and the corresponding lateral stress σ 2 is defined as the peak biaxial unloading lateral stress (σ ls ). The strength and deformation characteristics of numerical granite specimens containing different prefabricated defects under the three types of stress paths were compared, as shown in Figure 8. In Figure 8(a), for the same numbered specimens, the biaxial compression peak strains, the biaxial unloading peak strains, and the uniaxial compression peak strains are in decreasing order. The specimens exhibit the greatest variation in peak biaxial compression strains and the least variation in peak biaxial unloading strains. Overall, the peak specimen strains with fewer prefabricated defects were larger than those of specimens with more prefabricated defects along the same stress paths. The numerical biaxial compression results for specimens containing two fissures and specimens containing a single fissure and a single circular hole show that peak strain tends to decrease as the defect center moves, further away from the specimen central axis.

Strength and Deformation Properties.
The results of coupled biaxial loading and unloading were similar to biaxial compression. The numerical uniaxial compression results for the specimens containing two fissures and the specimens containing a single fissure and a single circular hole show that the peak strain of the specimens tended to increase and then decrease as the defect centers moved further away from the central axis of the specimen. However, the variation in peak strain for the specimens containing two fissures was small, which is consistent with the variation in uniaxial compressive strength in Figure 8(b). Figure 8(b) compares the compressive strengths under different stress paths, the numerical values of the uniaxial compressive strength, biaxial compressive strength, and biaxial unloaded compressive strength decreased as the number of prefabricated defects increased. The overall numerical results for the specimens containing two fissures and the specimens containing a single fissure and a single circular hole shows a decreasing trend in compressive strength as the defect centers moved further away from the central axis of the specimen. For the same numbered specimen, the biaxial compression strength, the biaxial unloading strength, and the uniaxial compression strength are in decreasing order. When the distance between the center of the prefabricated defects and the central axis of the specimen was fixed, the compressive strength of the specimens with two fissures was generally less than that of the specimens with one fissure and a single circular hole for all three types of stress paths.   The initial crack stress of the specimen under the coupled biaxial loading and unloading at the same number was larger than the uniaxial compression initial crack stress, indicating that initial crack generation under biaxial compression generally occurred later than uniaxial compression, and the lateral stress limited initial crack generation. The trend exhibited by uniaxial compression initial crack stress between specimens was similar to that of initial crack stresses under coupled biaxial loading and unloading, with both decreasing as the number of prefabricated defects increased. The initial crack stress results under the coupled biaxial loading and unloading show a decreasing trend as the defect centers moved further away from the central axis of the specimen. The uniaxial compression initial crack stress results show that the initial crack stress of the specimens with two fissures tend to decrease and then increase as the defect centers moved further away from the central axis of the specimen, while the initial crack stress of the specimens with one fissure and a single circular hole tend to decrease as the defect centers moved further away from the central axis of the specimen. From Figure 8(d), it can be observed that the peak lateral stress under coupled biaxial loading and unloading increased and then decreased as defect centers moved fur-ther away from the central axis of the specimen, but the variation in the specimens with two fissures was significantly larger than that of the specimens with one fissure and a single circular hole.

Cracking Characteristics under Coupled Loading and
Unloading. Figure 9 summarizes the axial stress-time curves and lateral stress, crack numbers, cumulative crack numbers of granite specimens simulated by PFC2D under biaxial unloading. The linear parallel bond model includes two mechanical behaviors of the linear model and parallelbond model. The linear model can provide an elastic mechanical behavior but cannot withstand tensile force, and its force-displacement law is where F l is the linear force, consisting of normal force F l n and shear force F l s , k l is the stiffness, consisting of normal stiffness k l n and shear stiffness k l s , and Δδ indicates the relative displacement between particles, including normal displacement Δδ n and shear displacement Δδ s , and then determine whether the contact is broken:  Figure 14: The crack coalescence process and particle contact force of specimen c-1 under coupled biaxial loading and unloading.

Lithosphere
Particle activity occurs where μ is the friction coefficient when the shear force exceeds the friction limit.
In the parallel-bond model, which can resist relative rotation and withstand certain tensile forces, the parallelbond part breaks when the applied stress reaches the strength limit, and the parallel-bond model degenerates to a linear model, as shown in Figure 10(a), and its forcedisplacement law is where the parallel-bond force F consists of parallel-bond normal force F n and parallel-bond shear force F s , k is the parallel-bond stiffness, consisting of parallel-bond normal stiffness k n and parallel-bond shear stiffness k s ; A is the bond area, the length direction is used as the diameter of the smallest particle D min , and the width is utilized as the thickness unit for the vertical plane amount, and then determine whether the bond model is broken: τ > τ c = c + σ tan shear failure: Where σ T is the tensile strength, τ c is the shear strength, c is the cohesion in the parallel-bond model, and φ is the angle of internal friction, as shown in Figure 10(c). When the parallel-bond normal stress σ > σ T , tensile cracks occur in the parallel-bond, and when the parallel-bond shear stress τ > τ c , shear cracks occur in the parallel-bond. When the normal stress or shear stress exceeds the corresponding parallel-bond normal strength and parallel-bond shear strength, tensile cracks, and shear cracks are generated, respectively, and the number of cracks and crack development process can be recorded and identified by fish language.
The effects of the different shapes, numbers, and relative positions of defects on the crack evolution process of the specimens were analyzed. A program was written to record crack distribution at each stage of the loading process. The calculation results show that the cumulative number of cracks grew slowly before unloading and increased faster from when unloading began to before peak strength, while the growth rate for the cumulative number of cracks after peak strength was significantly higher than that before peak strength, which showed that the primary growth period of the cumulative number of cracks was located in the postpeak damage stage after peak strength. The number of shear cracks where the specimens ultimately failed was significantly larger than that of tensile cracks, indicating that specimen failure was dominated by shear failure. Compared to specimens containing a single defect and two defects, the rapid rise over time of the cumulative crack number for specimens with nine prefabricated defects was shorter, the  Figure 15: The crack coalescence process and particle contact force of specimen l-2-1 under coupled biaxial loading and unloading.
13 Lithosphere peak phase of the number of cracks was more compact, and the curve of the total number of cracks rose more dramatically. Meanwhile, the ultimate numbers for cumulative cracks in specimens l-9 were greater than that of specimen c-9, but the first crack of specimen l-9 was initiated later. For specimens with two prefabricated defects, the ultimate cumulative number of total cracks, shear cracks, and tensile cracks in the specimen was lower, and the first crack was initiated earlier as the prefabricated defect center moved further from the center of the specimen. The ultimate numbers of cumulative cracks in the specimens with a single defect was significantly higher than that in specimens with nine prefabricated defects. During the period before peak strength, 16% of the ultimate cumulative crack numbers were generated in specimen l-1, while that of specimen c-1 remained at a low level, and the growth rate was significantly lower. The ultimate cumulative number of cracks after the rock failure was similar for both specimens, while the first crack of specimen l-1 was initiated later.
The ultimate damage modes for different prefabricated defect specimens under coupled biaxial loading and unloading conditions are listed in Figure 11. For the results of specimens with a single fissure and those with a single circular hole, the ultimate failure of both specimens was caused by crack coalescence on the left and right ends of the defect and the far-field cracks around the specimen. However, cracks near the left side of the defect in specimen l-1 extended upward, the cracks near the left side of the defect in the specimen c-1 extended downward before finally connecting with the cracks on the left corner of the specimen. The right side of the fissure in specimen l-1 developed a macroscopic fracture in the upper and lower directions. For specimens with two fissures, the wing cracks at the outer sides of the two fissures of specimen l-2-1 coalesced with each other, extending toward the bottom of the specimen. While the wing cracks at the inner sides of specimen l-2-2 and specimen l-2-3 coalesced, the crack density in the bridge area increased with the relative distance between the two fissures. For specimens with one fissure and a single circular hole, the cracks near the fissure and cracks near the circular hole of specimen c-2-1 did not coalesce directly, and no macroscopic cracks initiated from the right side of the circular hole. The wing cracks on the left side of the fissure of specimen c-2-2 eventually coalesced with the macroscopic cracks in the top right corner of the circular hole. While the wing cracks on the left side of the fissure of specimen c-2-3 eventually coalesced with the macroscopic cracks in the bottom right corner of the circular hole, a macroscopic tensile crack extended downward along the axial stress direction at the coalescence. For specimens with nine fissures, wing cracks near the fissure tips along the 45°line between the bottom left and top right corners of the specimen coalesced through each other. Similar to the results for specimen c-9, most of the cracks near the defects gradually coalesced with each other, while no obvious crack propagation behavior was observed near the fissure tips at the 3 Figure 16: The crack coalescence process and particle contact force of specimen l-2-2 under coupled biaxial loading and unloading.
14 Lithosphere upper and lower sides of the specimen central axis, the fissure at the lower end of the central axis of specimens l-9 produced two elongated cracks at the middle and lower right sides. The circular hole at the specimen c-9 center produced an X-shaped crack coalescence mode. The fractal dimension of crack propagation is an index parameter to quantitatively describe the complexity inherent to the crack propagation path. The crack propagation becomes more complex as the fractal dimension grows. This paper used the box counting method to calculate the fractal dimension D of crack propagation. The crack propagation map is covered with a square lattice of dimension r, the number of lattices containing cracks is counted and noted as NðrÞ, and the corresponding NðrÞ is recorded by varying the size of r. Taking the logarithm of r and NðrÞ, a linear fit is made to the computed data, and the fitted straight line slope is the fractal dimension D. The calculation formula is By analyzing crack propagation in the ultimate failure mode of the specimens, the fractal dimensions of crack propagation for specimens from l-1 to c-9 were calculated as 1. 588, 1.541, 1.561, 1.538, 1.482, 1.602, 1.598, 1.556,   1.536, and 1.464, as shown in Figure 12. The results indicate that the fractal dimension of crack propagation for the specimens containing two fissures tends to increase and then decrease as the defect center moved further away from the central axis of the specimen, while that of the specimens containing a single fissure and a single circular hole tends to decrease. In general, the fractal dimension of specimens containing a single defect was the largest, and that of specimens containing nine defects was the smallest, indicating that the specimens with a single defect produced the most complex secondary cracks during the failure process. It can be seen that different numbers of prefabricated defects, different relative positions of the defects, and different defect combinations significantly affected the fractal dimensions. Figures 13-22 recorded the crack coalescence and particle contact force processes over six stages during coupled biaxial loading and unloading. In the figures, the red points indicate cracks, and the brighter color represents the higher particle contact force. These specific points correspond to the points on the stress-time curve in Figure 10. We divided the process into three stages, namely Stage I: before peak strength, Stage II: at peak strength, Stage III: after peak strength. Figures 13 and 14 exhibit the crack coalescence processes and particle contact forces of specimens containing one fissure and specimens containing a single circular hole under  Figure 17: The crack coalescence process and particle contact force of specimen l-2-3 under coupled biaxial loading and unloading.
15 Lithosphere coupled biaxial loading and unloading. At Stage I before peak strength, the wing cracks were initiated from the right tip of the fissure and extended downward in specimen l-1 with a single fissure, and then the left tip of the fissure subsequently propagated a wing crack that extended upward. For specimen c-1, the cracks were initiated from the opening boundary and then extended along the diagonal of the circular hole. The particle contact force was larger near the fissure tips and the diagonal boundary of the circular hole. At Stage II at peak strength and Stage III after peak strength, the farfield cracks around the specimen gathered from the outside to the inside, gradually connecting with the cracks near the fissures and holes in the specimens. Figures 15 and 17 show the crack coalescence processes and particle contact forces for specimens containing two fissures under coupled biaxial loading and unloading. At Stage I, the wing cracks were initiated from the left tip of fissure A and the right tip of the fissure B of specimen l-2-1, then the right tip of fissure B subsequently propagated a vertical tensile crack, and the wing cracks on the left tip of fissure A continued to extend in an approximately vertical direction. The particle contact forces were greater near the two fissure tips, but the particle contact forces near the fissure tops and bottoms were weaker. Tensile cracks with directions perpendicular to each fissure were initiated nearby the right tips of two fissures in specimen l-2-2. Wing cracks initiated on the left tip of fissure B propagated toward the right tip of fissure A and gradually connected through. The particle contact forces were greater near the two fissure tips and in the bridge region and became greater with the increase in axial stress. For specimen l-2-3, the wing cracks were initiated from the inner tips of fissure A and fissure B. The particle contact forces were larger near the two fissure tips and in the bridge region, but the particle contact forces near the middle of fissure A and fissure B were weaker. The wing cracks at the two fissure A tips propagated continuously along the direction of the axial stress, and then connected with the crack propagation that roughly followed fissure B.
At Stage II when the axial stress reached peak strength, several cracks near the left tip of fissure A propagated and clustered towards the left tip of fissure B in specimen 1-2-1, while cracks near the right tip of fissure B propagated and clustered towards the right tip of fissure A. The farfield cracks around the specimen gathered from the outside while migrating to the inside. The propagation of far-field cracks in specimen l-2-2 was not obvious at this point. The wing cracks initiated from the left tip of fissure B extended downward to connect the wing cracks from fissure A, and the particle contact force in the bridge area became stronger. As for the specimen l-2-3, the wing cracks of fissure B propagated to the bridge area between the two fissures, and the wing cracks on the right end of fissure A extended  Figure 18: The crack coalescence process and particle contact force of specimen c-2-1 under coupled biaxial loading and unloading. 16 Lithosphere downward along the direction of the axial stress. The stress concentration area of the particle contact force in the bridge area became narrower, the particle contact force in the area below fissure A became significantly smaller, and the blank area without any contact force became larger. When the axial stress reached peak strength, fissure A and fissure B grew connected between each other overall. At Stage III after peak strength, the cracks near the two fissure tips in specimen l-2-1 finally coalesced with each other, and the far-field cracks in the top right corner propagated toward the fissure A and fissure B line. The far-field cracks in the bottom right corner gradually changed from their initial expansion direction to the vertical direction. The wing cracks near the left tip of fissure A changed its expansion direction to the vertical direction after peak strength coalesced with the far-field crack on the left side of the specimen. The particle contact force was significantly reduced in the failure area and was greater on the left and right sides near the fissures. For specimen l-2-2, the wing cracks initiated near the left tip of fissure A propagated in the vertical direction, and the cracks at the two tips of the same fissure clustered with each other. The far-field cracks on the left side of the specimen did not coalesce with fissure A, while the far-field cracks on the right side of the specimen developed toward the right tip of fissure B and aggregated with each other. The particle contact force was concentrated at the fissure tips, and the nearby area and the area without contact force rapidly became larger after the specimen failed. As for specimen l-2-3, the number of cracks near the outer tips of fissure A and fissure B increased sharply. The farfield cracks around the specimen gradually coalesced. Little particle contact occurred near the middle of fissure A, and the particle contact force appeared with stress concentrations along the direction of macroscopic crack coalescence. The far-field cracks around the specimen and the outer tip wing cracks of fissure A and fissure B fully coalesced, and the cracks in the bridge area also fully coalesced. Particle contact force concentration only occurred in the bridge area, fissure tips, and other small local areas. Figures 18-20 show the crack coalescence processes and particle contact forces of specimens containing one fissure and a single circular hole under coupled biaxial loading and unloading. At Stage I, the wing cracks were initiated from the right tip of fissure B of specimen c-2-1 and expanded in two directions: vertically upward and downward right. Some cracks parallel to fissure B were initiated from the left tip of the fissure. No obvious crack propagation behavior occurred near the circular hole. The particle contact force near the left and right sides of the fissure and circular hole was larger. Tensile cracks perpendicular to the  Figure 19: The crack coalescence process and particle contact force of specimen c-2-2 under coupled biaxial loading and unloading.
17 Lithosphere fissure were initiated near the tips in specimen c-2-2, and some cracks that were initiated from the left tip of the fissure propagated in a direction approximately parallel to the fissure. The tensile cracks on the right side of the circular hole expanded in a vertical upward direction, while the cracks on the left side of the circular hole expanded in a vertical downward direction. The particle contact force between the circular hole and the fissure became larger. For specimen c-2-3, cracks were initiated from the left sides of the hole and fissure, and then cracks were initiated from the right sides of the hole A and fissure B, in which the crack propagation direction was the same as the loading direction of axial stress. The particle contact forces at the hole and fissure top and bottom sides were weaker.
At Stage II when the axial stress reached peak strength, cracks on the top left side of the circular hole in specimen c-2-1 gradually converged, and far-field cracks were initiated around the specimen. Cracks on the right side of the circular hole in specimen c-2-2 converged with the left tip of the fissure, and far-field cracks were primarily initiated in the bottom left and top right corners of the specimen. For specimen c-2-3, the tensile cracks on the right side of the hole continued extending downward along the axial loading direction, the wing crack near the left tip of the fissure rapidly expanded to the bridge area, and the particle contact force in the rock area was more concentrated. Overall, the cracks connected with the defects of hole A and fissure B at peak strength stage of the rock specimen.
At Stage III after the axial stress was loaded to peak strength, cracks at the two tips each extended in opposite directions but eventually joined each other and converged with the far-field cracks for specimen c-2-1. Similarly, the cracks on the left end of the circular hole converged with the far-field cracks on the left end of the specimen, but no significant crack convergence occurred on the right end. The particle contact force was significantly greater at the points where cracks converged. For specimen c-2-2, a large number of cracks were generated on the right tip of the fissure after peak strength and coalesced with the far-field cracks on the right side of the specimen. The cracks at both fissure tips were interconnected, and the cracks on the left side of the circular hole extended in a vertical downward direction but did not connect directly with the far-field cracks on the left side of the specimen. Eventually, a large area of particle contact stress concentration formed on the left side of the circular hole. For specimen c-2-3, the wing cracks near the left side of the fissure coalesced with cracks near the right side of the hole, and the wing cracks near  18 Lithosphere the right side of the fissure continued extending downward. The particle contact force in the area above the fissure became obviously smaller, the blank area without contact force became larger, and the particle contact force stress concentration area in the bridge area became narrower. Then, the far-field cracks around the specimen gradually coalesced and connected with the outer sides of the hole and fissure. The tensile cracks on the right side of the hole continued expanding toward the specimen bottom, and particle contact force concentration only occurred in the bridge area and other small local areas. Figures 21 and 22 show the crack coalescence processes and particle contact forces of specimens containing nine fissures and specimens containing three circular holes and six fissures under coupled biaxial loading and unloading. At Stage I before peak strength, the cracks were initiated and gradually extended from the fissure tips and diagonal boundaries of holes and tended to connect between two adjacent rows up and down, and the particle contact forces near the areas above were greater. At Stage II at peak strength, some of the three rows of defects would connect, and the particle contact forces between the defects grew. At Stage III after peak strength, a large number of cracks were generated between defects and coalesced with the far-field cracks of the specimen.

Conclusions
In this study, uniaxial compression, biaxial compression, and coupled biaxial loading and unloading were numerically simulated for granite specimens containing different shapes and numbers of prefabricated defects. The main conclusions are summarized below: (a) For the same specimen, the peak strains for biaxial compression, biaxial unloading, and uniaxial compression are in a decreasing order. The uniaxial compressive strength, biaxial compressive strength, and biaxial unloading compressive strength decreased as the number of prefabricated defects increased. The overall results of the specimens containing two fissures and specimens containing one fissure and a circular hole show that compressive strength tends to decrease as the defect center moves further away from the central axis of the specimen, with biaxial compressive strength, biaxial unloading strength, and uniaxial compressive strength decreasing in order. The biaxial unloading of the initial crack stresses for specimens at the same number was greater than the uniaxial compression of initial crack stress, and they both decreased as the number of  Figure 21: The crack coalescence process and particle contact force of specimen l-9 under coupled biaxial loading and unloading.
19 Lithosphere prefabricated defects increased. As the defect center moved from the central axis of the specimen, the biaxial unloading initial crack stress weakened. The peak biaxial unloading lateral stress increased and then decreased as the defect centers moved further away from the central axis of the specimen (b) The relationship between the numerical axial stresstime curve, the crack distribution, and cumulative crack numbers under coupled biaxial loading and unloading were analyzed. The results show that the cumulative number of cracks grew slowly before unloading and increased faster from when unloading began to before peak strength, while the crack growth rate after peak strength was significantly higher than that before peak strength. The specimens containing nine prefabricated defects had a shorter rapid rise in cumulative crack numbers. For specimens with two prefabricated defects, as the prefabricated defect center move father from the center of the specimen, the ultimate cumulative number of total cracks, shear cracks, and tensile cracks in the specimen lessened, and the first crack was initiated earlier. The cumulative number of cracks in specimens containing a single defect was significantly more than that of other specimens (c) A detailed summary of crack initiation, crack propagation, and crack coalescence of granite specimens under coupled biaxial loading and unloading conditions was conducted with numerical simulation. Different shapes, number of defects, and relative defect positions large affected crack initiation, crack propagation, and crack coalescence. Generally, the cracks were initiated and gradually extended from defect boundaries and tended to connect between two adjacent defects, and the particle contact force nearby the defect boundaries was larger before peak strength. At peak strength, some defects would connect through cracks extending between the bridge areas, and the particle contact force between the defects grew. After peak strength, a large number of cracks were generated between defects and coalesced with the far-field cracks of the specimen

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare no conflicts of interest.  Figure 22: The crack coalescence process and particle contact force of specimen c-9 under coupled biaxial loading and unloading. 20 Lithosphere