# Elastoplastic Constitutive Model for Energy Dissipation and Crack Evolution in Rocks
## 5 Validation of the Theoretical Model
### Model Development and Numerical Simulation  
Based on the algorithm implementation of the theoretical elastoplastic constitutive model, as outlined in Table 1, a secondary development of the argillaceous silsttone model, as presented in Table 2, was conducted using C++ in the Microsoft Visual Studio 2010 development environment. The developed model was subsequently integrated into FLAC(30) 6.0 for a numerical simulation. To ensure the internal variables of the developed model are strictly consistent with the theoretical calculations [9], a single cubic element is established for the numerical simulations of the conventional triaxial compression at confining pressures of 1, 2, 3, 4, and 5 MPa, as shown in Figure 6. During the confining pressure loading phase, the element is subjected to equal confining pressures \(\sigma_{3}\) in the X, Y, and Z directions until stabilization. During the deviatoric loading phase, the bottom face is fixed, while the top face is subjected \begin{table}
\begin{tabular}{c c} \hline
**Model Parameter** & **Parameter Equation** \\ \hline Average elastic modulus \(\overline{E}\) & 6.51 GPa \\ Average Poisson’s ratio \(\overline{\mu}\) & 0.244 \\ Yield function \(f\) & Equation (3) \\ Cohesion \(c\) & Equation (44) \\ Internal friction angle \(\varphi\) & Equation (45) \\ Plastic shear strain \(\gamma_{p}\) & Equation (8) \\ Potential function \(g\) & Equation (9) \\ Dilatancy angle \(\psi\) & Equation (46) \\ Elastic strain energy stored during the deviatoric loading phase \(U^{e}-U_{0}\) & Equations (20) and (22) \\ Axial and circumferential work \(U_{1}+U_{3}\) & Equations (23) and (24) \\ Dissipated energy \(U^{d}\) & Equation (25) \\ Volumetric strain \(\varepsilon_{V}\) & Equation (27) \\ Elastic volumetric strain \(\varepsilon_{V}^{e}\) & Equation (28) \\ \hline \end{tabular}
\end{table} 
Table 2: Argillaceous silststone model. Figure 5: Relationship between dilatancy angle and plastic strain. to a downward rate of 0.002 mm/s until failure. The entire loading process records the following parameters: \(\gamma_{p}\), \(c\), \(\varphi\), \(\psi\), \(\sigma_{1}\), \(\sigma_{3}\), \(\epsilon_{1}\), \(\epsilon_{3}\), \(U_{1}+U_{3}\), \(U^{\epsilon}-U_{0}\), \(U^{d}\), \(\varepsilon_{V}\), \(\varepsilon_{V}^{\epsilon}\), and \(\varepsilon_{V}^{\epsilon}\). The processed data are illustrated in Figures 7-10. Figure 7 demonstrates that the simulated degradation of \(c\), \(\varphi\), and \(\psi\) with \(\gamma_{p}\) is consistent with the theoretical values of Equations (44)-(46). This consistency is because the degradation of parameters \(c\), \(\varphi\), and \(\psi\) in the developed model is programmed based on \(c(\gamma_{p})\), \(\varphi(\gamma_{p})\), and \(\psi(\gamma_{p})\), indirectly confirming the successful integration of these parameters into the model. Figure 8 shows that the simulated stress-strain curves at various confining pressures closely match the experimental data. This concordance is due to the variable \(\gamma_{p}\) in Equation (8), which incorporates deviatoric stress (\(\sigma_{1}-\sigma_{3}\)) and effectively captures both the plastic strain extent and the confining pressure effect. The dependent variables \(c(\gamma_{p})\), \(\varphi(\gamma_{p})\), and \(\psi(\gamma_{p})\) adjust the entire yielding process of the rock. Although there are slight deviations between the simulated and test data, they objectively represent the  
Figure 6: Loading diagram of the cubic element. Figure 7: Simulated relationships of strength parameters and dilatancy angles: (a) \(c,\varphi-\gamma_{p}\); (b) \(\psi-\gamma_{p}\). overall yield behavior of the rock under varying levels of confining pressure. This validates that the theoretical elastoplastic constitutive model effectively characterizes its nonlinear mechanical behavior. As depicted in Figures 9 and 10a, during the elastic phase (\(\gamma_{p}=0\)), \(U^{d}=0\), \(U_{1}+U_{3}\) overlaps entirely with \(U^{e}-U_{0}\). During yielding, \(U^{d}=(U_{1}+U_{3})-(U^{e}-U_{0})\).