## 4 Stress Distribution around Opening  
As stated above, the initiation position and sequence of the cracks, especially the elementary tensile cracks, rely on the surrounding stress around the opening under loads. To grasp the mechanism of crack development, a numerical study was conducted to simulate the stress states around the openings under uniaxial compressive loads. Since the crack occurred in stage III is related to the surrounding stress of the opening, it is meaningful to analyze the surround stress of the opening at the elastic deformation stage. Therefore, this problem is simplified into solving the stress distribution at the periphery of the opening in elastic plate under a constant uniaxial compressive stress. In this section, FLAC software was used for numerical simulation. Considering the difficulty in numerical modeling, the Midas/GTS software was \begin{table}
\begin{tabular}{c c c c c} \hline
**Sample No.** & \(\sigma_{\rm c}\) **(MPa)** & \(\sigma_{\rm i}\) **(MPa)** & \(\sigma_{\rm d}\) **(MPa)** & \(\sigma_{\rm p}\) **(MPa)** \\ \hline C-3 & 24.2 & 44.8 & 92.2 & 101.3 \\ S-2 & 14.0 & 21.2 & 55.2 & 78.5 \\ H-3 & 12.0 & 22.9 & 55.1 & 77.3 \\ D-1 & 10.2 & 18.8 & 53.6 & 68.7 \\ V-3 & 11.1 & 24.6 & 57.3 & 88.8 \\ \hline \end{tabular}
\end{table} 
Table 4: Threshold stress values involved in rock fracture subjected to uniaxial loading. used to build the model, and then imported it into FLAC software for post-processing. To reduce boundary effect, the sizes of the numerical model were determined as 150.0 m (length) \(\times\) 0.5 m (thickness) \(\times\) 150.0 m (height). They were more than ten times the maximum dimensions of the opening. To accurately reflect the true stress states, by means of grid seeding, the dimension of the grid near the opening was 0.1 m, while that away from the opening was 5 m.  
To study only the elastic stress distribution, both the shear and bulk moduli of the model were given a large value, and the elastic model was selected as the rock constitutive model. Moreover, the gravity of the model was neglected. The velocity of the surrounding surfaces was fixed along their axial direction, and the exerted compressive stress was 20 MPa. During modeling, twenty-four monitoring points with 30\({}^{\circ}\) intervals were uniformly arranged on the boundary of the opening in the counterclockwise direction. The polar angle of the starting point P1 was 15\({}^{\circ}\), while that of the last point P24 was 360\({}^{\circ}\). Figure 9a shows the layout of the monitoring points in the model with an inverted U-shaped opening. Based on the simulated stresses of the model along \(x\)-direction and \(y\)-direction, the tangential stresses of the monitoring points could be solved. In fact, if the shape of the two holes is exactly the same, the surrounding stress distributions are identical. Our previous research shows that the stress distribution round the inverted U-shaped opening can be acquired utilizing the complex variable approach [59]. As the shape of the opening in this study was the same as that in our previous research, the hoop stress concentration factors of the 24 monitoring points could be obtained directly. The good consistency between the numerical and theoretical results in Figure 9b indicates that the FLAC software was feasible and reliable for stress simulation. It is pointed out that the discrepancy might be related to the grid size. In fact, we could only monitor the stress of the element in FLAC software; that is, the obtained stress resides in the center of the element rather than the monitoring points. Thus, the smaller the dimension of the element, the closer the two points are, but the computer will take a long time to calculate [60]. Likewise, the hoop stress concentration factors of monitoring points on the boundaries of the two inverted U-shaped openings with different configurations are shown in Figure 10. Based on the stress distributions round the openings, the crack growth and coalescence mechanism were clearly revealed in the Discussion Section.