Introduction
In the previous part
of this discussion, I derived plane stress expressions for linear elasticity, the Drucker-Prager
yield function, and the associated flow rule.
Let us now review the approach used for finding the parameter \(\dot{\lambda}\) that is needed
in return algorithms, from a purely algebraic standpoint.
Warning: We are dealing only with perfect plasticity (no hardening) in this
series of articles.
Review of 3D plasticity
If we assume an additive decomposition of the elastic and plastic strains, the linear elastic
stress-strain relation can be expressed in tensor notation as
$$
\boldsymbol{\sigma} = \mathbf{C} : \boldsymbol{\varepsilon}^e
= \mathbf{C} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p)
$$
The flow rule can be written as
$$
\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial f}{\partial\boldsymbol{\sigma}} =: \dot{\lambda} \boldsymbol{n}
$$
and the consistency condition is
$$
\dot{f} = 0 \,.
$$
From the consistency condition, the flow rule, and the stress-strain relation, we have
$$
\dot{f} = \frac{\partial f}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}}
= \boldsymbol{n} :
\mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon}}^p)
= \boldsymbol{n} :
\mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\lambda} \boldsymbol{n}) = 0
$$
Solving for \(\dot{\lambda}\), we have
$$
\dot{\lambda} = \frac{\boldsymbol{n} :
\mathbf{C} : \dot{\boldsymbol{\varepsilon}}}{ \boldsymbol{n} :
\mathbf{C} : \boldsymbol{n}} \ge 0
$$
The elastic-plastic tangent modulus can be computed from the stress-strain relation:
$$
\dot{\boldsymbol{\sigma}} = \mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon}}^p) = \mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\lambda} \boldsymbol{n})
= \mathbf{C} : \left(\dot{\boldsymbol{\varepsilon}} - \frac{\boldsymbol{n} :
\mathbf{C} : \dot{\boldsymbol{\varepsilon}}}{ \boldsymbol{n} :
\mathbf{C} : \boldsymbol{n}} \boldsymbol{n}\right)
$$
or,
$$
\dot{\boldsymbol{\sigma}} = \left[\mathbf{C} -
\frac{
(\mathbf{C}: \boldsymbol{n}) \otimes \left(\boldsymbol{n} : \mathbf{C}\right)
}{
\boldsymbol{n} : \mathbf{C} : \boldsymbol{n}
}\right] : \dot{\boldsymbol{\varepsilon}} =: \mathbf{C}^{ep} : \dot{\boldsymbol{\varepsilon}}
$$
We will now specialize these results for plane stress conditions.
Plane stress plasticity
For plane stress Drucker-Prager non-hardening associated plasticity, the above relations
can be expressed in terms of matrices and vectors.
The linear elastic stress-strain relation becomes
$$
\boldsymbol{\sigma} = \mathbf{C}\, \boldsymbol{\varepsilon}^e =
\mathbf{C}\, (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p)
$$
or
$$
\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} =
\begin{bmatrix} \frac{4\mu(\lambda+\mu)}{\lambda+2\mu} & \frac{2\mu\lambda}{\lambda+2\mu} & 0 \\
\frac{2\mu\lambda}{\lambda+2\mu} & \frac{4\mu(\lambda+\mu)}{\lambda+2\mu} & 0 \\
0 & 0 & \mu
\end{bmatrix}
\begin{bmatrix} \varepsilon_{11} - \varepsilon_{11}^p \\
\varepsilon_{22} - \varepsilon_{22}^p \\
2(\varepsilon_{12} - \varepsilon_{12}^p) \end{bmatrix}
$$
The yield function is
$$
f(\boldsymbol{\sigma}) = \sqrt{\tfrac{1}{2} \boldsymbol{\sigma}^T \,\mathbf{P} \,\boldsymbol{\sigma}} - q(p) = 0
$$
and the flow rule is
$$
\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda}\,\boldsymbol{n} =
\dot{\lambda} \left(
\tfrac{1}{\sqrt{2}}
\frac{\mathbf{P}\,\boldsymbol{\sigma}}{
\sqrt{\boldsymbol{\sigma}^T\,\mathbf{P}\,\boldsymbol{\sigma}}
}
- \tfrac{1}{3} \frac{dq}{dp}\,\mathbf{I}
\right)
$$
$$\boldsymbol{n} = \begin{bmatrix} n_{11} & n_{22} & n_{12} \end{bmatrix}^T \,.$$
From the consistency condition, the flow rule, and the stress-strain relation, we have
$$
\begin{align}
\dot{f} & = n_{11}\dot{\sigma}_{11} + n_{12}\dot{\sigma}_{12} +
n_{21}\dot{\sigma}_{21} + n_{22}\dot{\sigma}_{22} =
\begin{bmatrix} n_{11} & n_{22} & 2n_{12} \end{bmatrix}
\begin{bmatrix} \dot{\sigma}_{11} \\ \dot{\sigma}_{22} \\ \dot{\sigma}_{12} \end{bmatrix} \\
& = \begin{bmatrix} n_{11} & n_{22} & 2n_{12} \end{bmatrix}
\begin{bmatrix} C_{11} & C_{12} & 0 \\ C_{12} & C_{22} & 0 \\ 0 & 0 & C_{33}
\end{bmatrix}
\begin{bmatrix} \dot{\varepsilon}_{11} - \dot{\lambda} n_{11} \\
\dot{\varepsilon}_{22} - \dot{\lambda} n_{22} \\
2(\dot{\varepsilon}_{12} - \dot{\lambda} n_{12}) \end{bmatrix}
= 0
\end{align}
$$
Solving for \(\dot{\lambda}\), we have
$$
\dot{\lambda} = \frac{\mathbf{n}^T\,\mathbf{C}\,\dot{\boldsymbol{\varepsilon}}}
{\mathbf{n}^T\,\mathbf{C}\,\mathbf{n}} \ge 0 \quad \text{where} \quad
\mathbf{n}^T = \begin{bmatrix} n_{11} & n_{22} & 2n_{12} \end{bmatrix} \,.
$$
The elastic-plastic tangent modulus is computed from the stress-strain relation
$$
\dot{\boldsymbol{\sigma}} =
\mathbf{C}\,(\dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon}}^p) =
\mathbf{C}\,(\dot{\boldsymbol{\varepsilon}} - \dot{\lambda} \boldsymbol{n})
= \mathbf{C}\,\left(\dot{\boldsymbol{\varepsilon}} -
\frac{\mathbf{n}^T\,\mathbf{C}\,\dot{\boldsymbol{\varepsilon}}}
{\mathbf{n}^T\,\mathbf{C}\,\mathbf{n}} \,\boldsymbol{n}\right)
$$
to get
$$
\dot{\boldsymbol{\sigma}} =
\left(\mathbf{C} -
\frac{\mathbf{C}\,\boldsymbol{n}\,\mathbf{n}^T\,\mathbf{C}}
{\mathbf{n}^T\,\mathbf{C}\,\mathbf{n}}
\right) \dot{\boldsymbol{\varepsilon}} = \mathbf{C}^{ep} \, \dot{\boldsymbol{\varepsilon}}
$$
Vectors \(\boldsymbol{n}\) and \(\mathbf{n}\) are related by
$$
\mathbf{n} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \end{bmatrix} \boldsymbol{n}
$$
Next we will see the return algorithm for this problem.
Forward Euler Return algorithm
Let us define a trial elastic stress rate
$$
\dot{\boldsymbol{\sigma}}^{\text{trial}} = \mathbf{C} : \dot{\boldsymbol{\varepsilon}}
$$
We can express the stress-strain rate relation as
$$
\dot{\boldsymbol{\sigma}} = \mathbf{C} : \dot{\boldsymbol{\varepsilon}}^e
= \mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon}}^p )
= \mathbf{C} : (\dot{\boldsymbol{\varepsilon}} - \dot{\lambda}\boldsymbol{n} )
= \dot{\boldsymbol{\sigma}}^{\text{trial}} - \dot{\lambda}\,\mathbf{C} : \boldsymbol{n}
$$
From the consistency condition
$$
\dot{f} = \boldsymbol{n} : \dot{\boldsymbol{\sigma}}
= \boldsymbol{n} : (
\dot{\boldsymbol{\sigma}}^{\text{trial}} - \dot{\lambda}\,\mathbf{C} : \boldsymbol{n})
= 0
$$
Therefore,
$$
\dot{\lambda} = \frac{\boldsymbol{n} : \dot{\boldsymbol{\sigma}}^{\text{trial}}}
{\boldsymbol{n} : \mathbf{C} : \boldsymbol{n}}
$$
And we get the stress rate in terms of the trial stress rate:
$$
\dot{\boldsymbol{\sigma}} =
\dot{\boldsymbol{\sigma}}^{\text{trial}} -
\left(\frac{\boldsymbol{n} : \dot{\boldsymbol{\sigma}}^{\text{trial}}}
{\boldsymbol{n} : \mathbf{C} : \boldsymbol{n}}\right)
\,\mathbf{C} : \boldsymbol{n}
= \left[\mathbf{I}^{(4)} - \frac{(\boldsymbol{n}:\mathbf{C})\otimes\boldsymbol{n}}
{\boldsymbol{n} : \mathbf{C} : \boldsymbol{n}}\right] :
\dot{\boldsymbol{\sigma}}^{\text{trial}}
$$
In the discrete form of the return algorithm we assume that we know the state at time \(t_n\) and
will calculate the state at time \(t_{n+1}\). To integrate the above equation we
go back to matrix notation and observe that
$$
\boldsymbol{\sigma}_{n+1}^{\text{trial}} = \boldsymbol{\sigma}_n +
\mathbf{C}\,(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_{n})
= \mathbf{C}\,(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_{n}^p)
$$
The actual stress that we want to compute is
$$
\boldsymbol{\sigma}_{n+1} = \mathbf{C}\,\boldsymbol{\varepsilon}_{n+1}^e
= \mathbf{C}\,(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_{n+1}^p)
$$
We integrate the flow rule using forward Euler to get
$$
\boldsymbol{\varepsilon}_{n+1}^p =
\boldsymbol{\varepsilon}_{n}^p + \dot{\lambda} \Delta t \, \boldsymbol{n}_{n} =
\boldsymbol{\varepsilon}_{n}^p + \Delta\lambda \, \boldsymbol{n}_{n}
$$
Therefore, we have
$$
\boldsymbol{\sigma}_{n+1}
= \mathbf{C}\,(\boldsymbol{\varepsilon}_{n+1} -
\boldsymbol{\varepsilon}_{n}^p - \Delta\lambda \, \boldsymbol{n}_{n})
= \boldsymbol{\sigma}_{n+1}^{\text{trial}} - \Delta\lambda \,\mathbf{C}\, \boldsymbol{n}_{n}
$$
We don’t know \(\Delta\lambda\) and \(\boldsymbol{\sigma}_{n+1}\) yet.
To find \(\Delta\lambda\) we use the consistency condition:
$$
f(\boldsymbol{\sigma}_{n+1}) = 0 =
\sqrt{\tfrac{1}{2} \boldsymbol{\sigma}_{n+1}^T \,\mathbf{P} \,\boldsymbol{\sigma}_{n+1}} - q(p_{n+1})
$$
At this stage we are ready to get into some unpleasant algebra.
We will explore the rest of the algebra and the backward Euler return algorithm in the next
article in this series after a short break.