# Coupled Elastic-Plastic Damage Modeling of Rock Based on Irreversible Thermodynamics
## 2 Basic Framework of the Elastoplastic Damage Coupled Model
### Elastoplastic Model  
The description of the material's plastic behavior involves combining the yield criterion, flow rule, and hardening function [52; 69]. Among these, the Drucker-Prager yield criterion is widely used for shale [14; 27; 87; 88; 89; 90; 91]. However, when dealing with strain hardening and strain softening behavior, traditional yield criteria are no longer suitable. They must be improved and optimized to allow for changes in the yield surface during dilation and contraction. Therefore, the modified yield criterion used in this article is:  
\[F^{p}(\sigma,\gamma^{p},\omega)=q+\alpha_{1}p-\kappa\gamma^{p\omega}=0 \tag{38}\]  
where \(q\) is the equivalent stress, expressed as \(q=\sqrt{1.55:S}\); \(p\) is the average stress, expressed as \(p=tr(\sigma)/3\); \(S\) is the deviatoric stress tensor; \(\alpha_{1}\) and \(\kappa\) are the Drucker-Prager strength parameters, representing internal friction and cohesion, respectively; \(\chi^{p\omega}\) is a function of \(\gamma^{p},f(\theta)\), and \(\omega\), reflecting the impact of plastic hardening, damage softening, and shale bedding angle on the evolution of the yield surface, expressed as:  
\[\chi^{p\omega}=f(\theta)\Big{[}h_{1}-(h_{1}-h_{0})\mathrm{e}^{(-p\gamma^{p})} \Big{]}(1-\omega)^{a} \tag{39}\]  
where \(h_{0}\) and \(h_{1}\) are the initial and final yield thresholds, respectively; \(b\) is the plastic hardening rate, reflecting the rate of hardening. The smaller the value of \(b\), the faster the hardening rate.  
Previous research has shown that post-peak strength degradation is related to confinement conditions; thus, it is necessary to introduce the parameter \(a\) to describe the influence of confining pressure on damage changes during triaxial experiments [26; 59; 83]. Moreover, shale and many other rocks exhibit significant directional dependence and anisotropy in mechanical behavior due to their naturally occurring bedding and joint planes [70; 92; 93; 94; 95]. To accurately describe this anisotropy caused by differences in bedding plane angles, this study incorporates a function \(f(\theta)\) related to the bedding plane angle into the yield surface evolution equation (Equation (31)). The magnitude of \(f(\theta)\) depends on the angle between the stress loading direction and the bedding plane (Figure 2). The specific expressions for parameters \(a\) and \(f(\theta)\) will be provided during the validation phase.  
Rock materials exhibit apparent volumetric contraction during triaxial loading [96], and traditionally associated flow rules overlook the energy dissipation and volumetric changes during plastic yielding. Therefore, this article adopts a non-associated flow rule, obtained from:  
\[G^{p}(\sigma,\gamma^{p},\omega)=q+\alpha_{2}p \tag{40}\]where \(\alpha_{2}\) is a parameter related to shear volume change.  
In the elastoplastic loading process without damage evolution, through the non-associated flow rule and consistency condition, it can be obtained:  
\[d\lambda^{p}=\frac{\partial F^{p}}{\partial\sigma}:D:d\varepsilon/H^{p} \tag{41}\]  
where \(H^{p}\) is the plastic modulus, expressed as:  
\[H^{p}=\frac{\partial F^{p}}{\partial\sigma}:D:\frac{\partial G^{p}}{\partial \sigma}-\frac{\partial F^{p}}{\partial\gamma^{p}}\sqrt{\frac{2}{3}\frac{ \partial G^{p}}{\partial S}:\frac{\partial G^{p}}{\partial S}} \tag{42}\]  
The incremental form of the stress-strain relationship is:  
\[d\sigma=D^{ep}:d\varepsilon \tag{43}\]  
\[D^{ep}=D-\frac{1}{H^{p}}(D:\frac{\partial G^{p}}{\partial\sigma})\otimes( \frac{\partial F^{p}}{\partial\sigma}:D) \tag{44}\]