First, I'm going to define some notation. I prefer the direct notation of modern continuum mechanics as presented in this book by Gurtin et al., so I won't be carrying around the indices that you use. Lastly, the following holds for small deformation, rate-independent plasticity with isotropic hardening.
- $\mathbf{T}$ = Cauchy stress tensor
- $\mathbf{T}_0 = \mathbf{T} - \frac{1}{3}\mathrm{tr}(\mathbf{T}) \mathbf{1}$ = Deviatoric part of Cauchy stress
- $\bar{\sigma} = \sqrt{\frac{2}{3} \mathbf{T}_0 : \mathbf{T}_0}$ = Equivalent tensile stress
- $\mathbf{N}^p = \sqrt{\frac{3}{2}}\frac{\mathbf{T}_0}{\bar{\sigma}}$ = Direction of plastic flow (assume co-directionality)
- $\mathbf{E} = \mathbf{E}^e + \mathbf{E}^p$ = Additive strain decomposition
- $\mathbf{T} = \mathbb{C}:\mathbf{E}^e = (2G)\mathbf{E}^e + \lambda\mathrm{tr}(\mathbf{E}^e)\mathbf{1}$ = Linearized Elasticity constitutive equation
- $\dot{\mathbf{E}}^p = \sqrt{\frac{3}{2}}\dot{\bar{\epsilon}}^p \mathbf{N}^p$
- $\dot{\bar{\epsilon}}^p = \sqrt{\frac{2}{3}} |\dot{\mathbf{E}}^p|$ = equivalent plastic strain rate
- $\bar{\epsilon}^p = \int_0^t \dot{\bar{\epsilon}}^p(t')\,dt'$ = accumulated plastic strain
- $Y(\bar{\epsilon}^p)$ = "flow resistance" = 1D tensile yield stress $\sigma_y$ (function of $\bar{\epsilon}^p$)
- $H(\bar{\epsilon}^p) = \frac{dY}{d\bar{\epsilon}^p}$ = strain-hardening rate
- $f(\mathbf{T}_0, Y(\bar{\epsilon}^p)) = \bar{\sigma} - Y(\bar{\epsilon}^p) \le 0 \; \forall \,\mathbf{T}, \,Y$ <-- Yield Function
Next, we write out the Kuhn-Tucker consistency condition at yield (which occurs at $f=0$):
- If, $f = 0$ then $\dot{\bar{\epsilon}}^p \dot{f} = 0$.
This may seem like a strong statement, and it is, but if you think about all of the possible types of loading from a point on the "yield surface" (defined as all stress states that give $f=0$), you will find this to be true.
Now, we apply this condition to find $\dot{\bar{\epsilon}}^p$.
First, compute $\dot{f}$.
\begin{align}
\dot{f} &= \sqrt{\frac{3}{2}} \dot{\bar{\sigma}} - \dot{Y}\\
&= \sqrt{\frac{3}{2}} \frac{\dot{\mathbf{T}}_0 : \mathbf{T}_0}{|\mathbf{T}_0|} - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p\\
&= \sqrt{\frac{3}{2}} \dot{\mathbf{T}}_0 : \mathbf{N}^p - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p \\
&= \left(\sqrt{\frac{3}{2}}\right)(2G)[\dot{\mathbf{E}} - \dot{\mathbf{E}}^p] : \mathbf{N}^p - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p \\
&= \left(\sqrt{\frac{3}{2}}\right)(2G)\dot{\mathbf{E}} : \mathbf{N}^p - (3G + H)\dot{\bar{\epsilon}}^p
\end{align}
Now, we have to consider three cases: 1) Elastic unloading, 2) Neutral Loading, 3) Plastic Loading
- Elastic Unloading: This is like releasing the stress on a body that is at its yield point. Thus, $\dot{\bar{\epsilon}}^p = 0$, and $\dot{\mathbf{E}}:\mathbf{N}^p < 0$. Therefore, $\dot{f} < 0$.
- Neutral Loading: This one is a bit trickier. Here, the loading is tangent to the yield surface. For this case, the strain rate is orthogonal to the direction of plastic flow, thus $\dot{\mathbf{E}}:\mathbf{N}^p = 0$. Here, still, we have not loaded the body plastically, so $\dot{\bar{\epsilon}}^p = 0$. In codes, this is usually handled without special consideration by either the above or below cases, with the same result.
- Plastic Loading: Here, we are loading the body plastically. This means that $\dot{\mathbf{E}}:\mathbf{N}^p > 0$. Since $f=0$ at the yield surface, and $f$ cannot have a positive value, this means that $\dot{f} = 0$. This is how we can find the value of $\dot{\bar{\epsilon}}^p$.
$$
\dot{\bar{\epsilon}}^p = \frac{\sqrt{3/2} (2G) (\dot{\mathbf{E}}:\mathbf{N}^p)}{3G + H(\bar{\epsilon}^p)}
$$
I glossed over some of the more tedius algebra and buried the assumption that $(3G + H)>0$, which is true for every case that I have encountered. You can go on from here to compute the so-called "elasto-plastic" modulus, which is necessary to do implicit integration in a finite element method.
