### Literature Review  
Discrete methods have been extensively explored to address the limitations of continuum-based numerical models in tunnel-stability analyses. Among these, lattice models [37; 51; 52; 53] represent materials as interconnected networks of mechanical elements such as springs or bars, where fractures propagate through element failure processes. While advantageous for simulating brittle fracture and microstructural mechanics [37], lattice models often impose fixed Poisson's ratios, rely on external failure criteria, and exhibit sensitivity to lattice topology, thereby affecting numerical robustness. Efforts to refine lattice-based models include the incorporation of beam elements and angular spring interactions to improve material response accuracy [38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50]. The Discretized Virtual Internal Bond (DVIB) method [45] introduced a more adaptable numerical framework by discretizing materials into bonded cell structures governed by hyperelastic potential functions. Unlike conventional lattice models, the DVIB enables the direct simulation of large deformations and fracture evolution without pre-defined failure criteria. However, early DVIB formulations relied on two-body interactions, resulting in fixed Poisson's ratios that restricted the ability of the model to represent anisotropic deformation patterns. [PERSON] et al. [54] addressed this constraint by integrating rotational effects via bond angles, improving Poisson's ratio representation in fracture mechanics applications. For tunnel stability analysis, numerical models must capture elastoplastic deformation and fracture propagation to ensure realistic stress redistribution within rock formations. Although DVIB extensions have introduced plasticity [50], the issue of a fixed Poisson's ratio persists, limiting the adaptability of modeling stress-dependent behaviors in underground environments. [PERSON] et al. [55] extended the Stillinger (SW) potential framework, enriching the constitutive relationship with microfrac mechanisms and improving the numerical accuracy of failure simulations. In addition to discrete models, multiscale simulation methods such as the quasi-continuum (QC) approach [56; 57], Coarse-Grained Molecular Dynamics (CGMD) [58], and bridge-scale (BS) models [59; 60] aim to bridge microstructural mechanics with macroscopic deformation phenomena. Although effective in representing material behavior across different scales, these techniques often incur high computational costs and suffer from interface artifacts, such as ghost forces, which affect numerical precision [61]. Continuum-enhanced frameworks, such as Extended Finite Element Methods (XFEM) [27] and meshless methods [29], resolve remeshing challenges, but they require higher computational complexity and increased degrees of freedom for accuracy [28]. Recent studies, including [PERSON] et al. [62], applied the plastic Stillinger-Weber DVIB (PSW-DVIB) model to tunnel deformation assessments, investigating lateral pressure coefficients, in situ stress effects, Poisson's ratio variations, and tunnel geometry-dependent responses. Their findings reinforced the applicability of the DVIB-based models for fracture evolution and stress redistribution simulations. However, their approach implemented ideal plasticity, which did not achieve explicit decoupling between plastic flow and lateral stiffness control, a critical factor in tunnel deformation modeling under complex geological conditions. This study builds upon a previously published work by [PERSON] et al. [63], which introduced an elastoplastic SW-DVIB model incorporating variable Poisson's ratios within the Stillinger-Weber potential framework. This study successfully embedded plasticity within the two-body interaction while maintaining the elastic behavior in the three-body component, enabling the simulation of fracture propagation and irreversible deformation in brittle materials.