An Experimental and Numerical Investigation on the Initiation and Interaction of Double Cracks in Rocks under Hydromechanical Coupling

The growth of double cracks is the main factor leading to progressive rock failure under hydromechanical coupling. The initiation modes and interaction behaviors of double cracks were investigated by using laboratory tests, and the inﬂuences of water pressure were analyzed. The maximum energy release rate criterion was modiﬁed to determine the crack growth characteristics. A numerical model was established and then veriﬁed by the test results. Based on the simulation, the distribution of stress ﬁelds and key fracture parameters of double cracks was investigated. Then, initiation characteristics and interaction behaviors of parallel and nonparallel cracks were quantitatively analyzed. The results indicate that the increase in water pressure leads to the crack initiation being inclined to the original surfaces and the growth length along the crack fronts tending to be uniform; the small tensile stress zones are formed close to the crack tips, and signiﬁcant compressive stress zones are formed at both sides of the crack surfaces; stress superposition and interaction occur when crack spacing is less than 2.5 a ; the interactive weakening eﬀect is mainly present in the inner side (rock bridge zone) of cracks, while a certain degree of interactive enhancement eﬀect exhibits in the outer sides; the cracks are much easier to initiate at the outer wing cracks when the spacing is less than the critical length (0.5 a ); and cracks with a dip angle of 45 ° are much easier to initiate at the endpoints of long axis. The research results provide certain theoretical guidance for the safety assessment of underground engineering.


Introduction
Discontinuities such as joints and cracks that are widely present in rock masses have a significant impact on the stability and safety of underground structures [1,2]. e unsteady growth and coalescence of joints and cracks are easily triggered due to the coupling effect of stress and seepage fields, especially in complex environments with high geostress and high groundwater pressure, thus leading to the rock failure and serious engineering accidents. erefore, the research on the growth and interaction mechanisms of cracks in rocks under hydromechanical coupling has received more and more attention in recent years.
To date, the problem concerning the growth of a single 3D crack has been studied extensively [3][4][5]. e previous work has laid a foundation for the prediction of the failure of rock containing a single crack [6][7][8]. However, studies and many engineering practices indicate that the initiation and growth modes of multiple cracks are quite different from that of the single crack due to the interaction between cracks, which is the main reason to cause large-scale rock instability [9,10]. It has been recognized that the spatial distribution of cracks (dip angle, spacing, and density) plays a vital role in the growth and interaction modes of double cracks [11]. Many experimental investigations and numerical simulations have been carried out and report that physicomechanical properties of rocks also have great influences on the crack initiation and coalescence behaviors [12][13][14]. Besides, the stress and seepage fields are significantly affected by the double cracks interaction behavior, which, in turn, changes the penetration of the rock bridge as well as the rock failure modes [15][16][17][18][19].
Currently, the growth and interaction of multiple cracks under the condition of pure ground stress have been well studied, while the problem considering the effect of hydromechanical coupling, which concerns the construction safety of deep underground engineering [20,21], still lacks an in-depth research. Notably, the growth and interaction mechanisms of multiple cracks with different water pressures and crack spatial distribution are not well understood, and the related quantitative analyses of rock failure under hydromechanical coupling are deficient. e laboratory tests were carried out using the rock-like specimen, each containing double 3D cracks, to investigate the initiation and interaction behavior of double cracks in rocks. Since the numerical method has a great advantage in dealing with complex 3D fracture problems of rocks, the numerical model was established based on the FEM (finite element method) code of FRANC3D V7.0 and verified by the test results. e maximum energy release rate criterion was modified by introducing the weight coefficient corresponding to three basic fracture modes. en, the distribution of stress fields and the key fracture parameters along crack fronts was quantitatively analyzed, and last, the initiation and interaction behavior of 3D cracks and the main influencing factors were clarified.

Specimen Preparation and Test Method.
e hydromechanical coupling test was carried out with specimens, each containing two cracks. e specimens were fabricated by using a transparent resin rock-like material, developed by Fu et al. [22], which had high transparency and is easy for the observation of cracking phenomena. Besides, the material exhibits good brittleness and similar mechanical properties to natural rocks at a low temperature.
Details of the specimen and the testing apparatus are shown in Figure 1(a). e external dimension of the rectangular specimen was L × W × H � 70 × 70 × 140 mm. Two elliptic hollow cracks (parallel or nonparallel) were fixed along the central axis of a specimen by a cotton thread, with a uniform long-axis length of 2a � 24 mm and a short-axis length of 2b � 12 mm. Each crack was created by sealing a pair of overlapped elliptic mica sheets with a thickness of 0.25 mm along their conjunct border. As shown in Figure 1(a), the asbestos washers were used to fill the conjunct border of the mica sheets to form hollow cracks. A small hole with a diameter of 2.0 mm was drilled on one side of the crack as a water injection hole. e crack dip angles (β 1 and β 2 ) were denoted as the angles between the crack surface and the horizontal plane, and the crack spacing s was defined by the vertical distance between two cracks from their centers.
eir values in the tests were β 1 , β 2 � 45°and s � 24 mm, respectively. Besides, two plastic hoses with an inner diameter of 1.5 mm, which connected to the upper and lower cracks through predrilled holes, respectively, were preembedded along the axis of the specimen as the water injection channel (see Figures 1(a) and 1(b)).
As mentioned in the work of Fu et al. and Li et al. [22,23], the resin specimens have properties similar to real rocks and high transparency at a temperature of −45°C.
erefore, the laboratory test was carried out after the resin specimens were placed in the freezer for more than 24 hours to get a stable temperature of −45°C.
e uniaxial compressive stress was applied to the specimen ends by using the rock mechanics testing system. Meanwhile, the water with different pressures (p w ) of 0 MPa (no water injection), 2.0 MPa, and 4.0 MPa, respectively, was injected into the hollow cracks from the preembedded hoses by the water injection system, as shown in Figure 1(c). Besides, the water was dyed red; therefore, during the test, the cracking behavior can be investigated by observing the diffusion of red water in the specimens. Figure 2 shows the initiation modes of cracks with different water pressures. Under the action of uniaxial compression and water pressure, the wrapping wing cracks were generated at the long-axis ends of prefabricated cracks, with large deflection angles to the original crack planes. e growth of the outer and inner wing cracks showed a certain degree of asymmetry due to the crack interaction effect, which is mainly reflected in the fact that the lengths of outer wing cracks were slightly longer than those of the inner ones. On the contrary, small lateral growth along the original crack planes was observed at the short-axis ends. However, the lateral growth was insufficient, especially in low water pressure or no water injection cases. e crack initiation angle is defined as the acute angle between the new and original crack surfaces. e spatial pattern of 3D crack growth is primarily controlled by the initiation angles and extension lengths at the endpoints of crack long and short axes. erefore, the initiation characteristics of cracks reflect the growth mode and interaction behavior between cracks to some extent. In order to accurately capture initiation characteristics of the double cracks, the crack initiation stress was measured through multiple pretests and test results indicated that the double cracks initiated when the axial compressive stress reached about 1/3 of the specimen peak strength (σ p � 98∼118 MPa). e ratio was much lower than the initiation stress threshold reported by the previous studies [24,25], which is ranging from 42%∼ 53%σ p . erefore, it is more accurate to observe the crack initiation characteristics when the axial pressure reaches about 1/3 of the peak strength of the specimen. In view of the symmetry of the specimen, the upper crack was mainly analyzed. e measured values of the initiation angles and lengths at the key positions are listed in Table 1, where A 1 and A 2 denote the long-axis endpoints of crack I, and B 1 denotes the short-axis endpoint (see Figure 1(a)).

Test Results and Analyses.
According to test results, the water pressure had significant impacts on the initiation mode and interaction behavior of the double cracks. e test results were measured by photos of the specimen taken when the initial stress was reached. First, the initiation angles at A 1 and A 2 gradually decreased with the increase in water pressure p w , while the initiation angles at B 1 were approximately equal to 0°, which were not listed in Table 1. Second, the growth lengths of the wing crack, defined by the distance from the original crack front to the newly formed crack front, decreased with the increase in water pressure, while the lateral extensions showed a certain increase. e ratio between ∆l B1 and ∆l A1 under different water pressures was 0.020, 0.059, and 0.127, respectively, indicating that the water pressure enhanced the lateral extensions; however, the growth of wing crack was still the leading cause of rock mass failure under hydromechanical coupling [26].
ird, the initiation characteristics of the inner sides (A 2 ) of double cracks were different from those of the outer sides (A 1 ) due to the interaction effect. As shown in Figure 2 and Table 1, the interaction mode was affected by the increase in water pressure. For each crack, the differences between the deflection angles of the two wing cracks increased, while the differences between the growth lengths of two wing cracks decreased.

Establishment of the Numerical Model.
In order to explore the internal mechanism of the interaction between cracks in the rock mass, a numerical model was established using the FEM code of FRANC3D (FRacture ANlysis Code for 3D) V7.0 developed by the America FAC Company, and the corresponding simulations were carried out. e numerical model had the same dimension, crack arrangement, and material properties with the resin specimen. e mesh of the model and crack surfaces is shown in Figure 3. e    Figure 3(b)) was analyzed. e position angle θ (0°≤ θ ≤ 180°) was introduced to define the position of any point at the crack front [1]. e boundary condition of the numerical model was consistent with that of the laboratory test. e vertical displacement of the model bottom was restrained, and the axial compressive stress was applied to the top surface of the model. Besides, the uniform stress was applied to the inner surfaces of cracks to simulate water pressure. e crack growth simulation was accomplished by repeating the following steps: (i) calculation of mixed-mode stress intensity factors (SIFs) and energy release rates (ERRs) based on the M-integration method [27]; (ii) prediction of deflection angle and growth length along the crack front using an appropriate fracture criterion; and (iii) generation of the new crack front. e maximum energy release rate criterion was modified by introducing the weight coefficient corresponding to three basic fracture modes due to their different contribution to the crack growth. e modified maximum energy release rate criterion was employed to determine the crack growth characteristics. e ERR at the crack front is expressed as follows [28]: where E and v are Young's modulus and Poisson's ratio of the material, respectively; K I , K II , and K III represent, respectively, the SIFs for modes I, II, and III fracture; and c I , c II , and c III are the weight coefficients corresponding to K I , K II , and K III , which denote the contribution of three basic fracture modes to the spatial growth of cracks. e growth pattern of the crack is primarily controlled by the angles and lengths at all points along the crack front. e crack initiation angle α and the growth length ∆l i at a certain point i along the crack front are determined as follows [26]: where the crack initiation angle α is obtained when G reaching its maximum value; K r I , K r II , and K r III represent the components of the tensile, slipping, and tearing stresses related to the resolved modes I, II, and III stress intensity factors [28]; η I , η II , and η III are weight coefficients corresponding to the components of tensile, slipping, and tearing stresses; K e i denotes the equivalent stress intensity factor at any point i along the crack front, and K e m represents their midvalue; ∆l m is the growth length corresponding to K e m ; and n is the power index. Besides, in this study, ∆l m is set to 2.4 mm (1/10 of the length of crack long-axis), and n is set to 2.0 to ensure the calculation consistency in this study. Figure 4 shows the initiation mode of double cracks obtained based on the modified energy release rate criterion, and the dashed line denotes the original crack front. e wrapping wing cracks, and a small amount of in-plane extensions were formed at the endpoints of crack long and short axes, respectively, which were similar to the initiation mode observed in the tests (see Figure 2), thus demonstrating the validity of the numerical model.

Verification of the Numerical Model.
Besides, the differences between test results and numerical simulations were quantitatively evaluated by analyzing the growth characteristics at the endpoints of crack long and short axes. Table 2 lists the initiation angles and growth lengths obtained by numerical simulation. When the water pressure increased from 0 to 2.0 and 4.0 MPa, the predicted initiation angle at A 1 decreased from 71.20°to 66.57°and 51.85°and that at A 2 decreased from 71.80°to 67.29°and 55.21°, respectively, while the angle at B 1 (a) remained 0°. e initiation angles had an average error of 0.67% compared with the test results. Besides, the numerically predicted ratios of lengths at the endpoints long and short axes ∆l B1 /∆l A1 were 0.020, 0.062, and 0.134 under different water pressures, which had an average error of 3.53% compared with the test results, further proving the accuracy of the modified criterion.

Case of Parallel Cracks considering the Different Spacings.
Parallel cracks exist commonly in rock masses [28,29], and the crack density controlled by the spacing is a crucial factor affecting the interaction between cracks. In order to investigate the influences of spacing between parallel cracks on their growth and interaction behaviors under hydromechanical coupling, the spacing s is set to 6, 12, 18, 24, 30, and 36 mm, respectively, that is, 0.5∼3.0 times of the half-length of crack long-axis (0.5a∼3.0a), while the water pressure and the crack dip angle remain constant as p w � 2.0 MPa and β 1 , β 2 � 45°. e distribution of stress and fracture parameters along the crack fronts is analyzed, and then, the initiation and interaction modes of parallel cracks with different spacings are clarified.

Distribution of Stress Field.
e interaction between rock cracks is mainly induced by the superposition and disturbance of stress fields around cracks [30]. Two different types of effects, including the interactive enhancement and the interactive weakening, may exhibit, depending on the relative geometric location of cracks. Figure 5 shows the distribution of the stress fields at the longitudinal profiles of the model with different crack spacings (tensile stress is positive). Under the combined action of uniaxial compression and water pressure, the small tensile stress zones (red zones) are formed close to the 3D crack tips, and significant compressive stress zones (blue zones) are formed at both sides of the crack surfaces. When the spacing is small (e.g., s � 6 mm or 12 mm), the tensile stress zones around the crack tips superimpose mutually with the compressive stress zones of the adjacent crack surface, resulting in a significant weakening effect. e superposition of stress fields continually decreases with the increase in crack spacing. When the spacing exceeds 30 mm (2.5a), the stress fields of upper and lower cracks are independent, and the crack interaction is negligible.
Besides, the maximum principal stresses σ max along the crack fronts directly reflect the stress concentration. As shown in Figure 5, the interactive weakening effect mainly presents in the inner sides of parallel cracks, while a certain degree of interactive enhancement effect exhibits in the outer sides of cracks. We suspect that, when the spacing is small, the tensile stress fields at the upper ends of long axes of both cracks overlap with each other, resulting in the interactive enhancement. However, interactive enhancement is weak due to the small range of tensile stress field and disappears when the spacing increases.

Distribution of SIFs and ERRs along 3D Crack Front.
Crack growth from 3D flaws is more complex than the prediction by 2D theories mainly because the I-II-III mixedmode fracture occurs at the 3D crack front, while the I-II mixed-mode fracture occurs at the 2D crack front [10,31,32]. According to the theory of fracture mechanics [28], the distribution and variation in the key fracture parameters, including SIFs and ERRs along the crack fronts, control the initiation modes of cracks and reflect the crack interaction effect to some extent. Figure 6 shows the distribution of SIFs (K I , K II , and K III ) and the ERRs along the front of crack I with different spacing s. Given the symmetry of the model, only the data of half of the crack front corresponding to the position angle θ ranging from 0°to 180°are plotted (see Figure 3(b)). As indicated in Figure 6(a), the K I increases first and then decreases with the increase in θ from 0°to 180°. Compared with the case of single crack [1], the distribution of K I along the double cracks front presents an apparent asymmetry due to the interaction effect. With the decrease in spacing, the K I values along the upper part of the crack front increase, while those of the lower part show a decreasing trend. e distribution of K I value is approximately a symmetry along the crack as the spacing increases to 36 mm, and the maximum value (0.343 MPa·m 1/2 ) occurs at the endpoint of crack short axis (θ � 90°). However, as the spacing decreases to 6 MPa, the K I ranges from 0.264 MPa·m 1/2 to 0.316 MPa·m 1/2 with the minimum value occurring at the crack long axis end in the rock bridge zone (θ �180°). e analysis indicates that the tensile stress of the crack fronts is weakened by the interactive weakening effect when the spacing is small. us, the contribution of mode I fracture to the crack propagation is reduced.   Advances in Materials Science and Engineering e distribution of slipping and tearing fracture modes has great influences on the 3D crack growth but is often neglected in many studies. As shown in Figures 6(b) and 6(c), the numerical results indicate that the distributions of K II and K III along the crack front are significantly affected by the position angle θ. e absolute value of K II decreases to 0  Advances in Materials Science and Engineering from θ � 0°to 90°and then increases approximately in a linear manner from θ � 90°to 180°, whereas the K III values show the opposite trend. Besides, the values of K II and K III are not sensitive to the crack spacing, whereas K II and K III along the lower part of the 3D crack front have a certain degree of reduction when s decreases to 6 mm (0.5a); the compressive stress besides the crack surfaces is superimposed and enhanced, increasing the frictional resistance of the relative slipping of the crack surfaces. e ERR is a comprehensive parameter taking into account the three basic fracture modes. As shown in Figure 6(d), the distribution of ERR with varying crack spacing has approximately parabolic shapes, with their maximum values appearing at θ � 0°or 180°and the minimum values at about θ � 90°, indicating that the cracks are initiated from the endpoints of long axis, which explains the wing crack pattern observed in the tests. Additionally, the ERR at the outer endpoint (θ � 0°) of the long-axis increases with the decrease in s, whereas the ERR at the inner endpoint (θ �180°) shows an opposite tendency.

Initiation Characteristics of Parallel Cracks.
e crack initiation mode reflects the interaction behavior between cracks and then influences the final rock failure under hydromechanical coupling. e initiation angle and length at the endpoints of long and short axes are two crucial indices for assessing the 3D crack initiation modes, and their variation tendencies with the crack spacing are plotted in Figures 7 and 8. Due to the symmetry of the numerical model, the upper crack (crack I) is mainly analyzed. Figure 7 shows the variation in initiation angle at the crack long-axis endpoints. With the increase in spacing, the initiation angles at the inner sides of parallel cracks (α A2 and α A3 ) present a decreasing tendency, while the angles at the outer sides (α A1 and α A4 ) increase slightly. e analysis shows that cracks initiate in the direction where the maximum local tensile stress occurs and then gradually propagate to the direction of far-field compressive stress [27]. In this study, the newly generated cracks at the inner endpoints of the long-axis propagate toward the neighboring crack as spacing decreases, due to the induction effect of its compressive stress fields. Besides, according to the theory of fracture mechanics [28,33], the crack with pure mode I fracture propagates along its original surface; therefore, the initiation angles at the outer endpoints decrease with the increase in crack spacing, due to the increasing contribution of mode I fracture (see Figure 6(a)). Additionally, when the crack spacing reaches a certain value, the initiation angles at inner and outer endpoints of the long-axis of parallel cracks are basically the same, just like the situation of single cracks [1]. Figure 8 shows the variation in growth length at different positions along crack fronts with spacing. e growth of inner endpoints of the long axis (Δl A2 and Δl A3 ) and endpoints of short axis (Δl B1 and Δl B2 ) is weakened when the spacing is small, while that of the outer endpoints of the long axis (Δl A1 and Δl A4 ) is promoted.
Analysis indicates that, when the crack spacing is less than 6 mm (0.5a), the tensile stress concentration of parallel cracks in the rock bridge zone is significantly reduced due to the interactive weakening effects, resulting in the decrease in growth length at the points of A 2 , A 3 , B 1 , and B 2 . However, under the combined action of uniaxial compression and water pressure, the stress concentrations at the points A 1 and A 4 are enhanced, thus improving the crack growth length and leading to that the cracks are easier to initiate at the outer sides of the crack. Similarly, when the spacing increases to 30 mm (2.5a), the growth length tends to a stable value, just like the situation of single crack.

Case of Nonparallel Cracks considering Different Dip
Angles. Nonparallel cracks widely exist in rock masses [34], and their spatial distribution and interaction behavior are much more complex than those of the parallel cracks. In the numerical simulation, the dip angle of the crack I β 1 is set to 15°, 30°, 45°, 60°, and 75°, respectively, while the dip angle of crack II and the crack spacing remain constant as β 2 � 45°a nd s � 12 mm. Based on simulation results, the distributions of stress state and fracture parameters are analyzed and then the initiation and interaction modes of nonparallel cracks with different dip angles are clarified. Figure 9 shows the stress distribution of nonparallel cracks with varying dip angles. e tensile and compressive stress zones are similar to those of the parallel cracks. However, the stress superposition form changes significantly with the increase in β 1 , and the tensile stress zone around the point A 2 deflects gradually to the right side of crack I surface.

Distribution of Stress Field.
e variation in σ max at the key positions of crack fronts with varying dip angles is shown in Figure 10. As depicted in Figure 10(a), the tensile stresses at the inner endpoints of long-axis are much smaller than those at the outer ones and σ A2 decreases significantly with the increase in dip angle β 1 , due to the interactive weakening effects. However, owing to the change in the direction of σ A2 , it is found that σ A2 turns to increase after β 1 exceeds 60°. Besides, σ A3 increases with    Advances in Materials Science and Engineering the increase in β 1 because point A 3 is gradually away from the compressive stress zone besides crack I. As shown in Figure 10(b), the σ max at the endpoints of crack short-axis (σ B1 , σ B2 ) increases significantly with β 1 . e horizontal projected area of the crack I is reduced due to the increase in its dip angle, thus resulting in the reduction of the stressed area of the axial compression, and then, the tensile stress at B 1 increases approximately linearly under the effect of water pressure. Additionally, the tensile stress at B 2 is larger than that of B 1 and continues to increase due to the interactive enhancement effect.

Distribution of SIFs and ERRs of 3D Crack Front.
e distribution of SIFs and ERRs along crack I fronts (0°≤ θ ≤ 180°) with varying dip angles is plotted in Figure 11. Figure 11(a) shows that the K I value is relatively average along the crack front when the dip angle is small, while the maximum value of K I gradually moves to the endpoints of the short-axis with the increase in β 1 . Besides, the K I values along the crack fronts decrease continually with the increase in β 1 due to the interactive weakening effect, especially at the lower endpoint (θ �180°) of crack long-axis. Figure 11(b) illustrates that the values of K II increase with β 1 in the range of 15°≤ β 1 ≤ 45°and then decrease when 45°≤ β 1 ≤ 75°. According to the principle of fracture mechanics, the maximum shear stress exists at the crack fronts when β 1 � 45°; therefore, K II reaches its maximum value at a dip angle of 45°. Besides, the K II value along the lower part of crack fronts when β 1 � 75°and 60°is slightly smaller than those of β 1 � 15°and 30°because the relative slipping deformation of crack surfaces is inhibited when the two cracks are close. e variation trend in K III with the dip angle is similar to K II , as shown in Figure 11(c). Figure 11(d) shows that, with the increase in β 1 , the ERR values first increase and then decrease with its maximum value appearing when β 1 � 45°, indicating that the initiation and growth of the crack with a dip angle of 45°are much easier than other angles. Comparing with the upper half of crack front (0°≤ θ ≤ 90°), the average ERR on the lower half (90°≤ θ ≤ 180°) is larger when Advances in Materials Science and Engineering β 1 � 15°and 30°and smaller when β 1 � 60°and 75°because of the significant weakening effects of the stress superposition. Figure 12 shows the variation in initiation angles at the endpoints of crack long-axis with different dip angles. e simulation results show that the initiation angles of the crack II (α A3 and α A4 ) do not change much; however, those of crack I (α A1 and α A2 ) increase with β 1 in the range of 15°≤ β 1 ≤ 45°, and then decrease when 45°≤ β 1 ≤ 75°. Besides, for each crack, the initiation angle of the inner endpoint of the long-axis is larger than those of the outer endpoint (α A2 > α A1 and α A4 > α A3 ). e variation in growth length at key positions of nonparallel cracks with different dip angles is shown in Figure 13. As indicated by Figure 13(a), the growth lengths at the crack I endpoints A 1 and A 2 (Δl A1 and Δl A2 ) increase with β 1 in the range of 15°≤ β 1 ≤ 45°and decrease with the increase in β 1 when 45°≤ β 1 ≤ 75°, which is similar to the variation trend in ERR. Besides, the growth length at the lower endpoint of crack II (Δl A4 ) almost does not change. However, the growth length of the upper endpoint Δl A3 decreases continually with the decrease in β 1 due to the interactive weakening effect. Crack II is less affected by the interaction effect when β 1 ≥ 60°.

Initiation Characteristics of Nonparallel Cracks.
As shown in Figure 13(b), the growth lengths at the short-axis endpoints (Δl B1 and Δl B2 ) of nonparallel cracks have different degrees of increase but are still much less than those at the endpoints of crack long-axis.

Discussion
From this paper, we can see that the initiation mode and interaction behavior of double cracks are closely related to the water pressure and crack distribution form. As reported by Zhao et al., Zhou et al., and Wang et al. [35][36][37], the mechanism of propagation and evolution of 3D cracks in rock mass under hydraulic coupling is more complicated than that under pure stress. erefore, the existence of water pressure changes the distribution form of stress intensity factor along the crack front and further influences the fracture characteristics, including the deflection angle and growth length. Besides, Zhou et al. [38] and Zhou et al. [39] investigated the initiation, propagation, and coalescence processes of double cracks in rocks and confirmed that the interaction effect of double cracks significantly affected their growth path, which was quite different from those of single crack [38,39]. erefore, the spatial distribution form of double cracks greatly affects the distribution of stress fields around cracks, and the interactive weakening or enhancement effect occurs due to the stress superposition, which has a significant influence on crack initiation and growth characteristics. Moreover, Wong et al., Lu [40][41][42][43][44], which is mainly reflected in the fact that the spatial propagation mode of 3D cracks is more complex, and there are intrinsic limits on 3D propagation of wing cracks, resulting in great change in rock failure mode. In our study, the contribution of mode II and mode III fracture to the 3D crack growth is investigated, and the influence of crack distribution and interaction characteristics is clarified in detail.    In the future, further research should be continued. More comprehensive and systematic studies could be conducted to investigate the interaction effect during the propagation and coalescence processes of double cracks based on the experimental and numerical methods. Further studies on the influence of crack interaction effect on the fracture and mechanical characteristics of rocks will be performed by considering the double cracks with other more complex distribution forms. In addition, in order to more truly reflect the cracking behavior of real rocks under hydromechanical coupling, different kinds of natural rocks, instead of rock-like materials, will be used to carry out experiments on the growth and interaction characteristics of double cracks.

Conclusions
Laboratory tests were performed by using specimens containing double 3D cracks to investigate the initiation and interaction behavior of double cracks under hydromechanical coupling. A numerical model was established and verified by the test results, and the corresponding numerical simulations were carried out. e maximum energy release rate criterion was modified to determine the crack growth characteristics. e distribution of stress fields, key fracture parameters, and crack initiation characteristics was studied in detail. e primary conclusions are as follows: (1) e increase in water pressure leads to the crack initiation being inclined to the original surface and the growth lengths along the crack front tending to be uniform. (2) Under the combined action of axial compression and water pressure, small tensile stress zones are formed close to the crack tips, and significant compressive stress zones are formed at both sides of the crack surfaces. Stress superposition and interaction occur when crack spacing is less than 2.5a. e interactive weakening effect is mainly present in the rock bridge zones, while the interactive enhancement effect is exhibited in the outer sides of cracks. (3) With the increase in spacing, the initiation angles at the inner sides of parallel cracks (α A2 and α A3 ) present a decreasing tendency, while the angles at the outer sides (α A1 and α A4 ) increase slightly. e growth lengths of the outer wing cracks increase, and cracks are much easier to initiate there when the spacing is less than 0.5a. (4) e stress distribution of nonparallel cracks under hydromechanical coupling is more complicated than that of the parallel cracks. With the increase in β 1 , the tensile stress region around the point A 2 deflects gradually to the right side of the crack I surface.
Besides, the maximum principal stresses at the endpoints of crack short-axis (σ B1 , σ B2 ) increase significantly with β 1 . (5) e initiation angles and lengths at the endpoints of crack I long-axis increase with dip angle β 1 when 15°≤ β 1 ≤ 45°and then decrease with β 1 when 45°≤ β 1 ≤ 75°. Cracks with a dip angle of 45°are much easier to initiate at the endpoints of the long-axis.