## 1 Introduction  
With the increasing global energy demands and rapid infrastructure development, the investigation of crack propagation and failure mechanisms in rock engineering has emerged as a pivotal research focus in geomechanics, hydraulic engineering, and energy development [1, 2, 3]. In practical engineering scenarios such as hydraulic fracturing, tunnel excavation, and the construction of underground energy storage caverns, the processes of crack initiation, propagation, and penetration in rock under coupled stress-seepage fields directly influence the safety and operational reliability of engineering structures [4]. Traditional fracture mechanics approaches--including linear elastic fracture mechanics (LEFM), stress intensity factor (SIF) methods, and extended finite element methods (XFEM)--exhibit inherent limitations in the efficient simulation of dynamic crack evolution and multi-fieldcoupling effects. These limitations principally stem from their reliance on predefined crack paths and the computational complexity associated with handling intricate boundary conditions [5]. As a continuum-based numerical technique, the phase field method (PFM) circumvents explicit crack surface tracking through the modeling of crack evolution as a continuous distribution of phase field variables via a diffuse crack model, demonstrating significant advantages in addressing complex crack propagation under multi-field coupling [6; 7; 8]. Originating from the energy variational principle proposed by [PERSON] et al. [9], the PFM unifies crack initiation, propagation, and branching within the framework of a continuous medium framework by formulating an energy functional that incorporates both crack surface energy and elastic strain energy. Recent advancements in its application to rock mechanics include the following: [PERSON] et al. [10] integrated thermodynamics and damage mechanics to develop a phase field model for brittle and quasi-brittle materials, resolving the length-scale dependency issue in traditional models; [PERSON] et al. [11] implemented a brittle material phase field fracture model in the COMSOL Multiphysics platform, clarifying through sensitivity analysis the influences of Young's modulus and critical energy release rate on crack trajectories. However, these elastic-based models neglect the critical role of plastic deformation in compressive-shear failure, leading to overly brittle predictions and an inability to capture residual strength caused by plastic slip [12]. Parallel research efforts in plastic damage modeling have focused on characterizing irreversible deformation mechanisms. Classical plasticity theories, such as the Drucker-Prager and Mohr-Coulomb criteria, have been widely adopted to simulate shear band localization and post-yield softening in rocks [13]. Although [PERSON] et al. [14] proposed an elastoplastic damage model for rocks, their approach relied on cohesive zone methods requiring predefined crack paths. While these studies have advanced the understanding of plastic deformation, they lack a unified framework to couple plasticity with fracture evolution, particularly in phase field formulations. A critical gap exists in the theoretical integration of plasticity theory with phase field fracture models. Existing compressive-shear phase field models predominantly rely on elastic damage theory, which only characterizes failure through strain energy decomposition and residual strength features. For instance, [PERSON] et al. [15] developed a hydromechanical coupled phase field model but omitted plastic strain energy dissipation. Similarly, [PERSON] et al. [16] simulated hydraulic fracturing in bedded shale using PFM yet treated the rock as purely elastic. These approaches fail to quantitatively assess the influence of plastic deformation on crack driving forces and propagation paths, significantly limiting their applicability to real-world scenarios where plastic slip dominates energy dissipation.