## 2 Establishment of Numerical Model with Non-Parallel Overlapping Flaws  
### Selection of Contact Model  
The core theory of PFC2D is the contact model of particle interaction, in which the parallel bond model (PBM) is commonly used to simulate the mechanical properties of rock-like materials ([PERSON], 2019; [PERSON], 2019). In the PBM, forces and moments can be transferred by applying a bond between particles, which is similar to the distribution of a group of springs with a certain stiffness centered on the contact plane of particles, as shown in Figure 1. After bonding, once the relative motion of particles occurs, the force and bending moment will be generated and the stress state of the medium will be changed. According to the different contact forces, the bonding forms between particles can be roughly divided into two categories: one is the contact bonding that can transfer the force between particles, and the other is the parallel bonding that can not only transfer the force but also the torque. The contact bond between particles can be decomposed along the normal and tangent directions of the contact surface. \[F_{\text{i}}=F^{\text{n}}n_{\text{i}}+F^{\text{s}}f_{\text{i}} \tag{1}\]  
where \(F^{\text{n}}\) and \(F^{\text{s}}\) are normal contact force and tangential contact force, respectively, and meet the conditions of \(F^{\text{s}}\leq\xi F^{\text{n}}\); \(\xi\) is the friction coefficient between particles; \(n_{\text{i}}\) and \(t_{\text{i}}\)are the unit normal vector and tangential vector along the contact surface, respectively. Figure 1: Cohesion model and its micromechanical behavior in the PBM ([PERSON], 2019; [PERSON], 2019). Each subsequent relative displacement and rotation increment (\(\Delta U^{\text{n}}\), \(\Delta U^{\text{s}}\), \(\Delta\theta^{\text{n}}\), \(\Delta\theta^{\text{s}}\)) produces an increment of elastic force and moment, which can be calculated from Equation (2) and is added to the current values. \[ \begin{cases}\Delta\overline{F}^{\text{n}}=\overline{k}^{\text{n}}A\Delta U^ {\text{n}}\\ \Delta\overline{F}^{\text{s}}=\overline{k}^{\text{s}}A\Delta U^{\text{s}}\\ \Delta\overline{M}^{\text{n}}=-\overline{k}^{\text{s}}J\Delta\theta^{\text{n}} \\ \Delta\overline{M}^{\text{s}}=-\overline{k}^{\text{n}}I\Delta\theta^{\text{s}} \end{cases}  \tag{2}\]  
where \(A\) is the area of the cross-section parallel to the bond; \(I\) and \(J\) are the moment of inertia and polar moment of inertia of the parallel bonded cross-section; \(\Delta U^{\text{n}}\) and \(\Delta U^{\text{s}}\) are relative displacement increments in normal and tangent directions respectively; \(\Delta\theta^{\text{n}}\) and \(\Delta\theta^{\text{s}}\) are the relative corner increments in the normal and tangent directions respectively. In addition, if the specimen is not subject to tensile shear failure during loading, the maximum tensile stress and tangential stress applied to the parallel bond shall also meet the following conditions:  
\[\overline{\sigma}^{\text{max}}=\frac{-\overline{F}^{\text{n}}}{A}+\overline{ \beta}\frac{|\overline{M}_{\text{b}}|\overline{R}}{I}<\overline{\sigma}_{ \text{c}} \tag{3}\]  
\[\overline{\tau}^{\text{max}}=\frac{\left|\overline{F}^{\text{s}}\right|}{A}+ \overline{\beta}\frac{|\overline{M}_{\text{d}}|\overline{R}}{J}<\overline{ \tau}_{\text{c}} \tag{4}\]  
where \(\overline{M}_{\text{b}}\) and \(\overline{M}_{\text{t}}\) are bending moment and torque of parallel bonding interface; \(\overline{R}\) is the particle radius used in the calculation; \(\overline{\sigma}_{\text{c}}\) and \(\overline{\tau}_{\text{c}}\) are ultimate compressive strength and shear strength, respectively; \(\overline{\beta}\) is the moment contribution factor. In PFC2D, the parallel bond model can clearly reflect the real failure situation of rock or rock-like materials. Therefore, the parallel bond model was used to simulate the samples with non-parallel overlapping flaws in this study.