Biswajit Banerjee

### Plane stress Drucker-Prager return algorithm

##### Introduction

Recently I’ve encountered questions on how the radial return algorithm works when applied to plane stress plasticity. There seems to be some confusion about the application of the plane strain constraint. Here’s the first part of my take on the issue based on Simo’s work from the mid-1980s.

Warning: I have used the same symbols for tensors, matrices, and vectors (Voigt notation). The meaning depends on the context.

##### Plane stress elasticity

In plane stress we assume that one of the eigenvalues of the stress tensor is zero.
The standard engineering convention is to assume that the plane defined by the basis vectors $\mathbf{e}_1$ and $\mathbf{e}_2$ contains the non-zero eigenvalues while the $\mathbf{e}_3$ vector indicates the direction in which the eigenvalue is zero. Therefore,

$$\boldsymbol{\sigma} = \begin{bmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} \equiv \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix}$$

Similarly, the engineering convention for the plane stress strain tensor leads to

$$\boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{12} & 0 \\ \varepsilon_{12} & \varepsilon_{22} & 0 \\ 0 & 0 & \varepsilon_{33} \end{bmatrix} \equiv \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \end{bmatrix}$$

where we have ignored $\varepsilon_{33}$ because it can be expressed as a function of $\varepsilon_{11}$ and $\varepsilon_{22}$.

For an isotropic material that is linear in its elastic range,

$$\boldsymbol{\sigma} = \lambda\, \text{tr}(\boldsymbol{\varepsilon})\, \mathbf{I} + 2 \mu\, \boldsymbol{\varepsilon}$$

and we can write

$$\begin{bmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} = \lambda ( \varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33}) \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} + 2\mu \begin{bmatrix} \varepsilon_{11} & \varepsilon_{12} & 0 \\ \varepsilon_{12} & \varepsilon_{22} & 0 \\ 0 & 0 & \varepsilon_{33} \end{bmatrix}$$

From the last row we can see that

$$\lambda ( \varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33} ) = -2\mu\varepsilon_{33} \quad \implies \quad \varepsilon_{33} = -\frac{\lambda}{\lambda + 2\mu}( \varepsilon_{11} + \varepsilon_{22} )$$

Therefore,

\begin{align} \begin{bmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} & = 2\mu \begin{bmatrix} \varepsilon_{11} - \varepsilon_{33} & \varepsilon_{12} & 0 \\ \varepsilon_{12} & \varepsilon_{22} - \varepsilon_{33} & 0 \\ 0 & 0 & 0 \end{bmatrix} \\ & = \frac{2\mu}{\lambda+2\mu} \begin{bmatrix} 2(\lambda+\mu)\varepsilon_{11} + \lambda\varepsilon_{22} & (\lambda+2\mu)\varepsilon_{12} & 0 \\ (\lambda+2\mu)\varepsilon_{12} & \lambda\varepsilon_{11}+2(\lambda+\mu)\varepsilon_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{align}

In engineering notation

$$\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_{22} \\ 2\varepsilon_{12} \end{bmatrix} \quad \equiv \quad \boldsymbol{\sigma} = \mathbf{C}\, \boldsymbol{\varepsilon}$$

Note carefully that the stiffness matrix $\mathbf{C}$ is not the same as the three-dimensional stiffness matrix for an isotropic linear elastic material with the zero stress rows and columns removed.

##### Drucker-Prager plasticity

The Drucker-Prager yield condition can be expressed as

$$f(\boldsymbol{\sigma}) = \sqrt{J_2} - q(p) = 0 \quad \text{where} \quad J_2 = \tfrac{1}{2} \boldsymbol{s}:\boldsymbol{s} ~,~~ p = \tfrac{I_1}{3} = \tfrac{1}{3} \text{tr}(\boldsymbol{\sigma})$$

and $q(p)$ is a scalar-valued function of the mean stress $(p)$ with $\boldsymbol{s} = \boldsymbol{\sigma} - p\, \mathbf{I} \,.$ The associated flow rule for the yield function $f(\boldsymbol{\sigma})$ is

$$\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda}\,\frac{\partial f}{\partial \boldsymbol{\sigma}}$$

The derivative of $f$ is

$$\frac{\partial f}{\partial \boldsymbol{\sigma}} = \tfrac{1}{\sqrt{2}} \frac{\partial}{\partial \boldsymbol{\sigma}}\left[\sqrt{\boldsymbol{s}:\boldsymbol{s}}\right] - \frac{dq}{dp}\,\frac{\partial p}{\partial \boldsymbol{\sigma}} = \tfrac{1}{2\sqrt{2}} \frac{1}{\sqrt{\boldsymbol{s}:\boldsymbol{s}}} \frac{\partial}{\partial \boldsymbol{\sigma}}(\boldsymbol{s}:\boldsymbol{s}) - \tfrac{1}{3} \frac{dq}{dp}\,\frac{\partial \text{tr}(\boldsymbol{\sigma})}{\partial \boldsymbol{\sigma}}$$

Now

$$\frac{\partial}{\partial \boldsymbol{\sigma}}(\boldsymbol{s}:\boldsymbol{s}) \equiv \frac{\partial}{\partial \sigma_{mn}}(s_{ij} s_{ij}) = 2 s_{ij} \frac{\partial s_{ij}}{\partial \sigma_{mn}} = 2 s_{ij} \left[\frac{\partial \sigma_{ij}}{\partial \sigma_{mn}} - \tfrac{1}{3} \frac{\partial \sigma_{kk}}{\partial \sigma_{mn}}\,\delta_{ij}\right]$$

or,

$$\frac{\partial}{\partial \boldsymbol{\sigma}}(\boldsymbol{s}:\boldsymbol{s}) \equiv 2 s_{ij} \left[\tfrac{1}{2}(\delta_{im}\delta_{jn}+\delta_{in}\delta_{jm}) - \tfrac{1}{3} \delta_{mn}\delta_{ij}\right] = 2 s_{mn} - \tfrac{2}{3} s_{kk} \delta_{mn} = 2 s_{mn} \equiv 2 \boldsymbol{s}$$

where we have used the observation that $s_{kk} \equiv \text{tr}(\boldsymbol{s}) = 0$. We also have

$$\frac{\partial \text{tr}(\boldsymbol{\sigma})}{\partial \boldsymbol{\sigma}} \equiv \frac{\partial \sigma_{kk}}{\partial \sigma_{mn}} = \delta_{mn} \equiv \mathbf{I}$$

Therefore,

$$\frac{\partial f}{\partial \boldsymbol{\sigma}} = \tfrac{1}{\sqrt{2}} \frac{\boldsymbol{s}}{\sqrt{\boldsymbol{s}:\boldsymbol{s}}} - \tfrac{1}{3} \frac{dq}{dp}\,\mathbf{I} \quad \implies \quad \dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda}\,\left[ \tfrac{1}{\sqrt{2}} \frac{\boldsymbol{s}}{\sqrt{\boldsymbol{s}:\boldsymbol{s}}} - \tfrac{1}{3} \frac{dq}{dp}\,\mathbf{I} \right]$$

Note that if $q(p) = 1/\sqrt{3} \sigma_y$ where $\sigma_y$ is the uniaxial stress yield stress, and there is no hardening, the Drucker-Prager yield condition reduces to the von Mises yield condition.

##### Plane stress Drucker-Prager yield function

Recall that in plane stress, the stress tensor can be simplified to the vector $\boldsymbol{\sigma} = [\sigma_{11} \, \sigma_{22} \, \sigma_{12}]^T$. The mean stress in plane stress can be written as

$$p = \tfrac{1}{3} (\sigma_{11} + \sigma_{22})$$

For the deviatoric stress,

$$\boldsymbol{s} = \begin{bmatrix} \sigma_{11} & \sigma_{12} & 0 \\ \sigma_{12} & \sigma_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} - \tfrac{1}{3} (\sigma_{11}+\sigma_{22}) \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

or,

$$\begin{bmatrix} s_{11} & s_{12} & 0 \\ s_{12} & s_{22} & 0 \\ 0 & 0 & s_{33} \end{bmatrix} = \tfrac{1}{3}\begin{bmatrix} 2\sigma_{11} - \sigma_{22} & \sigma_{12} & 0 \\ \sigma_{12} & -\sigma_{11} + 2\sigma_{22} & 0 \\ 0 & 0 & -\sigma_{11} - \sigma_{22} \end{bmatrix}$$

We can write the above in vector notation as

$$\begin{bmatrix} s_{11} \\ s_{22} \\ s_{12} \end{bmatrix} = \tfrac{1}{3} \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 3 \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} \quad \leftrightarrow \quad \boldsymbol{s} = \mathbf{T}\,\boldsymbol{\sigma} \,.$$

To express the Drucker-Prager yield function in terms of these vectors, using the matrix form of $\boldsymbol{s}$, we have

$$J_2 = \tfrac{1}{2} \boldsymbol{s}:\boldsymbol{s} \equiv \tfrac{1}{2} \text{tr}(\boldsymbol{s}^T \boldsymbol{s}) = \tfrac{1}{2} [ (s_{11}^2 + s_{12}^2) + (s_{12}^2 + s_{22}^2) + (s_{11} + s_{22})^2 ]$$

We can write this in matrix form as

\begin{align} J_2 & = \tfrac{1}{2}\begin{bmatrix} s_{11} & s_{22} & s_{12} \end{bmatrix} \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 0 \\ 0 & 0 & 2 \end{bmatrix} \begin{bmatrix} s_{11} \\ s_{22} \\ s_{12} \end{bmatrix} \\ & = \tfrac{1}{6}\begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{12} \end{bmatrix} \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 6 \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} =: \tfrac{1}{2} \boldsymbol{\sigma}^T \,\mathbf{P} \,\boldsymbol{\sigma} \end{align}

Therefore, the plane stress Drucker-Prager yield function can be written as

$$f(\boldsymbol{\sigma}) = \sqrt{\tfrac{1}{2} \boldsymbol{\sigma}^T \,\mathbf{P} \,\boldsymbol{\sigma}} - q(p) = 0$$

This form of the flow yield function is commonly seen in the plasticity literature.

##### Plane stress Drucker-Prager flow rule

Recall that the Drucker-Prager flow rule is

$$\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda}\,\left[ \tfrac{1}{\sqrt{2}} \frac{\boldsymbol{s}}{\sqrt{\boldsymbol{s}:\boldsymbol{s}}} - \tfrac{1}{3} \frac{dq}{dp}\,\mathbf{I} \right]$$

In matrix form,

$$\begin{bmatrix} \dot{\varepsilon}^p_{11} & \dot{\varepsilon}^p_{12} & 0 \\ \dot{\varepsilon}^p_{12} & \dot{\varepsilon}^p_{22} & 0 \\ 0 & 0 & \dot{\varepsilon}^p_{33} \end{bmatrix} = \dot{\lambda} \left( \frac{1}{\sqrt{2 \text{tr}(\boldsymbol{s}^T\boldsymbol{s})}} \begin{bmatrix} s_{11} & s_{12} & 0 \\ s_{12} & s_{22} & 0 \\ 0 & 0 & s_{33} \end{bmatrix} - \tfrac{1}{3} \frac{dq}{dp}\,\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \right)$$

$$\boldsymbol{s}:\boldsymbol{s} = \text{tr}(\boldsymbol{s}^T\boldsymbol{s}) = \boldsymbol{\sigma}^T \,\mathbf{P} \,\boldsymbol{\sigma}$$

Therefore, in vector form,

$$\begin{bmatrix} \dot{\varepsilon_{11}}^p \\ \dot{\varepsilon_{22}}^p \\ \dot{\varepsilon_{12}}^p \end{bmatrix} = \dot{\lambda} \left( \frac{1}{\sqrt{2 \boldsymbol{\sigma}^T\,\mathbf{P}\,\boldsymbol{\sigma}}} \mathbf{T}\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} - \tfrac{1}{3} \frac{dq}{dp}\,\begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} \right)$$

If we use the standard engineering notation for strain, the flow rule becomes

$$\begin{bmatrix} \dot{\varepsilon_{11}}^p \\ \dot{\varepsilon_{22}}^p \\ 2\dot{\varepsilon_{12}}^p \end{bmatrix} = \dot{\lambda} \left( \frac{\mathbf{P}\,\boldsymbol{\sigma}}{\sqrt{2 \boldsymbol{\sigma}^T\,\mathbf{P}\,\boldsymbol{\sigma}}} - \tfrac{1}{3} \frac{dq}{dp}\,\begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} \right)$$

Therefore, the plane stress Drucker-Prager flow rule in vector form is

$$\dot{\boldsymbol{\varepsilon}}^p = \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)$$

#### Remarks

We will explore the return algorithm for the plane stress Drucker-Prager plasticity model in the next article in this series.