## 2 Establishment and Calculation Process of a Numerical Model
### Setting of Constitutive Model for Numerical Simulation  
In this study, the strain softening model is selected as the mechanical model for the numerical simulation depicted in Figure 1. This model builds on the Mohr-Coulomb (M-C) framework with non-correlated shear and correlated tensile flow rules. Unlike the M-C model, the parameters of the strain softening model undergo softening or hardening after the specimen yields. The M-C yield criterion determines rock mass failure as follows [39]:  
\[\left\{ \begin{array}{l}f_{s}=\sigma_{1}-\sigma_{3}\frac{1+\sin\phi}{1-\sin \phi}-2c\sqrt{\frac{1+\sin\phi}{\Gamma-\sin\phi}}\\ f_{t}=\sigma_{3}-\sigma_{t}\end{array} \right. \tag{1}\]  
where \(\sigma_{1}\) and \(\sigma_{3}\) are the maximum and minimum principal stresses, respectively; \(c\) and \(\varphi\) are the cohesion and internal friction angle, respectively. When \(f_{s}>0\), the material will undergo shear failure. Under normal stress conditions, the tensile strength of rock mass is very low, so the tensile strength criterion (\(f_{s}>0\)) can be used to determine whether the rock mass has undergone tensile failure [40]. In the strain softening model, the incremental form of principal stress and strain is expressed by Hooke's Law as follows:  
\[\left\{ \begin{array}{l}\Delta\sigma_{1}=\alpha_{1}\Delta e_{1}^{e}+\alpha_{ 2}\big{(}\Delta e_{2}^{e}+\Delta e_{3}^{e}\big{)}\\ \Delta\sigma_{2}=\alpha_{1}\Delta e_{2}^{e}+\alpha_{2}\big{(}\Delta e_{1}^{e}+ \Delta e_{3}^{e}\big{)}\\ \Delta\sigma_{3}=\alpha_{1}\Delta e_{3}^{e}+\alpha_{2}\big{(}\Delta e_{1}^{e}+ \Delta e_{2}^{e}\big{)}\end{array} \right. \tag{2}\]  
Among these, \(\Delta e_{1}^{e}\), \(\Delta e_{2}^{e}\), and \(\Delta e_{3}^{e}\) represent the incremental expressions of the principal strain, while \(\alpha_{1}\) and \(\alpha_{2}\) are material constants defined by the shear modulus \(G\) and bulk modulus \(K\), as follows:  
\[\left\{ \begin{array}{l}\alpha_{1}=K+\frac{4}{3}G\\ \alpha_{2}=K-\frac{2}{3}G\end{array} \right. \tag{3}\]  
The shear yield function expression for the strain softening model is given as follows:  
\[F^{\text{sy}}=\sigma_{1}-\sigma_{3}N_{\varphi}+2c\sqrt{N_{\varphi}} \tag{4}\] \begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \multirow{2}{*}{**Layers**} & \multirow{2}{*}{**Density \(\rho\)/kg/m\({}^{3}\)**} & \multicolumn{2}{c}{**Elastic Modulus**} & \multirow{2}{*}{**Poisson’s Ratio**} & \multirow{2}{*}{**Cohesion \(\sigma\)/MPa**} & \multicolumn{2}{c}{**Internal Friction**} & \multicolumn{1}{c}{**Tensile Strength**} \\  & & \(E\)/GPa & & & & **Angle \(\varphi^{\text{sy}}\)** & \(\sigma_{\text{t}}\)/MPa \\ \hline Silstone & 2800 & 9.95 & 0.24 & 3.8 & 31 & 0.34 \\ Fine sandstone & 2750 & 10.17 & 0.25 & 3.5 & 32 & 0.43 \\ Mudstone & 1780 & 1.7 & 0.35 & 1.2 & 33 & 0.20 \\ Coal & 2350 & 1.7 & 0.21 & 2.3 & 25 & 0.28 \\ Medium-grained sandstone & 2900 & 16 & 0.30 & 5.3 & 28 & 0.58 \\ Shaft & 2400 & 35.2 & 0.24 & 8.6 & 52.0 & 4.10 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Mechanical parameters for numerical calculation of strata and shaft [8]. where \(N_{\varphi}=\frac{1+\sin\varphi}{1-\sin\varphi}\). The equation for tensile yield strength is expressed as follows:  
\[F^{\mathrm{ty}}=\sigma_{1}-\sigma_{3} \tag{5}\]  
It is generally assumed that the total strain increment can be decomposed into elastic and plastic components. The flow law for plastic yielding is expressed as follows:  
\[\Delta e_{i}^{p}=\lambda\frac{\partial g}{\partial\sigma_{i}} \tag{6}\]  
where \(i=1\), 3: \(\lambda\) represents the plasticity coefficient to be determined, \(g\) denotes the potential function, and its shear yield component is \(g^{\mathrm{s}}\). The corresponding non-correlation equation is given as follows:  
\[g^{\mathrm{s}}=\sigma_{1}-\sigma_{3}N_{\psi} \tag{7}\]  
where \(N_{\psi}=\frac{1+\sin\varphi}{1-\sin\psi}\), \(\psi\) is the shear dilation angle.