Biswajit Banerjee

### The plane stress return algorithm

Part 2 of the Drucker-Prager return algorithm

##### 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.

#### Remarks

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.