# A Three-Dimensional Elastoplastic Constitutive Model for Geomaterials
## 7 Numerical Implementation  
### Elastoplastic Constitutive Relations' General Manifestation  
The current elastic zone's limit is set by the yield surface in the stress space. A stress point within this yield surface is in an elastic state and only has elastic properties. The point's stress state on the yield surface is in a plastic state with elastic or elastoplastic characteristics. A loading process that is determined by the loading criterion can be obtained if the stress state tends to move away from the yield surface.  
The loading condition in the proposed constitutive model is established using the loading criterion put out by [PERSON] [50] and [PERSON] et al. [51]. The stress tensor is determined using the generalized Hooke's law as follows:  
\[\sigma_{ij}=D^{e}_{ijkl}\Big{(}\epsilon_{kl}-\epsilon^{p}_{kl}\Big{)} \tag{21}\]  
where \(\sigma_{ij}\) represents the stress tensor; \(\epsilon_{kl}\) represents the strain tensor; \(\epsilon^{p}_{kl}\) represents the plastic strain tensor; and \(D^{e}_{ijkl}\) is the fourth-order elastic material tensor, which has the form [52]:  
\[D^{e}_{ijkl}=\frac{E}{2(1+v)}\Big{(}\delta_{ik}\delta_{jl}+\delta_{il}\delta_{ jk}\Big{)}+\frac{vE}{(1+v)(1-2v)}\delta_{ij}\delta_{kl} \tag{22}\]  
where \(v\) is [PERSON]'s ratio, \(E\) is [PERSON]'s modulus, and \(\delta_{ij}\) is the Kronecker symbol.  
One may write the yield function \(f\) as follows:  
\[f\Big{(}\sigma_{ij},H\Big{)}=F\Big{(}\epsilon_{kl},\epsilon^{p}_{kl},H\Big{)} \tag{23}\]  
where \(H\) is hardening parameter. Taking the partial differential of Equation (21) gives the following:  
\[\frac{\partial F}{\partial\epsilon_{ij}}=\frac{\partial f}{\partial\sigma_{kl }}D^{e}_{ijkl} \tag{24}\]  
Equation (25), which provides the loading criteria, provides their physical meanings (which are shown in Figure 9):  
\[\frac{\partial F}{\partial\epsilon_{ij}}d\epsilon_{ij}=\begin{cases}\frac{ \partial f}{\partial\sigma_{ij}}D^{e}_{ijkl}d\epsilon_{kl}>0&\text{Loading}\\ \frac{\partial f}{\partial\sigma_{ij}}D^{e}_{ijkl}d\epsilon_{kl}=0&\text{Neutral loading}\\ \frac{\partial f}{\partial\sigma_{ij}}D^{e}_{ijkl}d\epsilon_{kl}<0&\text{ Unloading}\end{cases} \tag{25}\]  
Figure 9: Loading criteria.  
According to [PERSON]'s law, the relationship between the elastic strain increases and the stress increment during the pure elastic deformation process is as follows:  
\[d\sigma_{ij}=D_{ijkl}^{e}d\epsilon_{kl}^{e} \tag{26}\]  
The whole strain increment can be split into elastic and plastic strain increments, according to the elastic-plastic theory:  
\[d\epsilon_{ij}=de_{ij}^{e}+de_{ij}^{p} \tag{27}\]  
where \(de_{ij}^{p}\) stands for an increment in the plastic strain \(\epsilon_{ij}^{p}\), and can be computed through the associated flow rule with the normality condition, as  
\[de_{ij}^{p}=\lambda\frac{\partial f}{\partial\sigma_{ij}} \tag{28}\]  
with \(\lambda\) as a scalar representing the plastic strain increment's magnitude, which is called the plastic scalar factor.  
Combining Equations (26) and (27), the formula for the stress-strain increment relationship is as follows:  
\[d\sigma_{ij}=D_{ijkl}^{e}\Big{(}d\epsilon_{ij}-d\epsilon_{ij}^{p}\Big{)} \tag{29}\]