# Study on Elastoplastic Damage Coupling of Reservoir Mudstone Considering Permeability Change
## 3 Elastoplastic Damage Coupling Mechanical Model
### Construction of the Mechanical Model  
Based on the existing research [25], the elastic free energy \(\varphi^{e}\) and plastic free energy \(\varphi^{p}\) of Helmholtz are, respectively:  
\[\varphi^{e}=\frac{1}{2}\epsilon^{e}:\mathbf{D}(\omega):\epsilon^{e} \tag{17}\]  
\[\varphi^{p}=(1-\omega)A\bigg{[}\gamma_{p}-b\Big{(}1-\alpha_{p}^{0}\Big{)}ln \frac{\gamma_{p}+b}{b}\bigg{]} \tag{18}\]  
where \(\mathbf{D}(\omega)\) is the fourth-order elastic stiffness tensor after material damage, and the bulk modulus \(K\) and shear modulus \(G\) can be expressed as the function of the damage variable \(\omega\)[26]:  
\[K(\omega)=(1-\omega)K_{0}\,\ G(\omega)=(1-\omega)G_{0} \tag{19}\]  
where \(K_{0}\) and \(G_{0}\) are the volume modulus and shear modulus of the material without damage, respectively. The damage driving force is obtained:  
\[Y=-\frac{\partial\varphi}{\partial\omega}=Y^{e}+Y^{p} \tag{20}\]  
\[Y^{e}=-\frac{1}{2}\epsilon^{e}:\mathbf{D}(\omega):\epsilon^{e} \tag{21}\]  
\[Y^{p}=A\bigg{[}\gamma_{p}-b\Big{(}1-\alpha_{p}^{0}\Big{)}ln\frac{\gamma_{p}+b }{b}\bigg{]} \tag{22}\]  
where \(\mathbf{D}^{\prime}(\omega)=\partial\mathbf{D}(\omega)/\partial\omega\). Based on the damage evolution of mudstone, the damage evolution equation of mudstone is defined as [27]:  
\[f^{\omega}=Y^{e}+\frac{Y^{p}}{Y^{p}+1}-\Big{(}m\omega+Y^{0}\Big{)}=0 \tag{23}\]where \(m\) is the damage parameter, which is related to the damage evolution, and \(Y^{0}\) is the damage threshold, which is related to the damage starting time. This paper takes \(Y^{0}=0\), which means deformation and damage simultaneously occur. When \(f^{\omega}<0\), there is no damage. When \(f^{\omega}=0\), damage occurs, and the damage evolution satisfies the consistency condition:  
\[\dot{f}^{\omega}=\frac{\partial Y^{e}}{\partial\epsilon^{e}}:\dot{\epsilon}^{e }+\frac{\partial f^{\omega}}{\partial Y^{p}}\frac{\partial Y^{p}}{\partial \gamma_{p}}\dot{\gamma_{p}}+\frac{\partial f^{\omega}}{\partial\omega}\dot{ \omega}=0 \tag{24}\]  
The average stress \(p\) and deviatoric stress \(q\) are used to describe the limit yield state. The failure of mudstone conforms to non-linear yield characteristics (Figure 6), and the non-linear plastic yield function corresponding to mudstone is obtained:  
\[ \begin{array}{c}f^{p}=q^{2}+A\alpha_{p}(p-C_{0})p_{0}=0\\ p=\frac{\alpha_{p}}{2},\;q=\sqrt{3J_{2}}\\ J_{2}=\frac{1}{2}S_{ij}S_{ij}\;,\;S_{ij}=\sigma_{ij}-\frac{\alpha_{p}}{3} \delta_{ij}\end{array}  \tag{25}\]  
where \(A\) is the internal friction coefficient, \(C_{0}\) is the cohesive coefficient, and \(\alpha_{p}\) is the hardening function, which can be expressed as a power function monotonically increasing with the plastic internal variables \(\gamma_{p}\):  
\[ \begin{array}{c}\alpha_{p}=(1-\omega)\left[a_{p}^{0}-\left(a_{p}^{0}-1 \right)\frac{\gamma_{p}}{g+\gamma_{p}}\right]\\ \gamma_{p}=\sqrt{\frac{2}{8}\epsilon_{ij}^{p}\epsilon_{ij}^{p}}\;,\;\;e_{ij}^ {p}=\epsilon_{ij}^{p}-\frac{\epsilon_{p}^{p}}{3}\delta_{ij}\end{array}  \tag{26}\]  
where \(\alpha_{p}\) is the initial value of the hardening function, and \(b\) is the hardening rate parameter. The variation of the yield surface with plastic deformation and damage evolution can be described by \(\alpha_{p}\): before the peak stress, the yield surface expands with the increase in plastic deformation, showing hardening characteristics; after the peak stress, the yield surface decreases with the increase in damage, showing softening characteristics; and when the damage increases to the maximum, the mudstone reaches the residual strength. In the loading process, when \(f^{p}<0\), it is only elastic deformation.