## 4 Numerical Simulation Methods
### Deviatoric Stress  
Under hydrostatic pressure conditions, the three principal stresses of the surrounding rock are mutually equal, i.e., \(\sigma_{1}=\sigma_{2}=\sigma_{3}\). After roadway excavation, the stress equilibrium  
Figure 8: Schematic diagram of unit energy relationship.  
state of the original rock is disrupted, leading to new stress distribution characteristics where \(\sigma_{1}\geq\sigma_{2}\geq\sigma_{3}\), resulting in:  
\[\begin{bmatrix}\sigma_{1}-\sigma_{m}&0&0\\ 0&\sigma_{2}-\sigma_{m}&0\\ 0&0&\sigma_{3}-\sigma_{m}\end{bmatrix}+\begin{bmatrix}\sigma_{m}&0&0\\ 0&\sigma_{m}&0\\ 0&0&\sigma_{m}\end{bmatrix}=\begin{bmatrix}\sigma_{1}&0&0\\ 0&\sigma_{2}&0\\ 0&0&\sigma_{3}\end{bmatrix} \tag{8}\]  
Here, we introduce the mean stress \(\sigma_{m}\), defined as \(\sigma_{m}=1/3(\sigma_{1}+\sigma_{2}+\sigma_{3})\). The deformation of an object can be primarily attributed to two components: changes in volume and changes in shape. Specifically, the second term in Equation (8) represents the spherical stress tensor that causes volume changes, while the first term describes the deviatoric stress tensor that induces shape changes. To analyze deviatoric stress in the FLAC\({}^{\text{3D}}\) model, after the model has completed its run, pointer functions are utilized to capture the principal stresses \(\sigma_{i}(i=1,2,3)\) of the elements that are mutually perpendicular. Deviatoric stress is then calculated using Equation (9).  
\[\sigma t=\sigma_{1}-\sigma_{m}=\sigma_{1}-\frac{1}{3}(\sigma_{1}+\sigma_{2}+ \sigma_{3}) \tag{9}\]  
Using Fish language to incorporate the calculation of dissipated energy into the FLAC\({}^{\text{3D}}\) model, as well as including the calculation of deviatoric stress and plastic zone, follows a specific operational procedure as illustrated in Figure 9.  
Figure 9: Fish language operation process.  
After the model construction, the entire element mesh is first traversed to assign a unique identifier to each element. Subsequently, the endpoint of the computation is set, and a single cycle operation time step \(m\) is defined. Before computation, the program records the initial stress state \(\sigma_{t}\) of each element. The iterative computation is then initiated, recording the stress change \(\sigma_{t+\triangle t}\) and strain increment \(\epsilon_{t}\) for each element after \(m\) time steps. Based on these data, the energy value for each element is further computed. Once the computation reaches the predefined endpoint, the program automatically assigns the calculated energy values to the corresponding elements and generates a dissipated energy density cloud map. Subsequently, relevant data are exported, and the computation is terminated. If the computation has not reached the predetermined endpoint, the program automatically enters the next cycle of computation, repeating this process until reaching the computation endpoint.  
Upon reaching the endpoint of model operation, pointer functions are utilized to extract the principal stresses \(\sigma_{t}(i=1,2,3)\) of elements that are mutually perpendicular. Deviatoric stress is then calculated, and the calculated deviatoric stress is assigned to the corresponding elements, followed by printing the deviatoric stress distribution map. Additionally, pointer functions are used to extract information regarding element failure, calculating the number of failed elements and the failure ratio within the study area, with the computed data being exported.