Unloading-Induced Crack Propagation of Two Collinear Unequal Length Flaws in Brittle Rocks

The propagation and coalescence of numerous discontinuous joints significantly contribute to landslide instability during excavation unloading. The tip expression of stress intensity factors of two collinear unequal length cracks in a typical rock mass under unloading conditions was calculated based on the superposition principle and fracture mechanics to determine the mesoinfluence law of intermittent joint interaction in the slope under the action of excavation. The effects of many factors on this interaction were also analyzed theoretically. Unloading tests were conducted on rock-like specimens with two collinear unequal length cracks in addition to numerical simulation and theoretical analysis. The results show decreased interaction between the two cracks with increased crack distance, increased influence of the main crack on a secondary crack with increased length of the main crack, and decreased influence of the secondary crack on the main crack with decreased length of the secondary crack. Wing tensile cracks first appear at the tip of flaws, and the propagation of these cracks occurs with the generation of secondary tensile cracks and shear cracks during unloading. Propagation and coalescence between cracks lead to tension and shear mixed failure of a rock bridge, and tensile cracks appear near the unloading surface. The axial initiation and peak stress of a crack increase with increased flaw distance, and the theoretical calculations were confirmed by lateral unloading test results.


Introduction
Numerous discontinuity cracks in rock slopes control the stability and failure mode of a rock slope. At the end of a discontinuity plane in a rock mass, high-stress concentration occurs during of rock slope excavation and unloading. This stress leads to crack propagation and coalescence and can further evolve into sudden local or global instability failure of the macroscopic upper slope, presenting a significant threat to the surrounding environment and operators [1,2]. Therefore, determining the crack propagation mechanism and interaction law of a rock slope under excavation unloading has important guiding significance for slope excavation design, construction, and long-term safe operation.
There have been several attempts to study the mechanism of propagation and evolution and the law of interaction between multiple cracks in the rock mass. Many effective methods have been proposed, but each has limitations. For multicrack interaction analysis and solution, the most widely used method is the singular integral equation method [3], however, when the crack distance is changed arbitrarily, the calculation will become quite complex. Huang and Karihaloo [4] simplified the singular integral equation method to obtain an approximate solution of crack interaction. Horii and Nemat-Nasser [5] analyzed the interactions between multiple cracks in infinite plates based on a pseudo-traction method, but this approach cannot be used to analyze closely spaced cracks. Kachanov [6] further analyzed the interactions of multiple cracks and proposed the expression of stress intensity factors at the tips of cracks based on a simplified pseudo-traction method, however, for a large crack distance, this method gives a significant calculation error. Li et al. [7] improved the accuracy of the calculation method of Kachanov to analyze of the interactions of closely spaced cracks and then applied the improved method to the analysis of compression-shear cracks. Kuang and Chen [8] used an alternating iterative method to deduce the stress intensity factors at the tips of cracks. Qing and Yang [9] proposed a new method to study multicrack interaction in infinite plates that incorporated aspects of both the alternating iteration method and the classical Kachanov method.
In addition, several studies have examined the propagation and evolution mechanism and interactions between multiple flaws in rock mass using physical tests and numerical simulation approaches. Zhang et al. [10][11][12][13][14] investigated the crack propagation and coalescence of nonisometric flaws in brittle rocks by digital imaging and AE techniques under uniaxial loading. Yang et al. [15][16][17] investigated the crack coalescence behaviour of sandstone samples containing two coplanar and two nonparallel flaws during deformation. Similarly, Cheng et al. [18] studied the crack propagation modes and failure modes of rock specimens with three flaws during deformation. Tang et al. [19,20] developed a numerical method, Realistic Failure Process Analysis (RFPA), to simulate the progressive failure of rocks. Zhou et al. [21] proposed a new meshless numerical simulation method that incorporated generalized particle dynamics (GDP) and used this approach to simulate the crack initiation, propagation, and coalescence in rock specimens with prefabricated flaws under uniaxial compression. Wong and Einstein [22] used a grain flow program (PFC) to simulate the propagation and coalescence modes of prefabricated flaws in rock specimens under uniaxial compression. Chen et al., Wang et al., and Zhao et al. [23][24][25] studied the crack propagation of rocklike specimens with one or two flaws by experimental and numerical methods under conditions of uniaxial and biaxial compressions. Similar studies include Zhou et al. [26], Wang et al. [27], Lin et al. [28], Zhao et al. [29], and Wong and Xiong [30].
The above studies on the propagation, evolution, and interactions cracks in rock a mass mainly focused on loading mechanics, and only few studies examined the interactions of rock mass cracks under unloading conditions. Zhao et al. [31] studied the influence of the unloading rate on the strain burst characteristics of Beishan granite under true triaxial unloading conditions. Wang et al. and Zheng et al. [32,33] addressed the influence of stress-strain, strength, and deformation characteristics, as well as the mechanical properties of artificial intermittent joint rock samples by conventional triaxial unloading and reloading of rock mass after unloading damage. Chen et al. [34] studied the expansion model of a jointed rock bridge under the unloading condition using a true triaxial test system.
Most previous investigations of multiple cracks in a rock mass under unloading conditions concentrated on qualitative analysis of the influence of different geometric distributions of cracks on the mechanical properties of a rock mass using experimental tests and numerical simulation, but there remains insufficient quantitative theoretical research. This problem is complicated because the stress state of the rock mass is complex under the unloading condition, and it is difficult to directly obtain the stress intensity factor at the tip of the crack. In this study, the stress intensity factors at the tips of two collinear cracks of unequal length were determined under unloading conditions based on the superposition principle and the theory of crack mechanics. The contributions of different geometric layouts and the sizes of double cracks on the process of propagation and the law of interaction were assessed based on the ratio of the stress intensity factor. Finally, the theoretical calculation was verified using unloading experiments of rock-like materials and RFPA numerical simulation. The results provide an important reference for the study of rock mass instability mechanism induced by excavation and suggest strategies for slope design and excavation scheme. In mine rock slope engineering, excavation unloading is a main contributor to slope instability. After a slope is excavated, the in situ stress in the horizontal direction is suddenly unloaded to zero in a certain direction, and this change in stress state will lead to differential spring-back deformation of the rock mass perpendicular to a certain depth of the excavation surface. This deformation will gradually decrease from the excavation surface to the slope rock mass, and differential spring-back deformation will cause tensile stress around the crack surface, which will change the stress state of the rock mass from a compressive shear stress state to a tensile shear stress state [35]. According to Huang [36], stress field redistribution occurs in the rock mass of the slope after excavation. The secondary stress field can be divided into three regions: the area of reduced stress (failure area), the area of increased stress, and the original rock stress area, as shown in Figure 1. Local or complete slope landslide is often caused by the propagation and coalescence of intermittent cracks in the stress reduction area. Here, we selected two intermittent collinear cracks of unequal length in the stress reduction area as the focus of the study and then discussed the stress state of the cracks and the interaction between the cracks. As described previously, before excavation, a rock mass in the stress reduction area is subjected to vertical in situ stress σ 1 and horizontal in situ stress σ 3 , as shown in Figure 1(a). After excavation, the stress path has changed. The vertical direction is still affected by σ 1 , but in the stress reduction area, the horizontal direction changes from compressive stress to tensile stress, and the stress state of the rock mass changes from a compressive shear stress state to a tensile shear stress state, as shown in Figure 1(b).

Solution of Stress Intensity Factors at the Tips of Single
Crack. The stress state of a cracked rock mass in the slope changes from a state of compression-shear stress to a state of tension shear-stress during excavation and unloading. It is difficult to directly obtain the stress intensity factor at the tip of a crack based on two-dimensional fracture mechanics. According to the superposition principle, for the on-line elastic range, the stress intensity factor at the tip of a crack under two or more different loads is equal to the sum of the stress intensity factor under each load alone under the same boundary conditions [37][38][39][40]. Therefore, we generalized the model of an unloading rock mass in slope excavation and performed analysis according to two-dimensional fracture mechanics.

Geofluids
This was done from the point of view of the plane strain, and superposition was applied to analyze the stress state of the cracked rock mass under the unloading condition [41]. The stress state for a cracked rock mass under unloading condition was decomposed into the uniaxial compressive stress state B and the uniaxial tensile stress state C. The principle of decomposition is shown in Figure 2, where α is the crack dip angle in rock mass under uniaxial compression, β is the crack dip angle in the uniaxial tensile state, and α + β = 90°. According to the decomposition principle, the stress intensity factor of the crack tip in stress state A can be calculated as follows: where K A is the stress intensity factor at the tip of the crack under unloading condition, K B is the stress intensity factor at the tip of the crack under uniaxial compression, and K C is the stress intensity factor at the tip of the crack under uniaxial tension.

Stress Intensity Factors at the Tips of Crack under
Uniaxial Compression. The crack will close under uniaxial compression stress, and the rock mass will experience shear sliding failure along the crack surface. The normal stress on closed cracks can be neglected, and the crack tip only has K B II . The stress intensity factor at the tip of crack can be expressed as [42]: where λ is the half-length of crack, μ is the coefficient of friction in the crack plane, and K B II is the mode II stress intensity factor at the tip of the crack.

Stress Intensity Factors at the Tips of Crack under
Uniaxial Tension. The stress intensity factor at the tip of a crack under uniaxial tension can be expressed as [36]: where K Ι and K II are the mode I and II stress intensity factors at the tip of a crack under uniaxial tension, respectively.

Stress Intensity Factors at the Tips of Crack under
Unloading Conditions. Based on the above results, the stress intensity factor at the tip of a crack under unloading condition described in (2) and (3) and substituted into (1) can be     3 Geofluids obtained as follows: where K A Ι and K A II are the mode I and II stress intensity factors at the tip of cracks under unloading condition, respectively.

Solution of Stress Intensity Factors at the Tips of Two
Collinear Cracks of Unequal Length. The process to solve for the stress intensity factors at the tips of two collinear cracks is the same as that for a single crack under unloading condition. The stress state of double cracked rock mass under unloading condition can be decomposed into the uniaxial compressive stress state M and the uniaxial tensile stress state N, with α + β = 90°. As shown in Figure 3, the stress intensity factors at the tips of crack unloading condition can be obtained as follows: where K D is the stress intensity factor at the tip of the crack under unloading condition, K M is the stress intensity factor at the tip of the crack under uniaxial compression, and K N is the stress intensity factor at the tip of the crack under uniaxial tension.

Stress Intensity Factors at the Tips of Crack under
Uniaxial Tension. The mechanical model of a cracked rock mass under uniaxial tension is shown in Figure 4. There are two collinear, unequal cracks in an infinite plane, which are subjected to tensile stress σ 3 at infinity, and the angle between the cracks and the horizontal direction is β. The coordinate system is established as shown in Figure 4, and the midpoint of the tip line in the two cracks is the origin of the coordinates. The distance between the cracks is 2d, the length of the main crack is 2l 1 ða 2 − dÞ, and the length of the secondary crack is 2l 2 ða 1 − dÞ.
According to the definition of stress intensity factor, composite failure of mode I-II will occur on the two collinear crack tips of unequal length under uniaxial tension, and the stress intensity factors at the inner and outer tips of crack can be expressed as follows [43]:   Geofluids Þare the mode I stress intensity factors at the inner and outer tips of the main crack and secondary crack, respectively, under uniaxial tension, and K N II ðdÞ, K N II ða 2 Þ, K N II ð−dÞ, and K N II ð−a 1 Þ are the mode II stress intensity factors at the inner and outer tips of the main crack and secondary crack, respectively, under uniaxial tension.The D, E, and F can be expressed as: The expression of I 2 ðn, kÞ is: where F(k) and Πðn, kÞ are the complete elliptic integrals of the first kind and the third kind, respectively, and the simplified expressions can be written as: The module k of the first kind of the complete elliptic integral F (k) and the third kind of complete elliptic integral can be determined by the following formula.

Stress Intensity Factors at the Tips of Crack under
Uniaxial Compression. Under uniaxial compression, the crack will close under normal compressive stress, so there will only be mode II stress intensity factors at the tips of the crack. Considering the model of compression-shear collinear cracks shown in Figure 5, the stress intensity factors at the inner and outer tips of the crack can be expressed as follows: where K M II ðdÞ, K M II ða 2 Þ, K M II ð−dÞ, and K M II ð−a 1 Þ are mode II stress intensity factors at the inner and outer tips of the main and the secondary crack, respectively, under uniaxial compression, τ c is the cohesion on the crack surface, and μ is the friction coefficient on the crack surface.

Stress Intensity Factors at the Tips of Crack under
Unloading Conditions. According to the previous superposition principle, the stress intensity factors at the tips of the crack under unloading condition can be obtained as follows: where K D Ι ðdÞ, K D Ι ða 2 Þ, K D Ι ð−dÞ, and K D Ι ð−a 1 Þ are the mode I stress intensity factors at the inner and outer tips of the main crack and the secondary crack, respectively, under unloading condition, and K D II ðdÞ, K D II ða 2 Þ, K D II ð−dÞ, and K D II ð−a 1 Þ are the mode II stress intensity factors at the inner and outer tips of the main crack and the secondary crack, respectively, under unloading condition. The explicit expression of stress intensity factors at the tips of crack K D under unloading conditions can be obtained from the expressions of stress intensity factors at the tips of the crack under uniaxial tension ((6)~(9)) or under uniaxial compression ((15)~(16)).

Interaction Analysis and Discussion of Two Collinear Cracks of Unequal Length
The interaction between cracks can increase, decrease, or have no effect on the stress intensity factors at the tips of cracks. In this study, we defined the ratio of stress intensity factors as = F = K D /K A , which represents the ratio of the 5 Geofluids stress intensity factors at the tips of a crack when there are two cracks to the stress intensity factors when there is only a single crack. We then used the F of size to characterize the degree of influence between the two cracks. According to the literature [44], crack propagation is only related to the specific parameters of a particular crack. For main and secondary cracks of equal length, the cracks have virtually no effect on each other, so F = 1:05. For a ratio of stress intensity factor at the tips of crack of F ≤ 1:05, there is no effect of interaction on crack propagation. Determining the stress intensity factors at the tips of the crack revealed equal stress intensity factors for mode I and mode II at the inner and outer tips of the crack, therefore, we only analyzed the variation of the stress intensity factor ratio of mode I with the distance and size of cracks.
3.1. The Influence of Crack Distance. The influence of crack distance on the interaction of cracks can be analyzed by keeping the length of cracks constant, that is, l 1 = 2 l 2 = 0:8, and changing the crack distance. Table 1 shows the calculated results of the mode I stress intensity factor ratio at the crack tips for different crack distance values. As shown in Table 1, compared with a single crack, there are larger stress intensity factors at the tips of two cracks. In addition, the F decreases with the increase of crack distance, indicating decreased interaction between cracks with increased crack distance. The variation of F of the outer tip of the main crack is smallest with the variation of crack distance, which indicates little impact on the outer tip of the main crack by the interaction between cracks. The variation of the inner tip F of the secondary crack is the largest with the variation of crack distance, which indicates that the inner tip of the secondary crack is most impacted by the interaction between cracks. When the crack distance is close to or exceeds the length of the secondary crack, that is, 2d ≥ 0:8, the main crack approaches F ≤ 1:05, which indicates that the main crack propagation is little impacted by the secondary crack. When the crack distance is close to or larger than the length of the main crack, that is, 2d ≥ 1:6, the secondary crack approaches F ≤ 1:05, indicating that secondary crack propagation is little impacted by the main crack.
3.2. The Influence of the Main Crack. The influence of the length of the main crack on crack interaction can be analyzed by keeping the length of the secondary crack and crack distance constant, that is, d = 2 l 2 = 0:2, and changing the length of the main crack. Table 2 shows the calculated results of the mode I stress intensity factor ratio at the crack tips for different lengths of the main crack. As shown in Table 2, when the crack distance is equal to the length of the secondary crack, the F of the secondary crack increases almost linearly with the increase of the length of the main crack. The inner tip F of the secondary crack is slightly larger than that of the outer tip, but the F variation of the main crack is small, which indicates that the influence of the main crack on the secondary crack increases linearly with the increase of the length of the main crack. When the length of the main crack exceeds the crack distance, that is, 2l 2 ≥ 0:4, the tip F in the main crack becomes larger than 1.05, which indicates that the main crack begins to affect the inner tip propagation of the secondary crack. When the length of the main crack is equal to 1.5 times the length of the crack distance, that is, 2l 2 ≥ 0:6, the outer tip F of the secondary crack becomes larger than 1.05, which indicates that the main crack begins to affect the outer tip propagation of the secondary crack. There is less variation of F of the main crack, which indicates that the main crack is less affected by the secondary crack.  Table 2: Relationship between the ratio F of stress intensity factor and the length of the main crack (l 1 ).
The length of main crack l 1 (mm) The ratio F of stress intensity factor-K D 3.3. The Influence of the Secondary Crack. The influence of the length of the secondary crack on crack interaction was analyzed by keeping the length of the main crack and crack distance constant, that is, l 2 = 2d = 0:3, and changing the length of the secondary crack. Table 3 shows the calculated results of the mode I stress intensity factor ratio at the crack tips for different lengths of the secondary crack. As shown in Table 2, when the length of the secondary crack is equal to the distance between cracks, the F of the main crack decreases with the decrease of the length of the secondary crack and the decrease of the inner tip is larger than that of the outer tip. This indicates a decreased effect on the propagation of the main crack by the secondary crack with decreased length of the secondary crack. However, there is a greater influence of the inner tip of the main crack compared to the outer tip. When the length of the secondary crack is equal to the length of the main crack, the F values at the inner and outer tips of the two cracks are equal, indicating the same influence between the cracks. When the length of the secondary crack is equal to the crack distance, that is, 2l 2 ≤ 0:3, the F of the main crack is less than 1.05, which indicates the propagation of the main crack is not impacted by the secondary crack.

Testing System.
Based on the slope shown in Figure 1, the unloading test was based on the stress path of the slope before and after excavation, from the point of view of the plane strain (independent of intermediate principal stress). We simulated the process of slope excavation and unloading using a biaxial loading and unloading test and analyzed the propagation and evolution of nonpenetrating collinear intermittent cracks of unequal length in the rock specimen during unloading. The test was carried out using the MS-500 triaxial impact rockburst test system of the State Key Laboratory of Deep Geotechnical Mechanics and Underground Engineering of China, University of Mining and Technology (Beijing). This system can realize the loading and unloading from three sides independently, and the three-dimensional stress data can be recorded in real-time [45]. The test process system is shown in Figure 6, and a SA-5 high-speed camera is used to collect images of the specimen during the process of failure.

Specimen Preparation.
Because it is difficult and expensive to generate a slit in rock specimens, rock-like materials were used to fabricate flawed specimens to verify the accuracy of the above theoretical analysis. The rock-like material used here was composed of 425 ordinary Portland cement, standard sand, and water, with a mass ratio of 1 to 2.35 to 0.5. According to the relevant standards of the International Society of Rock Mechanics, the specimen size was machined into 110 mm × 110 mm × 30 mm samples, as shown in Figure 7. The flaws were made of high strength thin steel sheet with a thickness of 0.5 mm. After mixing the material evenly at the appropriate proportion, the material was poured into a mold to prepare a test block with the predetermined mix ratio. The steel sheet was extracted before the initial solidification of the material, and the block was maintained at room temperature for 12 hours to demould, and then maintained in a curing room for 28 days. After Table 3: Relationship between the ratio F of stress intensity factor and the length of secondary crack l 1 .
The length of secondary crack l 2 (mm) The ratio F of stress intensity factor-K The location of specimen The system of image acquisition Triaxial testing machine The control system of test Day lighting lamp Figure 6: The true-triaxial unloading test system. 7 Geofluids removal from the mold, the flatness and penetration of the flaws were assessed and then subjected to grinding for the preparation of the specimen. The mechanical parameters of the resulting specimens are shown in Table 4. As previously described [46], the mechanical parameters of the standard specimens used here are similar to those of sandstone. To eliminate any boundary effect to the maximum extent, the prefabricated flaw was generated in the middle of the specimen. The final test scheme and crack size are shown in Table 5.

Testing Procedures.
To effectively simulate the process of slope excavation, biaxial unloading was used. Considering that the test is a plane strain problem, the stress direction was set as shown in Figure 8(a), with the specific stress path as shown in Figure 8(b). Each group of tests was repeated three times independently to reduce the discreteness of the test results. The specific operation steps of the test are as follows: (1) First, the axial stress σ 1 and the horizontal stress σ 3 were simultaneously loaded to the specimen at 0.004 mm/s by displacement control, and the specimen was then fixed (2) The axial stress σ 1 and horizontal stress σ 3 were simultaneously loaded to σ 1 = σ 3 = 15 MPa at 0.004 mm/s by displacement control, and the horizontal displacement was then locked (3) With constant horizontal displacement, the axial displacement was loaded to σ 1 = 30 MPa at 0.004 mm/s by displacement control, and then the axial displacement was locked (4) The stress control mode was used to quickly unload σ 3 to 0, and the axial displacement continued at the rate of 0.004 mm/s until the specimen was destroyed and the test was complete

Analysis of Stress-Strain Curve and Mechanical
Parameters. As shown in Figure 9, the specimen of rocklike material showed an obvious brittle stress drop. The axial stress-strain curve was not smooth, and many jumps were observed due to crack propagation in the specimen. The curve compaction stage of the specimen was long, with an obviously shortened linear elastic stage, and the multijump phenomenon of the curve corresponds well to the crack propagation. As shown in Figure 10, the average axial peak stresses of specimens 1 to 4 were determined as 55.   Table 5: Flaw geometric parameters in specimens.

Specimen number
The 8 Geofluids flaw initiation stress and peak stress with increased flaw distance. This is mainly due to decreased interaction between flaws with increased distance between flaws, which makes the initiation of the flaw tip more difficult. Figure 11 shows the process of crack propagation for specimen 1. During unloading, wing tensile cracks are produced at the inner and outer tips of the secondary flaw and the inner tip of the main flaw.

The Analysis of Crack Propagation.
The outer tip of the main flaw exhibited delayed initiation. With continuous unloading, wing tensile cracks propagated along the direction of the maximum principal stress, shear cracks occurred in the intermediate rock bridge area, and the propagation of wing tensile cracks and shear cracks resulted in tension and shear coalescence of the rock bridge. Further, secondary tensile cracks occur with the propagation of wing tensile cracks at the tip of the secondary flaw. Finally, the cracks coalesce with each other to form a local detachment zone, resulting in failure of the specimen. Secondary tensile cracks appeared near the unloading surface. Figure 12 shows the process of crack propagation in specimen 2 during unloading. Wing tensile cracks appeared at the tips of the main and secondary cracks. Shear crack was produced in the middle rock bridge area, and the secondary flaw and the outer tip of the main flaw delayed the initiation of the crack. The wing tensile cracks expanded along the direction of the maximum principal stress propagation, and secondary tensile cracks appeared in the center near the unloading surface. Finally, the shear cracks at the outer tip of the secondary flaw coalesce with each other and with the tensile crack to form a local detachment zone, resulting in failure of the specimen. Figure 13 shows the process of crack propagation for specimen 3 during unloading. The inner and outer tip of the main flaw and the outer tip of the secondary flaw formed wing tensile cracks along the direction of the maximum principal stress propagation, and the cracks at the two tips of the main flaw also exhibited upward propagation. With continuous unloading, the outer tips of the two flaws continued in the direction of the maximum principal stress propagation, with no coalescence in the middle rock bridge area. Secondary tensile cracks occurred with the propagation of wing tensile cracks at the tip of the main flaw. Finally, the cracks coalesce to form a local detachment zone, resulting in failure of the specimen. Figure 14 shows the process of crack propagation of specimen 4 during unloading. The inner and outer tips of the secondary flaw germinated wing tensile cracks along the direction of the maximum principal stress propagation. With continuous unloading, secondary tensile cracks occurred with the propagation of wing tensile cracks at the outer tip of the main flaw. Finally, the tensile cracks coalesce with each other to form a local detachment zone, resulting in failure of the specimen. Secondary tensile cracks appeared near the unloading surface.
(1)The Influence of Flaw Distance. As shown in Figure 11, when the flaw distance is smaller than the length of the secondary flaw, there is a greater influence of the main flaw on the tip of the secondary flaw, especially the inner tip. This results in the preferential initiation of the inner and outer tips of the secondary flaw during unloading. However, the secondary flaw has only a weak influence on the main flaw and has the weakest effect on the outer tip of the main flaw, resulting in the initiation of the inner tip of the main flaw before the outer tip during unloading. As shown in Figure 12, the influence of the main flaw on the secondary flaw gradually weakens with the increase of flaw distance,  9 Geofluids which delays the initiation of the outer tip of the secondary flaw. Similarly, the influence of the secondary flaw on the main flaw also weakens with the increase of flaw distance, which then delays the initiation of the outer tip of the main flaw. As shown in Figure 13, the flaws do not influence each other with the continuous increase of flaw distance. The above test results are consistent with the results presented in Section 3.1.
(2)The Influence of the Main Flaw. As shown in Figure 13, when the length of the main flaw is equal to the flaw distance, the flaws have no influence on each other. As shown in Figure 12, when the length of the main flaw is equal to twice the crack distance, the main flaw has a significant influence on the secondary flaw, which results in the initiation of the inner tip of the secondary flaw during unloading. However, the secondary flaw has little influence on the main flaw, which delays the initiation of the outer tip of the main flaw. As shown in Figure 14, when the length of the main flaw is equal to three times the flaw distance, the inner and outer tips of the secondary crack initiate first during unloading, and the main flaw is less affected by the secondary flaw. The above test results are consistent with the results presented in Section 3.2.
(3)The Influence of the Secondary Flaw. As shown in Figure 11, when the length of the secondary flaw is larger than the flaw distance, the influence between flaws is obvious, and the inner and outer tip of the secondary flaw and the inner tip of the main flaw initiate first during unloading. However, the influence on the outer tip of the main crack is the smallest. As shown in Figure 12, when the length of the secondary flaw is equal to the flaw distance, the secondary flaw has little influence on the main flaw, which delays the initiation of the outer tip of 10 Geofluids the main flaw. As shown in Figure 13, when the length of the secondary flaw is less than the flaw distance, the flaws have no effect on each other. The above test results are consistent with the analysis results presented in Section 3.3.
According to the failure process of specimens 1 to 4 shown in Figures 11-14, it can be seen that the flaws control the failure of specimens under unloading conditions. Wing cracks occur at the flaw tips during unloading, and the inter-action of the two flaws result in tension-shear mixing coalescence of the rock bridge. Secondary shear cracks occur with the propagation of wing cracks and tension cracks form near the unloading surface under further unloading. The failure of the specimens is caused by the coalescence of the cracks. The main flaw contributes more to the failure process than the secondary flaw. The main and secondary flaws propagate and coalesce with the tensile cracks near the unloading surface, resulting in specimen failure.

Principle of Numerical Analysis.
To verify the theoretical and experimental results, rock crack failure analysis software, RFPA 2D , was applied. RFPA 2D was developed by the Research Center of Rock Crack and Instability of Northeast University and can simulate the damage and failure process of materials. In RFPA 2D , FEM is used to calculate and analyze the stress field, and heterogeneity is considered by assigning different material properties to the mesoscopic elements of a solid or structural model. Based on elastic damage mechanics, the constitutive relationships of these mesoscale elements in tensile and shear modes are investigated to determine whether these elements reach their damage thresholds as judged by the maximum tensile criterion or the Mohr-Coulomb criterion [19,20].

Heterogeneity Consideration.
Rock can be a heterogeneous material with spatial variation in mechanical properties. To consider the heterogeneity of a rock mass, the RFPA 2D analysis program was used to discretize the numerical model into many mesoscale elements. The mechanical properties and physical parameters of each mesoscale element, including elastic modulus, Poisson's ratio, strength, and density, are discrete and assumed to conform to the Weibull distribution. As a result, the relationship between mesoscopic and macroscopic mechanical properties can be established. Here, the probability density function φðαÞ is given by: where α is a material (rock) mechanical property parameter (such as the elastic modulus, strength, or Poisson's ratio); α 0 is a scale parameter denoting the average value of a material property; and m is a shape parameter that defines the shape of the distribution function. The physical meaning of m reflects the homogeneity of the material (rock) medium. Figure 15 shows the Weibull distribution for material properties with different homogeneity indices m. When the homogeneity index m has a relatively high value, the mechanical properties and physical parameters of mesoscale elements in the rock are in a very small range, indicating that the rock is relatively homogenous. In contrast, when the rock's properties are widely distributed with a relatively lower m index, there is greater inhomogeneity of the rock's mesoscale element properties. The homogeneity index is considered an effective approach to assign material properties for inhomogeneous rock and concrete [46][47][48].    [49], considering the process of material damage, and established a set of damage models: where σ is Cauchy stress; σ e is effective stress; E is Young's moduli; and D is the damage variable. When D > 1, the material is not damaged and this is an ideal state. When D = 1, the material is completely damaged. When 0 < D < 1, the material shows some degree of damage.In the calculation model, the inhomogeneous material is composed of many mesoscale elements with different parameters. After investigating the stress level of these elements in the uniaxial stress state, the following four basic states are determined [50,51], as shown in Figure 16: (1) Elastic state. When σ 1 − σ 3 ð1 + sin ϕ/1 − sin ϕÞ < f c or σ 3 > −f t , the element is in an elastic state. Where f c is compressive strength, f t is tensile strength, and ϕ is the friction angle (2) Damage state. When σ 1 − σ 3 ð1 + sin ϕ/1 − sin ϕÞ ≥ f c , because the element damage is induced by compressive stress, the elastic modulus is weakened according to the elastic-brittle damage constitutive equation with residual strength The evolution of damage variables is defined as follows: where f c is residual strength.
When σ 3 ≤ −f t , because element damage is induced by compressive stress, the elastic modulus is weakened according to the elastic-brittle damage constitutive equation with residual strength. The evolution of damage variables is defined as follows: where f c is the tensile residual strength (3) Fracture state. When the tensile strain reaches the ultimate tensile strain ðε − ε tu Þ, the element completely loses its bearing capacity and stiffness, and its elastic modulus is small, resulting in a crack element (4) Crack closure state. When the compressive strain of the crack element reaches the ultimate compressive strain, the crack is in a closed state, and the elastic modulus of the element increases with the increase of compressive stress  Table 6. The stress path of the numerical experiment is the same as that of the biaxial unloading physical experiment.

Numerical
Results. The initiation and propagation processes of cracks in Models 1-4 under unloading are shown in Figure 14. As shown in Figures 17-20, when the distance between flaws is smaller than the length of the secondary flaw (Figure 17), the main flaw has a significant impact on the secondary flaw. The influence of the main flaw on the secondary flaw gradually weakens with increased distance between flaws ( Figure 18), delaying the initiation of the outer tip of the secondary flaw. The flaws do not affect each other with the continuous increase of flaw distance ( Figure 19). When the length of the main flaw is equal to three times the distance between flaws (Figure 20), the inner and outer tips of the secondary crack initiate first unloading, with less effect on the main flaw of the secondary flaw.

Discussion
Analysis of crack propagation by unloading test and RFPA 2D numerical simulation revealed that the propagation law of flaws was similar in the unloading test to that modeled by numerical simulation. There was weakened interaction 13 Geofluids between flaws with the increase of crack distance. The influence of the main flaw on the secondary flaw increased linearly with the increase of the length of the main flaw.
As shown in Table 7, the mechanical parameters of the test and the numerical simulation are basically the same and the experimental path is the same, suggesting the accuracy of the numerical simulation. As shown in Figure 21, the numerical simulation results were very similar to the experimental results. The failure of the specimen exhibited symmetrical initiation and propagation along both tips of the flaws. The flaw tips first initiate wing tensile cracks during unloading, and the rock bridge is connected by tension-shear mixed coalescence. In model 3, the rock bridge does not exhibit coalescence because the length of the rock bridge is equal to the length of the main flaw and there is almost no effect between the two flaws, as seen in both the theoretical analysis and the test results. There are some differences between the simulation results and the test results. In the experiment, some secondary tensile cracks appear on both sides of the unloading surface, but this was not observed in the numerical simulation. Unloading is relatively stable during numerical simulation, but in physical tests, due to the oscillation of the oil pressure system, there may be some disturbance, which leads to the emergence of more cracks.  14 Geofluids Another explanation of the difference between simulation and experimental values is that the uniformity of the model can be guaranteed for the numerical simulation, but real specimens have heterogeneity that is observed in physical tests. The law of interaction between flaws in the numerical simulation is basically consistent with the results of physical tests and theoretical calculations.

Conclusions
In this study, the propagation and interaction laws of two unequal collinear cracks under unloading conditions were theoretically analyzed based on the superposition principle and fracture mechanics theory. Unloading test and RFPA 2D numerical simulation were performed to validate the theoretical analysis results. The conclusions can be summarized as follows: (1) Crack distance has an important influence on the crack interaction law. The influence between cracks decreases with increased crack distance. When the crack distance is equal to the length of the secondary crack, the main crack propagation is barely impacted by small crack. When the crack distance is equal to the length of the main crack, secondary crack propagation is little affected by the main crack  (2) Crack length also has an important influence on the crack interaction law. The influence of the main crack on the secondary crack increases linearly with increased main crack length. When the main crack length is larger than the crack distance, the main crack begins to affect the inner propagation of the secondary crack. With the decrease of the secondary crack length, the main crack propagation is impacted by the secondary crack. When the secondary crack length is equal to the crack distance, main crack propagation is not impacted by the secondary crack (3) Crack interaction was verified by unloading tests of rock-like material and RFPA 2D numerical simulation. Wing tensile cracks occurred at the tips of the flaw during unloading, and secondary tensile cracks and shear cracks occur with the propagation of the wing tensile cracks under further unloading. Finally, the cracks coalesce with each other to form a local detachment zone, resulting in specimen failure. Additionally, secondary tensile cracks form near the unloading surface (4) The flaw distance and size have significant effects on the crack initiation load and peak load of specimens. The axial crack initiation and peak load of the specimens increases with the increase of crack distance. In addition, the theoretical calculation results are consistent with the test and numerical simulation results, which verifies the rationality of the theoretical results

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 no conflict of interest.  16 Geofluids