**1. Introduction**

Numerous discontinuities on different scales occur in rocks and rock masses, such as flaws, joints, cracks, and faults, which are induced during the formation of rock masses and successive tectonic processes. The existence of joints, cracks, or faults in rock masses has two mechanical effects [1–4]: (1) The strength and stiffness of the rock mass will be decreased; and (2) they will generate new fractures that may further develop and connect with other flaws and then may lead to the degradation of the strength and stiffness of the rock or rock mass, respectively. Therefore, the mechanism of crack propagation and coalescence is of grea<sup>t</sup> interest in geomechanics [5,6].

To further study the laws of cracks and the development of single cracks, physical experiments were conducted [7–18]. The observed crack patterns under uniaxial compression can be categorized as follows: (1) Wing cracks initiated at the tips of single, pre-existing flaws propagating along the direction of axial stress; and (2) secondary shear or mixed mode tensile-shear cracks initiated later from the tips of single pre-existing flaws propagating as quasi-coplanar or oblique cracks.

Crack coalescence mechanisms were studied by other researchers with rock-like materials containing two or more parallel flaws in a compression state [1,19–27]. It was observed that flaw angle, bridging angle, non-overlapping length, and frictional coefficient have significant effects on the crack initiation modes and the patterns of crack coalescence. Moreover, it was shown that the peak strength and stiffness of rocks and rock-like material depends on the number and geometry of flaws [28].

Previous work provides a good understanding of the coalescence patterns obtained from specimens with one or several single, non-intersecting flaws. However, joint patterns with more than one joint set are common in nature. Pollard and Aydin [29] defined several types of joint intersection geometries, which can be classified as orthogonal (+intersections) and non-orthogonal (X intersections). Both types can be divided into three groups according to the continuity of the joints at the intersections: (1) "all continuous joints"; (2) "continuous and discontinuous sets"; and (3) "all discontinuous sets".

Lee and Jeon [30] used three materials (PMM, Diastone and Hwangdeung granite) to study the coalescence of a horizontal flaw with an underneath lying, inclined flaw. They stated that the geometry of a double-flaw can improve the understanding of crack propagation and coalescence because echelon type of cracks can be initiated by fracture planes that are not parallel to each other. Zhang et al. also studied the crack coalescence of gypsum with two non-parallel flaws under di fferent flaw angles and di fferent bridge angles. Their results showed that the stress distribution in the bridge area of the non-parallel flaws is more complicated than that of the parallel flaws. This di fference a ffects the crack initiation stress as well as the coalescence pattern. Their studies provided fundamental understanding of the coalescence of non-parallel flaws, but observations made from two non-parallel flaws can only provide a limited understanding of the complex behavior of the rock mass with several crossing discontinuities. Besides this, numerical simulation also plays an important role in crack coalescence research [1,22,31–33]. One of the popular tools nowadays is the discrete element method [18,34–40].

In the present study, numerical simulation results are shown in respect to peak strength, crack initial stress, crack coalescence stress, and coalescence patterns of rock-like material with two intersecting flaws under uniaxial compression using a bonded-particle model (BPM).

### **2. Numerical Simulation Procedure**

We used the bonded particle model (BPM), which is one of the contact models available in the two-dimensional particle flow code (PFC2D) software for numerical simulation. Turco [41] suggested a micro-mechanical model, which completely is described by the definition of the strain energy of the interaction between two nearby grains without consideration of the dynamical e ffects. This model can also e ffectively handle large deformation processes, even if quasistatic simulations require the solution of remarkable numerical challenges because they might present instabilities and bifurcation phenomena [42]. In the present study, numerical simulation focuses on the peak strength, crack initial stress, crack coalescence stress, and coalescence patterns of rock-like material with two intersecting flaws under uniaxial compression, thus the bonded particle model (BPM) could satisfy this study.

In PFC (such as that in BPM), the model specimens are created by combining a large number of circular particles with varying diameters. These particles are interconnected using parallel bonds, as shown in Figure 1. Tensile and shear stresses acting on the parallel bonds are calculated during simulations and if the maximum tensile stress exceeds the tensile strength of the parallel bonds or the maximum shear stress exceeds the shear strength of the parallel bonds, the parallel bonds will break and a microscopic tensile or shear crack will form (Figure 1).

To resemble the real cracking processes, discrete cracks that are close enough are connected with each other by a unique technique [18]. The symbol "a" stands for the centroids' distance form a crack to the other one nearby, and the symbol "c" stands for the length of any of these two cracks (Figure 2). If a/c is less than or equal to 1, the crack, and the other one nearby, are regarded as a single continuous crack, and then the centroids of the two micro-cracks will connect with each other as a macroscopic crack trace line [18].

**Figure 1.** Illustration of the parallel bond model.

**Figure 2.** Constructing a macro-crack based on connecting centroids of micro-cracks [18].

The geometry and parameters of the model specimens were chosen in such a way that the model duplicates the mechanical behavior of the gypsum specimens used in the experimental investigation of Wong and his subsequent numerical simulations [43]. The model specimens were created with the dimensions 50 mm in width and 100 mm in height, which are consistent with those of the specimens used for the laboratory testing by Wong [25]. The model specimens created in this manner comprised about 19,650 particles. The porosity of the model specimens was found to be 0.01.

The simulation program consisted of two stages. First, model specimens free of pre-existing flaws (i.e., intact) were simulated. Micro-parameters for the particles were selected in order to obtain a consistent macro-scale mechanical behavior with that of the gypsum specimens used by Wong [25]. Table 1 outlines these micro parameters assigned for the particles and parallel bonds of the numerical specimens independently. The next section reports the comparison of the mechanical properties of the numerical specimens of the present study and those of the specimens used for laboratory testing by Wong [25].

The second stage of the study involved a series of numerical simulations on model specimens containing two intersecting flaws with varied geometry to explore the influence of various geometries on the mechanical behavior of model specimens; flaw length L and thickness were fixed at 20 and 1 mm for both flaws, respectively, while the intersection angle α of the two flaws varied at 90◦, 60◦, 45◦, and 30◦ (Figure 3). The inclination angle of one flaw (angle β of flaw B as shown in Figure 3b) was also varied in the next stage, while α was kept unchanged. To better distinguish the crack coalescence types and for the clarity of the explanations, the bridging zone (the zone surrounding the intersection of flaws) was divided into four zones, as shown in Figure 3c (I, II, III, and IV).


**Table 1.** Micro-parameters used for the model specimens.

**Figure 3.** Geometries of the model containing two intersecting flaws: (**a**) Flaws are numbered A and B; (**b**) flaw inclination angle β, intersection angle α and flaw length L; and (**c**) four bridging zones.

For all the simulations in the present numerical study, specimens were vertically loaded in a velocity-controlled manner and no confining pressure was applied (i.e., uniaxial compression). A sufficiently low loading rate −0.05 m/s was applied to ensure that the specimens remained in a quasi-static equilibrium throughout the test.

### **3. Results and Discussion**

### *3.1. Mechanical Properties of Intact Models*

We first determined the macro-scale mechanical properties of the model specimens and compared them with the corresponding properties reported by Wong [43]. Table 2 and Figure 4 show this comparison. It can be seen that the numerical model has very closely reproduced the mechanical behavior of gypsum specimens used in the laboratory investigation by Wong [43]. Following this verification, we simulated the mechanical behavior of specimens containing two intersecting flaws and the forthcoming sub-sections report the results of those simulations.


**Table 2.** Lab test results [43] vs. numerical simulation results.

**Figure 4.** The stress-strain curves: Lab test vs. numerical simulation results.

### *3.2. Mechanical Behavior of Specimens Containing Flaws*

### 3.2.1. Effect of the Flaw Geometries on Peak Strength

The influence of flaw geometry was investigated in regards to the strength of the sandstone samples, gypsum, and marble containing single crack, parallel cracks, or non-parallel cracks in uniaxial loading [20,21,44,45]. Their results revealed that flaw geometry has a key effect on the strength. Thus, variations of peak strength against β for different values of α are shown in Figure 5.

**Figure 5.** Variations of peak strength against β for different values of α.

From Figure 5, the following general trends can be deduced:


To explain this, the peak strengths of the model containing single crack against β for different initial angles are shown in Figure 6.

**Figure 6.** The peak strengths of the model containing a single flaw against β for different flaw initial inclination angles. -1 -2 -3 -4 -5 curves corresponding to five flaw initial inclination angles in Figure 7.

**Figure 7.** Five flaw initial inclination angles of the model containing a single flaw. The red dotted line represents the final state of rotation for the flaw. The blue arrow represents the direction of rotation.

The overall strength of a rock mass containing two sets of discontinues is given by the lowest strength envelope in the individual strength curves, according to Bray's suggestion [46]. It is clear from Figure 6 that the strength curves of the model containing intersect cracks had a similar trend with the lowest strength envelope in the individual strength curves. For example, the lowest strength envelope in the individual strength curves -1 -5 show the strength of the model containing two intersect cracks (α = 90◦); the peak strength first increases up to about β = 45◦ and then decreases with the increase of the inclination angle. It is worth noting that the lowest strength envelope in the individual strength curves -1 -4 are a little bit different (when β = 60◦) from the trend of α = 60◦ in Figure 5. One plausible explanation is that it was due to the two cracks effecting each other. Another plausible explanation is that it was due to a random distribution of particle sizes.

### 3.2.2. Crack Initiation and Coalescence Behavior

The crack evolution of models containing two intersecting flaws (α = 90◦, 60◦, 45◦, and 30◦) is shown in Table 3. Five stages corresponded to five points on the stress–strain curve in Figure 8: Document crack initiation, propagation, coalescence, as well as crack type at peak stress or post-peak stress level. To shorten and avoid repetition, only selected constellations of α and β are documented in detail.

**Table 3.** The crack evolution pattern at different stages of loading with different flaw geometry. The cracks corresponding to five points on the stress–strain curve in Figure 8. Microscopic tensile and shear cracks are shown in white and red, respectively. Original flaws are in yellow.

**Figure 8.** Axial stress (black line) and total crack number (blue line) versus axial strain with different values of α and β. Five points, A, B, C, D, and E, on the stress–strain curve are monitored to observe the crack initiation, propagation, and coalescence.

For flaw geometry with α = 90◦ and β = 30◦, macroscopic wing cracks initiate from the tips of the A flaw. With an increasing load, the wing cracks extend, and macroscopic tensile cracks appear in the tip area of the B flaw as secondary cracks. When axial stress reaches 18.8 MPa, the first coalescence is achieved in bridge zone II by the extension of the secondary crack starting from the right tip of the B flaw. During a further load increase up to the peak strength, a downward extension of the macroscopic tensile crack, starting from the left tip of the B flaw, reaches the left tip of the A flaw, and the second coalescence occurs in bridge zone IV. The test specimen is broken by the development of cracks during the post-peak stress stage. The coalescence cracks in the bridge region become progressively wider during the cracking processes.

For flaw geometry with α = 60◦ and β = 0◦, and when axial stress has reached 4.8 MPa, macroscopic wing cracks start to propagate from the tips of the A flaw. With a further load increase, the wing cracks extend and some microscopic cracks appear in the tip area of the B flaw (secondary cracks). At 18 MPa vertical stress, the first coalescence is achieved in bridge zone III by the extension of the secondary crack initiated from the right tip of the B flaw. At peak stress of 19.5 MPa, coalescence cracks become wider and some new microscopic tensile and shear cracks propagate downwards or upwards from the

tips of the B flaw. Beyond peak stress, some secondary cracks extend quickly towards the edge of the model, and finally the model fails.

For flaw geometry with α = 45◦ and β = 45◦, macroscopic tensile cracks start from the tips of the A flaw. With an increasing load, wing cracks extend and some microscopic tensile cracks are initiated at the tips of the B flaw. With a further stress increase, secondary tensile and shear cracks propagate towards the wing crack initiated at the top tip of the A flaw. At 21.1 MPa, the first coalescence occurs in bridge zone I. At peak stress of 23.9 MPa, the secondary cracks propagate further upwards from the right tip of the B flaw towards the top tip of the A flaw. Secondary coalescence occurs in bridge zone II. Additional cracks are induced close to the B flaw tip regions in the post-peak stage.

For flaw geometry with α = 30◦ and β = 0◦, at 5.4 MPa, a macroscopic wing crack starts to propagate from the tips of the A flaw. Further load increases lead to the initiation of micro tensile cracks at the right tip of the B flaw. When the axial stress increases further up to 17.6 MPa, the first coalescence is achieved in bridge zone III by the extension of a secondary crack initiated from the right tip of the B flaw. At peak strength, coalescence occurs in bridge zone I. Meanwhile, some mixed mode tensile–shear cracks appear and extend downwards near the right tips of the A and B flaws. In the post-peak region at a stress of 17.4 MPa, some secondary cracks become wider and extend to the edge of the model and the model fails.

### *3.3. Crack Initiation Stress*

In the numerical studies performed by Zhang et al., wing cracks start to propagate simultaneously at the tips of the flaws. However, in this study, the macro-cracks initiated at the tips of the two intersecting flaws do not occur simultaneously. Therefore, the stress state when a macro-crack starts to develop from either the A flaw or B flaw is called crack initiation (CI) stress.

Specimens with di fferent flaw geometries have di fferent peak strengths; therefore, CI stress cannot be compared directly. Thus, the CI stress σ*ci* is normalized by the respective peak strength σ*cn* (see Figure 9). From Figure 9, the following general trends can be deduced:


**Figure 9.** <sup>σ</sup>*ci*/<sup>σ</sup>*cn* versus flaw inclination angle β for di fferent values of α.

To better understand the effect of α and β in the ratio <sup>σ</sup>*ci*/<sup>σ</sup>*cn*, the differences (Δ <sup>σ</sup>*ci*/<sup>σ</sup>*cn*) between the maximum and minimum values of <sup>σ</sup>*ci*/<sup>σ</sup>*cn* are calculated, as shown in Figure 9. From Figure 10a, the difference of <sup>σ</sup>*ci*/<sup>σ</sup>*cn* is first decreased down to α = 45◦ and then increased with the increase of α, which indicates that the inclination angle β had less effect on the ration <sup>σ</sup>*ci*/<sup>σ</sup>*cn* when α = 45◦.

**Figure 10.** Variations of Δ <sup>σ</sup>*ci*/<sup>σ</sup>*cn* against α (**a**) and β (**b**).

From Figure 10b, it can be deduced that the difference of <sup>σ</sup>*ci*/<sup>σ</sup>*cn* also decreases down to β = 30◦ and then increases with the increase of β, which indicates that the intersection angle α had less of an effect on the ration <sup>σ</sup>*ci*/<sup>σ</sup>*cn* when β = 30◦.

### *3.4. Crack Coalescence*

Many crack coalescence types have been observed in physical experiments [1,21,26] and have been classified as seven different crack coalescence types, based on their geometry and propagation mechanism using a high-speed camera. Zhang [40] classified five types of linkages observed between two non-parallel flaws using parallel bonded-particle models. However, these crack coalescence types are based on the coalescence between two parallel or non-parallel flaws. Our study reveals that two intersect flaws show different coalescence behavior. Table 4 shows the crack coalescence patterns for two intersecting flaws with different values for α and β.

The different crack coalescence patterns observed in the case of two intersecting flaws with different geometries are shown in Figure 11 and are summarized in Table 5. Two typical patterns can be distinguished:



**Table 4.** Crack coalescence patterns for two intersecting flaws with different values of α and β. Microscopic tensile and shear cracks are shown in white and red, respectively. Original flaws are in yellow.

**Figure 11.** Observed crack coalescence patterns for two intersecting flaws (left: Simulation results, right: Simplified sketch). White and red segments represent microscopic tensile and shear cracks, respectively. The arrows represent the crack propagation direction.


**Table 5.** The coalescence patterns of different geometries.

### *3.5. Crack Coalescence Stress*

The coalescence stress σ*co* is the vertical stress that is observed when the first coalescence of the two intersecting flaws occurs. The ratio between crack coalescence stress σ*co* and peak strength σ*cn* varies, as shown in Figure 12. The crack coalescence stress reaches 85%–100% of the peak stress.

**Figure 12.** Ratio between coalescence stress σ*co* and peak strength σ*cn* for different values of α and β.

For two cases of α = 90◦ and one case of α = 60◦, α = 45◦, and α = 30◦, the ratios of their cases are at or just less than 1.0, which means that their model coalesced near or at the peak stress. In contrast, for four cases of β = 45◦, three cases of β = 90◦, and one case of β = 0◦, the ratios of their cases are significantly less than 1.0, which means that the coalescence occurs more easily. In addition, for any intersection angle α, the minimum values of <sup>σ</sup>*co*/<sup>σ</sup>*cn* are observed at β = 45◦, which indicates that coalescence occurred more easily with that geometry.
