## 4 Discussion  
In this study, an elastoplastic Discretized Virtual Internal Bond (DVIB) model embedded within a modified Stillinger-Weber (SW) potential was used to investigate tunnel deformation in brittle and quasi-brittle geomaterials. The model introduces a bilinear two-body force-displacement law that governs plastic deformation, permitting bond stretching beyond the yield prior to rupture. In conjunction, a three-body angular stiffness component controls the lateral expansion and volumetric deformation, thereby capturing the Poisson's ratio effects. Together, these features allow for the independent regulation of normal and shear responses, effectively representing the coupling between the plastic flow and anisotropic fracture evolution in brittle rocks. Numerical simulations showed that the SW-DVIB model successfully captured plastic deformation, crack initiation and propagation, and stress redistribution under varying in situ stress conditions and Poisson ratios. The radial displacements of the tunnel wall increased with both vertical and lateral stress magnitudes and aligned well with excavation-induced damage patterns. Specifically, displacements increased from 0.45 mm at \(\
u=0.17\) to 1.32 mm at \(\
u=0.35\), confirming the amplifying role of Poisson's ratio in deformation response. The integration of fracture energy calibration ensured mesh-independent failure behavior, lending physical realism and numerical robustness to the model. Figure 24: Tunnel deformation with diverse Poisson’s ratios: (**a**) \(\
u=0.2\); (**b**) \(\
u=0.25\); (**c**) \(\
u=0.3\). (Nodal displacements have been magnified by a factor of 20 for enhanced visualization). Figure 23: Displacement contour of model with diverse Poisson’s ratios: (**a**) \(\
u=0.2\); (**b**) \(\
u=0.25\); (**c**) \(\
u=0.3\). A significant contribution of this work is the explicit inclusion of pre-crack inclination effects via Compact Shear Cracked Brazilian Disk (CSCBD) simulations. As detailed in Section 3.5, the model accurately captured the dependence of fracture toughness on the crack angle, corroborating the findings of [PERSON] et al. [69]. Their machine learning-based study showed a minimum fracture toughness at specific crack angles owing to the combined influence of tensile and shear stresses. The SW-DVIB model independently reproduced this behavior through physics-based simulations, identifying a critical angle corresponding to the minimum failure load. Moreover, the model captured the transition from Mode I tensile fracture at low crack inclinations to Mode II shear-dominated fracture at higher angles, demonstrating its capability to simulate complex crack coalescence without predefined damage paths. Earlier plastic models, such as those developed by [PERSON] et al. [48] and [PERSON] et al. [49], lacked angular interactions, limiting their ability to represent anisotropic deformation and ductile flow, which are both essential for an accurate tunnel stability assessment. In contrast, the SW-DVIB model incorporates softening behavior through decoupling of the plasticity (two-body) and lateral stiffness (three-body) components. This allows for greater flexibility in the post-yield response tuning. Unlike the PSW-DVIB model [62], which assumes ideal plasticity, the present model allows for independent manipulation of the normal and shear components, thus reflecting a more accurate representation of the interactions governing brittle fracture. Traditional continuum approaches, such as the Finite Element Method (FEM) and Finite Difference Method (FDM), typically require external remeshing or damage-tracking algorithms to simulate the fracture evolution. Similarly, lattice-based methods often impose fixed Poisson's ratios and depend on pre-defined failure criteria [51; 52; 37; 53]. The SW-DVIB model mitigates these limitations by enabling both plasticity and fracture to evolve naturally through interparticle interactions.