# High-Precision Elastoplastic Four-Node Quadrilateral Shell Element
## 4 Plastic Analysis  
The deformation rate tensor \(\dot{\epsilon}\) is decomposed into the elastic \(\dot{\epsilon}^{\text{e}}\) and plastic \(\dot{\epsilon}^{\text{P}}\) parts:  
\[\dot{\epsilon}=\dot{\epsilon}^{\text{e}}+\dot{\epsilon}^{\text{P}} \tag{61}\]  
Let \(\sigma_{ij}\) be the stress tensor, \(\alpha_{ij}\) be the back stress tensor, and \(s_{ij}\) be the deflection stress tensor after taking into account the back stress, i.e.,  
\[s_{ij}=\sigma_{ij}-\alpha_{ij}-\delta_{ij}(\sigma_{mm}-\alpha_{mm})/3 \tag{62}\]  
where \(\delta_{ij}\) is the Kronecker delta. For ideal elastoplastic or isotropic hardening, the back stress \(\alpha_{ij}=0\). In \(J_{2}\) flow theory, the equivalent stress \(\overline{\sigma}\) is defined as a function of the second invariant \(J_{2}\) of the tensor \(s_{ij}\):  
\[\overline{\sigma}=\sqrt{3J_{2}}=\sqrt{1.5s_{ij}s_{ij}} \tag{63}\]  
The yield condition is  
\[f=\sqrt{3J_{2}}-\sigma_{y}(\vec{\epsilon}^{\text{P}})=\overline{\sigma}- \sigma_{y}(\vec{\epsilon}^{\text{P}})=0 \tag{64}\]  
where \(\sigma_{y}(\vec{\epsilon}^{\text{P}})\) represents the yield stress of the material, which typically depends on the internal variable \(\vec{\epsilon}^{\text{P}}=\int\text{d}\vec{\epsilon}^{\text{P}}\). \(\text{d}\vec{\epsilon}^{\text{P}}\) denotes the increment in the equivalent plastic strain. \(f<0\) corresponds to elasticity and \(f=0\) represents cumulative equivalent plastic anisotropy. The ideal elastic-plastic models, isotropic hardening, kinematic hardening, and mixed hardening can be expressed in a unified form by assuming a linear relationship between \(\sigma_{y}(\vec{\epsilon}^{\text{P}})\) and \(\sigma_{y}(\vec{\epsilon}^{\text{P}})\) (i.e., linear strengthening). Additionally, a parameter \(\beta\) is introduced to facilitate this relationship. \[\sigma_{y}=\sigma_{y0}+\beta E_{\text{P}}\vec{\epsilon}^{\text{P}} \tag{65}\]  
\[\alpha_{ij}^{\text{VJ}}=\frac{2}{3}(1-\beta)E_{\text{P}}\dot{\epsilon}^{\text {P}}_{ij} \tag{66}\]  
where \(E_{\text{P}}=\text{d}\sigma_{y}/\text{d}\vec{\epsilon}^{\text{P}}=0\) is the ideal elastoplastic case; if \(E_{\text{P}}\
eq 0\), \(\beta=1\) corresponds to isotropic hardening; if \(E_{\text{P}}\
eq 0\), \(\beta=0\) corresponds to kinematic hardening; if \(E_{\text{P}}\
eq 0\), \(\beta\in(0,1)\) corresponds to mixed hardening. According to the correlated plastic flow law, the plastic part of the deformation rate tensor \(\dot{\epsilon}\) is  
\[\dot{\epsilon}^{\text{P}}_{ij}=\dot{\lambda}\frac{\partial f}{\partial\sigma_ {ij}}=\frac{\dot{\lambda}}{2}\sqrt{\frac{3}{J_{2}}}\frac{\partial J_{2}}{ \partial\sigma_{ij}}=\frac{3s_{ij}}{2\overline{\sigma}}\dot{\lambda} \tag{67}\]  
where \(\dot{\lambda}\) is the scalar plastic flow rate. The Von Mises equivalent plastic deformation rate \(\dot{\vec{\epsilon}}^{\text{P}}\) is defined as  
\[\dot{\vec{\epsilon}}^{\text{P}}=\left(\frac{2}{3}\dot{\epsilon}^{\text{P}}_{ij }\dot{\epsilon}^{\text{P}}_{ij}\right)^{0.5}=\frac{\dot{\lambda}}{\overline{ \sigma}}\left(\frac{3s_{ij}s_{ij}}{2}\right)^{0.5}=\dot{\lambda} \tag{68}\]  
The Jaumann rate \(\dot{\alpha}^{\text{VJ}}_{ij}\) of the back stress is obtained:  
\[\dot{\alpha}^{\text{VJ}}_{ij}=\frac{2}{3}(1-\beta)E_{\text{P}}\dot{\epsilon}^{ \text{P}}_{ij}=\frac{2}{3}(1-\beta)E_{\text{P}}\frac{3s_{ij}}{2\overline{ \sigma}}\dot{\lambda}=(1-\beta)E_{\text{P}}\frac{s_{ij}}{\overline{\sigma}}\dot{\lambda} \tag{69}\]  
\[{}^{t}\overline{\sigma}=\overline{\sigma}^{*}-3G\mathrm{d}\lambda=\overline{\sigma} ^{*}-3G\mathrm{d}\bar{\mathrm{e}}^{\mathrm{P}} \tag{76}\]  
This illustrates the relationship between \({}^{t}\overline{\sigma}\), \(\overline{\sigma}^{*}\), and \(\mathrm{d}\bar{\mathrm{e}}^{\mathrm{P}}\). For linear reinforcement, where \(\sigma_{\mathrm{y}}\) is linear to \(\overline{\mathrm{e}}^{\mathrm{P}}\), \(\mathrm{d}\bar{\mathrm{e}}^{\mathrm{P}}\) can be calculated using the \(t\) time yield function.