Two dimensional plasticity with radial return

Question

Note - This is a duplicate of the question I asked over on Physics which hasn't received any attention.

What is the proper method to scale deviatoric stresses back to the yield surface in plane stress 2D?

I cannot get the radial return method to produce the desired results under plain stress conditions in 2D. I'll show this using an example in both 3D and 2D in the hope that my mistake might become apparent.

Three Dimensions

Assume an arbitrary trial Cauchy Stress tensor of:

$\sigma_{tr} = \begin{bmatrix} 90 & 20 & 0 \\ 20 & 90& 0 \\ 0 & 0& 0 \end{bmatrix}$

Then the hydrostatic component is given by: $\sigma_{tr,hyd} = \begin{bmatrix} 50 & 0 & 0 \\ 0 & 50& 0 \\ 0 & 0& 50 \end{bmatrix}$

And the deviatoric component by: $\sigma_{tr,dev} = \begin{bmatrix} 40 & 20 & 0 \\ 20 & 10& 0 \\ 0 & 0& -50 \end{bmatrix}$

We will use the Von Mises criterion to determine whether the stresses are still in the yield surface. This calculated as (for general plane stress):

$\sigma_v = \sqrt{\sigma_{11}^2-\sigma_{11}\sigma_{22}+\sigma_{22}^2+3\sigma_{12}^2} = 86.6025$

For the sake of argument let's assume a yield stress of $\sigma_y=80$. In which case the deviatoric stresses need to be scaled back by a factor of $\alpha={\sigma_y}/{\sigma_v}=0.9238$. So the new deviatoric stress is given by:

$\sigma_{dev}=\alpha*\sigma_{tr,dev} = \begin{bmatrix} 36.95 & 18.48 & 0 \\ 18.48 & 9.24& 0 \\ 0 & 0& -46.19 \end{bmatrix}$

And the scaled back

$\sigma=\sigma_{tr,hyd}+\sigma_{dev} = \begin{bmatrix} 86.95& 18.48 & 0 \\ 18.48 & 59.24& 0 \\ 0 & 0& 3.81 \end{bmatrix}$

If we calculate the corresponding Von Mises stress we get $\sigma_{v}=80$ indicating that the scaling was correct.

Two Dimensions

I haven't repeated a lot of the text here for brevity. Check page 157 here for confirmation of hydrostatic stress in 2D, it's just $(\sigma_{11}+\sigma_{22})/2$. Assuming plane stress then:

$\sigma_{tr} = \begin{bmatrix} 90 & 20 \\ 20 & 90 \\ \end{bmatrix}$ $\sigma_{tr,hyd} = \begin{bmatrix} 75 & 0 \\ 0 & 75 \\ \end{bmatrix}$ $\sigma_{tr,dev} = \begin{bmatrix} 15 & 20 \\ 20 & -15 \\ \end{bmatrix}$ $\sigma_v = 86.6025$

The problem is that scaling the deviatoric stresses ($\alpha$ is the same as before) results in:

$\sigma_{dev}=\alpha*\sigma_{tr,dev} = \begin{bmatrix} 13.86 & 18.48 \\ 18.48 & -13.86 \\ \end{bmatrix}$

And then:

$\sigma=\sigma_{tr,hyd}+\sigma_{dev} = \begin{bmatrix} 88.86& 18.48 \\ 18.48 & 61.14 \\ \end{bmatrix}$

Calculating the Von Mises stress gives a value of $\sigma_{v}=85$, which is outside the yield surface. What is the problem here? What is the proper method to scale the deviatoric stresses back to the yield surface in 2D?

Update I am slowly coming around to the ideal that this radial return method is only appropriate for full 3D or plane strain (in which case the $\sigma_{33}\neq0$). Page 25 of this seems to also indicate that this is the case, but I am still unsure.

Answer 1

Let's look at this problem mathematically a la Simo.

Let stress tensor ($\boldsymbol{\sigma}$) and its deviator ($\boldsymbol{s}$) in plane stress space be $$ \boldsymbol{\sigma} = \left[\sigma_{11}\, \sigma_{22}\, \sigma_{12}\right]^T \quad \text{and} \quad \boldsymbol{s} = \left[s_{11}\, s_{22}\, s_{12}\right]^T $$ Note that $\sigma_{33} = \sigma_{13} = \sigma_{23} = s_{23} = s_{13} = 0$. Even though $s_{33} \ne 0$, we can ignore $s_{33}$ because $s_{11} + s_{22} + s_{33} = 0$.

Let $\bar{\mathsf{T}}$ be the matrix that takes $\boldsymbol{\sigma}$ to $\boldsymbol{s}$, i.e., $\boldsymbol{s} = \bar{\mathsf{T}}\,\boldsymbol{\sigma}$. Also recall the convention for writing the plane stress strains ($\boldsymbol{\varepsilon}$) and their deviators ($\boldsymbol{e}$): $$ \boldsymbol{\varepsilon} = \left[\varepsilon_{11}\,\varepsilon_{22}\, 2\varepsilon_{12}\right]^T \quad \text{and} \quad \boldsymbol{e} = \left[e_{11}\, e_{22}\, 2\varepsilon_{12}\right]^T $$ You can show, after some algebra and making sure that the convention of the strain tensors is respected, that if we use a slightly modified version of $\bar{\mathsf{T}}$ (let's call it $\mathsf{T}$): $$ \mathsf{T} = \frac{1}{3}\begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 6 \end{bmatrix} $$ we can get relatively simple expressions for the yield condition and flow rule for plane stress.

The yield function for von Mises plasticity in 3D is (in tensor notation) $$ f = \sqrt{\frac{3\, \boldsymbol{s}:\boldsymbol{s}}{2}} - \sigma_Y $$ In plane stress, and using the matrices we have defined earlier (and a bit of algebra), you can show that the yield function can be written as $$ f = \sqrt{\frac{3\, \boldsymbol{\sigma}^T\,\mathsf{T}\,\boldsymbol{\sigma}}{2}} - \sigma_Y $$ You will get yielding at your trial stress values.

You will have to compute the return using the plane stress flow rule: $$ \dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda}\,\mathsf{T}\,\boldsymbol{\sigma} \quad \text{where} \quad \boldsymbol{\varepsilon}^p = \left[\varepsilon_{11}^p\,\varepsilon_{22}^p\, 2\varepsilon_{12}^p\right]^T $$ These modifications are needed for you to reach the right point on the yield surface in your return process.

Update

If the linear elastic relation is $$ \boldsymbol{\sigma} = \mathsf{C}\,\boldsymbol{\varepsilon}^e $$ where $\mathsf{C}$ is the plane-stress $3\times 3$ stiffness matrix, we can use the flow rule to show that $$ \boldsymbol{\sigma_{n+1}} = \left[\mathsf{C}^{-1} + \Delta\lambda\mathsf{T}\right]^{-1} \mathsf{C}^{-1} \boldsymbol{\sigma_{n+1}}^{\text{trial}} $$ We solve for $\Delta\lambda = \dot{\lambda} \Delta t$ by using the consistency condition $f(\sigma_{n+1}) = 0$. The resulting equation is a quadratic in $\Delta\lambda$ and therefore a simple scaling relation cannot be used for the plane stress radial return.