## 2 Formulation and Theoretical Derivation of the Elastoplastic Phase Field Fracture Model
### Elastic Strain Energy  
For quasi-brittle rock materials, under the small-strain assumption, the total strain \(\varepsilon\) can be decomposed into elastic strain \(\varepsilon^{e}\) and plastic strain \(\varepsilon^{p}\)[11]  
\[\varepsilon=\varepsilon^{e}+\varepsilon^{p} \tag{5}\]  
During the elastic deformation stage, to address tension-compression asymmetry in strain energy-driven damage evolution, the elastic strain tensor \(\varepsilon^{e}\) is decomposed into tensile \(\varepsilon^{e}_{+}\) and compressive \(\varepsilon^{e}_{-}\) components:  
\[\varepsilon^{e}_{\pm}=\sum_{a=1}^{d}\langle\varepsilon^{e}_{a}\rangle_{\pm} \,\mathbf{n}_{a}\otimes\mathbf{n}_{b} \tag{6}\]  
In this formula, \(\varepsilon^{e}_{a}\) denotes the principal elastic strains; \(\mathbf{n}_{a}\) and \(\mathbf{n}_{b}\) represent the principal strain directions, respectively; and \(\langle\varepsilon^{e}_{a}\rangle_{\pm}\) is defined via the Macaulay operator \(\langle\varepsilon^{e}_{a}\rangle_{\pm}=(\varepsilon^{e}_{a}\pm|\varepsilon^{ e}_{a}|)/2\). In rock deformation and failure processes, material damage is predominantly attributed to local elastic tensile strains. Thus, the total elastic strain energy density during fracture can be expressed as [24]  
\[\psi(\varepsilon)=g(\phi)\psi^{e}_{+}+\psi^{e}_{-} \tag{7}\]  
In this formula, \(\psi^{e}_{+}\) and \(\psi^{e}_{-}\) correspond to the elastic strain energy densities under tensile and compressive stress states, respectively. At complete fracture \(\phi=1\), tensile elastic strain energy dissipates progressively. To model this behavior, a phase field damage degradation function is introduced:  
\[g(\phi)=(1-\delta)(1-\phi)^{2}+\delta \tag{8}\]  
In this formula, \(\delta\) is a small regularization parameter to avoid numerical model singularity when \(\phi=1\). In engineering practice, rocks commonly experience compressive-shear stress states, exhibiting distinctive compressive-shear failure characteristics. Conventional phase field theories predominantly address tensile fracture and struggle to accurately predict compressive-shear crack propagation [25; 26; 27]. Theoretical analyses demonstrate that strain energy decomposition directly governs crack trajectories and propagation modes. To accurately model rock failure under prevalent compressive-shear conditions in underground engineering, the total strain energy is partitioned into three physically interpretable components: tensile strain energy--drives crack nucleation at stress concentrations; tensile-shear strain energy--governs mixed-mode crack growth under combined stresses; compressive-shear strain energy--regulates shear band formation and stable crack propagation. This decomposition accounts for the asymmetric tension-compression response of rocks: tensile stresses initiate cracks at notch tips, while compressive-shear stresses promote shear localization and post-peak strength evolution. To characterize post-failure residual strength, a compressive-shear residual strength parameter \(\eta\) is incorporated, modifying the elastic strain energy density function in Equation (7) to  
\[\psi(\varepsilon)=g(\phi)(\psi_{+}^{e,I}+\psi_{+}^{e,\Pi})+\{g(\phi)+\eta[1-H( tr(\varepsilon))][1-g(\phi)]\}\psi_{-}^{e} \tag{9}\]  
In this formula, \(\psi_{+}^{e,\Pi}\) and \(\psi_{+}^{e,\Pi}\) represent the volumetric and deviatoric tensile strain energy densities, respectively; \(\psi_{-}^{e}\) denotes the compressive-shear elastic strain energy density; \(\eta\) is the residual strength parameter; and \(H(x)\) is the Heaviside step function  
\[H(x)= \begin{cases}0,&x<0\\ 1,&x\geq 0\end{cases}  \tag{10}\]  
The mathematical expressions for strain energy components are  
\[\left\{ \begin{array}{l}\psi_{+}^{e,\Pi}=\frac{1}{2}\lambda\langle\text{tr }(\varepsilon)\rangle_{+}^{2}\\ \psi_{+}^{e,\Pi}=\mu\text{tr}\left(\langle\varepsilon\rangle_{+}^{2}\right) \\ \psi_{-}^{e}=\frac{1}{2\mu}\tau^{2}\end{array} \right.