## 4 Model Establishment and Instability Evaluation Method
### Boundary Condition Setting  
During the mining process, the initial stress field is generally influenced by both the tectonic stress field and the stress field generated by self-weight. In the FLAC\({}^{3\text{D}}\) model, the internal stress is equilibrated exclusively through redistribution mechanisms. The self-weight of the rock mass is applied to the model in the form of body force. The initial  
Figure 11: Arrangement of seismic sources and anchor.  
Figure 12: Boundary conditions at the source.  
Figure 10: Numerical model geometry.  
stress field obtained by stress redistribution is usually expressed as the self-weight stress field, which is inconsistent with the actual situation. As a result, following the rapid stress boundary fitting technique proposed by [PERSON] et al. [27], the numerical model does not impose a velocity boundary condition. Instead, it applies a stress boundary condition to the model's surface based on the measured ground stress data. The force applied is considered to be the stress acting on the exterior of the model. The initial stress field obtained from the rapid stress boundary fitting technique can be regarded as a combination of the tectonic stress field and the vertical stress field. Finally, the displacement velocity generated after the initial stress field is balanced is cleared, the corresponding static boundary conditions are set, and the roadway excavation is carried out. The relationship between horizontal and vertical in situ stress applied on the surface of the model is defined as follows. The ratio is calculated from the in situ stress test results of 4-5 coal seams:  
\[\sigma_{zz}^{0}=\gamma H \tag{13}\]  
\[\sigma_{yy}^{0}=1.243\sigma_{zz}^{0} \tag{14}\]  
\[\sigma_{xx}^{0}=1.184\sigma_{zz}^{0} \tag{15}\]  
In the formula, \(\gamma\) is the bulk density of the material, and \(H\) is the depth from the surface.  
When conducting dynamic calculations in the excavated roadway, it is crucial to accurately replicate the natural energy dissipation. Rayleigh damping is implemented by utilizing a specified damping ratio at two particular frequencies. Within this frequency spectrum, the model's actual damping ratio is observed to be lower than the true attenuation effect present at the site. When the frequency exceeds this range, the damping will increase significantly, and attenuation inconsistent with the field vibration wave will occur. [PERSON] et al. [19] conducted a validation of the damping scheme using three Maxwell components, as proposed by earlier researchers, as an alternative to Rayleigh damping. The findings suggest that, in the majority of seismic deformation analyses, this methodology is effective in producing a relatively stable damping effect within the designated frequency range. In dynamic analysis, it is essential to modify the static boundary conditions employed in static analysis to viscous boundary conditions in order to minimize the reflection and refraction of vibration waves [28]. Under viscous boundary conditions, the normal and tangential viscous forces produced by the damper are determined using the following equation:  
\[t_{n}=-\rho C_{p}v_{n} \tag{16}\]  
\[t_{s}=-\rho C_{s}v_{s} \tag{17}\]  
In the formula, \(v_{n}\) is the normal velocity component at the model boundary, and \(v_{s}\) is the tangential velocity component.