## 2 Visco-Elastic-Plastic Constitutive Model for Squeezing Rock
### Strain Integration Scheme of the Rheological Model
A complete implicit method with the characteristics of unconditional stability is used to solve this model. Using Taylor series expansion on \(\left\{\hat{e}_{vp}^{t+1}\right\}\) in Equation (30), we achieve:  
\[\left\{\Delta\varepsilon_{vp}\right\}=\left\{\hat{e}_{vp}^{n}\right\}\Delta t_ {n}+[C^{n}]\Delta\sigma^{n} \tag{26}\]  
\[[C^{n}]=\Delta t_{n}[H^{n}],\;[H^{n}]=\frac{\partial\hat{e}_{vp}}{\partial \sigma^{n}}=\frac{1}{\eta_{3}}\bigg{\{}\phi\cdot\frac{\partial a^{T}}{ \partial\sigma}+\frac{d\phi}{dF}a\cdot a^{T}\bigg{\}}. \tag{27}\]  
The Hermitian matrix derived from the Drucker-Prager yield criterion:  
\[[H]=\frac{\phi}{\eta_{3}}\cdot\left(\frac{1}{{{{{{{{{{{{{{{{{{{{{{{ }}}}}}}}}}}}}}}}}}}^{2}}{ \Omega_{1}}-\frac{1}{{{{{{{{{{{{{{{{}}}}}}}}}}}}}^{3}}}\Omega_{2}}\right)+ \frac{1}{\eta_{3}}\frac{d\phi}{dF}\cdot\left(Q_{3}+\frac{1}{\sqrt{{{{{{{{ {{{{\left({{{{{{{\left}}}}}}}}}}}}}}}^{2}}}Q_{4}+\frac{1}{{{{{{{{{{{{ {\left({{{{{}}}}}}}}}}}{}}}}}^{2}}}Q_{2}}\right)\right) \tag{28}\]  
where  
\[Q_{1}=\left[ \begin{array}{cccccc}\frac{4}{3}&-\frac{2}{3}&-\frac{2}{3}&0&0&0 \\ &\frac{4}{3}&-\frac{2}{3}&0&0&0\\ &\frac{4}{3}&0&0&0\\ &&4&0&0\\ \text{\it sysm}&&&4&0\\ &\end{array} \right] \tag{29}\]  
\[Q_{2}=\left[ \begin{array}{cccccc}4(\sigma_{x}^{\prime})^{2}&4\sigma_{x}^{ \prime}\sigma_{y}^{\prime}&4\sigma_{x}^{\prime}\sigma_{z}^{\prime}&8\tau_{xy} \sigma_{x}^{\prime}&8\sigma_{x}^{\prime}&8\tau_{xz}\sigma_{x}^{\prime}\\ &\left(\sigma_{y}^{\prime}\right)^{2}&\sigma_{y}^{\prime}\sigma_{z}^{\prime} &8\tau_{xy}\sigma_{y}^{\prime}&8\tau_{yz}\sigma_{y}^{\prime}&8\tau_{xz}\sigma _{y}^{\prime}\\ &&(\sigma_{z}^{\prime})^{2}&8\tau_{xy}\sigma_{z}^{\prime}&8\tau_{yz}\sigma_{z}^ {\prime}&8\sigma_{z}^{\prime}\\ &&&16\tau_{xy}^{2}&8\tau_{yz}\sigma_{x}\tau_{xz}\\ \text{\it sysm}&&&16\tau_{yz}^{2}&8\tau_{yz}\tau_{xz}\\ &&&&16\tau_{zz}^{2}\end{array} \right] \tag{30}\]  
\[Q_{3}=\left[ \begin{array}{cccccc}\alpha^{2}&\alpha^{2}&\alpha^{2}&0&0&0\\ &\alpha^{2}&\alpha^{2}&0&0&0\\ &&\alpha^{2}&0&0&0\\ &&0&0&0\\ \text{\it sysm}&&&0&0\\ &\end{array} \right] \tag{31}\]  
\[Q_{4}=\left[ \begin{array}{cccccc}4\alpha\sigma_{x}^{\prime}&\alpha\left( \sigma_{x}^{\prime}+\sigma_{y}^{\prime}\right)&\alpha(\sigma_{x}^{\prime}+ \sigma_{z}^{\prime})&4\alpha\tau_{xy}&4\alpha\tau_{yz}&4\alpha\tau_{xz}\\ &4\alpha\sigma_{y}^{\prime}&\alpha\left(\sigma_{y}^{\prime}+\sigma_{z}^{ \prime}\right)&4\alpha\tau_{xy}&4\alpha\tau_{yz}&4\alpha\tau_{xz}\\ &&&4\alpha\sigma_{z}^{\prime}&4\alpha\tau_{xy}&4\alpha\tau_{yz}&4\alpha\tau_{xz }\\ &&&0&0&0\\ \text{\it sysm}&&&0&0\\ \text{\it sysm}&&&0&0\end{array} \right] \tag{32}\]  
According to the visco-elastic-plastic incremental constitutive relationship of the rock materials, the stress and strain correlation can be expressed as  
\[\left\{\Delta\sigma\right\}_{n}=\left[\hat{D}\right]\left(\left\{\Delta \varepsilon\right\}_{n}-\left\{\Delta\varepsilon_{vc}\right\}_{n}-\left\{ \hat{\varepsilon}_{cp}\right\}_{n}\cdot\Delta t_{n}\right) \tag{33}\]  
where the visco-elastic-plastic tangent modulus is \(\left[\hat{D}\right]=\left([D]^{-1}+[C]\right)^{-1}\). When the stress level is low (\(\phi(F)\leq 0\)), \(\left\{\hat{\varepsilon}_{vp}\right\}\) is \(0\) and results in \(\left[\hat{D}\right]\) being equal to \([D]\). According to the principle of virtual work, the incremental balance equation of the system can be derived as  
\[\sum_{n}\int_{v_{n}}[B]_{n}^{T}\{\Delta\sigma\}_{n}dv-\left\{\Delta R\right\}_{ n}=0 \tag{34}\]where \(\left[B\right]_{n}\) and \(\left\{\Delta R\right\}_{n}\) are the geometric and the external load increment matrix at \(\Delta t_{n}\).