# Elastoplastic Constitutive Model for Energy Dissipation and Crack Evolution in Rocks
## 4 Construction and Algorithm Implementation of the Theoretical Model  
This section constructs the corresponding theoretical elastoplastic constitutive model based on the analyses presented in Sections 2 and 3, as shown in Table 1. Additionally, the three-dimensional finite difference scheme of the theoretical model is derived, facilitating the implementation of the numerical algorithm within FLAC\({}^{\text{3D}}\). The three-dimensional finite difference format of the theoretical model is as follows: Firstly, the next step stress \(\sigma_{n}^{N}\) is obtained from the current step stress \(\sigma_{n}^{c}\) and strain increment \(\Delta\varepsilon_{i}\). The elastic guess stress tensor \(\sigma_{n}^{I}\) can be expressed as follows [45]:  
\[\sigma_{n}^{I}=\sigma_{n}^{c}+\text{S}_{n}(\Delta\varepsilon_{i}),n=1,2,3,i=1, 2,3 \tag{30}\]  
The components of \(\text{S}_{n}(\Delta\varepsilon_{i})\) are as follows:  
\[ \begin{cases}\text{S}_{1}(\Delta\varepsilon_{1},\Delta\varepsilon_{2},\Delta \varepsilon_{3})=\delta_{1}\Delta\varepsilon_{1}+\delta_{2}(\Delta\varepsilon _{2}+\Delta\varepsilon_{3})\\ \text{S}_{2}(\Delta\varepsilon_{1},\Delta\varepsilon_{2},\Delta\varepsilon_{3 })=\delta_{1}\Delta\varepsilon_{2}+\delta_{2}(\Delta\varepsilon_{1}+\Delta \varepsilon_{3})\\ \text{S}_{3}(\Delta\varepsilon_{1},\Delta\varepsilon_{2},\Delta\varepsilon_{3 })=\delta_{1}\Delta\varepsilon_{3}+\delta_{2}(\Delta\varepsilon_{1}+\Delta \varepsilon_{2})\end{cases}  \tag{31}\]where \(\delta_{1}=E(1-\mu)/[(1-2\mu)(1+\mu)]\), \(\delta_{2}=E\mu/[(1-2\mu)(1+\mu)]\). The values of \(\delta_{1}\) and \(\delta_{2}\) are determined by Equation (2). Meanwhile, the plastic shear strain \(\gamma_{p}\) in its total form is given by Equation (8), and, numerically, \(\gamma_{p}\) is computed in an incremental form:  
\[\gamma_{p}=\gamma_{p}^{c}+\Delta\gamma_{p} \tag{32}\]  
\[\Delta\gamma_{p}=\sqrt{3}\Delta e^{p}\big{(}1+N_{\psi}\big{)}/\sqrt{1+N_{\psi }+N_{\psi}{}^{2}} \tag{33}\]  
\[\Delta e^{p}=\sqrt{\left[\left(\Delta e_{1}^{p}-\Delta e_{m}^{p}\right)^{2}+ \left(\Delta e_{m}^{p}\right)^{2}+\left(\Delta e_{3}^{p}-\Delta e_{m}^{p} \right)^{2}\right]/2} \tag{34}\]  
\[\Delta e_{m}^{p}=\left(\Delta e_{1}^{p}+\Delta e_{3}^{p}\right)/3 \tag{35}\]  
where \(\gamma_{p}^{c}\) and \(\Delta\gamma_{p}\) are the current step plastic shear strain and plastic shear strain increment, respectively. The total and incremental forms of \(\gamma_{p}\) are numerically equal. The internal friction angle \(\varphi\) and cohesion \(c\) within the yield function \(f\) are determined by \(c\big{(}\gamma_{p}\big{)}\) and \(\big{(}\gamma_{p}\big{)}\), as given in Equation (3). The dilatancy angle \(\psi\) in the potential function \(g\) is determined by \((\gamma_{p})\), as given in Equation (9). Secondly, substituting \(\sigma_{n}^{I}\) into Equation (3), if \(f\big{(}\sigma_{n}^{I}\big{)}<0\), the material is in the elastic stage and \(\sigma_{n}^{N}=\sigma_{n}^{I}\); if \(f\big{(}\sigma_{n}^{I}\big{)}\geq 0\), the material reaches the yield stage, and the strain increment \(\Delta\varepsilon_{i}\) splits into elastic and plastic parts:  
\[\Delta\varepsilon_{i}=\Delta e_{i}^{e}+\Delta e_{i}^{p} \tag{36}\]  
The stress increment for the current step \(\Delta\sigma_{n}\) can be expressed as follows:  
\[\Delta\sigma_{n}=\mathrm{S}_{n}(\Delta\varepsilon_{i}^{e}) \tag{37}\]  
The next step stress \(\sigma_{n}^{N}\) can be expressed as follows:  
\[\sigma_{n}^{N}=\sigma_{n}^{c}+\Delta\sigma_{n} \tag{38}\]  
By substituting Equations (10), (30), (36), and (37) into Equation (38), we obtain the following:  
\[\sigma_{n}^{N}=\sigma_{n}^{I}-\mathrm{S}_{n}\bigg{(}\lambda\frac{\partial g} {\partial\sigma_{i}}\bigg{)} \tag{39}\]