# A Three-Dimensional Elastoplastic Constitutive Model for Geomaterials
## 7 Numerical Implementation
### The Mixed Complementarity Problem and its GSPC Algorithm
\begin{array}{l}\text{Step0}:\beta_{I}=\beta_{E}=1;k=0;\\ \text{Step1}:\overline{\lambda}=\max[\lambda-\beta_{I}f_{1}(\lambda,\sigma),0]; \\ \overline{\sigma}=\sigma-\beta_{E}f_{E}(\overline{\lambda},\sigma)\\ \text{If}\left\|\overline{\lambda}-\lambda\right\|_{\infty}\leq\varepsilon_{ \lambda}\left\|\overline{\lambda}\right\|_{\infty}\text{ and }\left\|\overline{\sigma}-\sigma\right\|_{\infty}\leq \varepsilon_{\sigma}\left\|\overline{\sigma}\right\|_{\infty}\text{, then }\lambda= \overline{\lambda};\sigma=\overline{\sigma};\text{break;}\\ r_{\lambda}=\frac{\beta_{I}\left\|f_{1}(\lambda,\sigma)-f_{1}(\lambda, \sigma)\right\|_{2}}{\left\|\lambda-\lambda\right\|_{2}}\\ \text{while }r_{\lambda}>\
u\\ \beta_{I}=\frac{\beta_{I}}{2}\beta_{I}\min\left(1,\frac{1}{r_{\lambda}}\right) ;\overline{\lambda}=\max[\lambda-\beta_{I}f_{1}(\lambda,\sigma),0];\\ r_{\lambda}=\frac{\beta_{I}\left\|f_{1}(\lambda,\sigma)-f_{1}(\lambda, \sigma)\right\|_{2}}{\left\|\lambda-\lambda\right\|_{2}};\\ \text{end(while)}\\ r_{\sigma}=\frac{\beta_{E}\left\|f_{E}(\overline{\lambda},\sigma)-f_{E}( \overline{\lambda},\overline{\sigma})\right\|_{2}}{\left\|\sigma-\overline{ \sigma}\right\|_{2}};\\ \text{while }r_{\sigma}>\
u\\ \beta_{E}=\frac{2}{3}\beta_{E}\min\left(1,\frac{1}{r_{\sigma}}\right); \overline{\sigma}=\sigma-\beta_{E}f_{E}(\overline{\lambda},\sigma);\\ r_{\sigma}=\frac{\beta_{E}\left\|f_{E}(\overline{\lambda},\sigma)-f_{E}( \overline{\lambda},\overline{\sigma})\right\|_{2}}{\left\|\sigma-\overline{ \sigma}\right\|_{2}};\\ \text{end(while);}\\ d_{\lambda}(\lambda,\overline{\lambda})=(\lambda-\overline{\lambda})-\beta_{I} \left[f_{1}(\lambda,\sigma)-f_{1}(\overline{\lambda},\sigma)\right];\\ \alpha=\frac{(\lambda-\overline{\lambda})^{T}d_{\lambda}(\lambda,\overline{ \lambda})}{\left\|d_{\lambda}(\lambda,\overline{\lambda})\right\|_{2}^{2}}; \lambda=\lambda-\gamma ad_{\lambda}(\lambda,\overline{\lambda});\\ \text{If }r_{\lambda}\leq\mu\text{ then }\beta_{I}=1.5\beta_{I};\\ d_{\sigma}(\sigma,\overline{\sigma})=(\sigma-\overline{\sigma})-\beta_{E} \left[f_{E}(\overline{\lambda},\sigma)-f_{E}(\overline{\lambda},\overline{ \sigma})\right];\\ \alpha=\frac{(\sigma-\overline{\sigma})^{T}d_{\sigma}(\sigma,\overline{ \sigma})}{\left\|d_{\sigma}(\sigma,\overline{\sigma})\right\|_{2}^{2}}; \sigma=\sigma-\gamma ad_{\sigma}(\sigma,\overline{\sigma});\\ \text{If }r_{\sigma}\leq\mu\text{ then }\beta_{E}=1.5\beta_{E};\\ \text{Step2}:k=k+1;\text{go to Step 1}\\ \end{array} \]  
The GSPC algorithm has the following characteristics: (1) The GSPC converges for any size of the strain increment \(\Delta\varepsilon\), and the implicit return mapping algorithm calculates the stress \(\sigma\) and the plastic multiplier \(\lambda\) using the Newton technique. When the strain increment is small enough, the convergence can be ensured because the Newton technique is locally convergent. (2) The GSPC algorithm's computational efficiency is higher than the return mapping algorithm because it is simple to program and the second derivative of the potential function does not need computing. The return mapping algorithm needs to first determine the stress path and yield surface corner points, as well as whether those yield surfaces are active, before forming the plastic potential function's Jacobian matrix. The Hessian matrix of the plastic potential function must be determined to derive the Jacobian matrix. Usually, the return mapping algorithm requires a lot of trial and error for this method to provide accurate results. Except for \(\phi_{b}\) being replaced with the expansion angle \(\psi_{b}\), the expression of \(f\) and the plastic potential function \(g\) are consistent. When using the GSPC integral algorithm to update the stress and calculate the plastic multiplier, the first-order derivative of the yield function \(F_{GMC}\) to stress is required (refer to Appendix A). ABAQUS provides users with many user subroutine interfaces, thus allowing users to write the required user subroutines according to the corresponding problems encountered in the actual project.