Biswajit Banerjee

### Plane stress forward Euler Drucker-Prager

##### Introduction

In Part 2 of this series, we saw that a forward Euler return algorithm leads to the discretized equations

$$\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{n+1}^{\text{trial}} - \Delta\lambda \,\mathbf{C}\, \boldsymbol{n}_{n}$$

where

$$\boldsymbol{n}_n = \begin{bmatrix} n^n_{11} \\ n^n_{22} \\ n^n_{12} \end{bmatrix} = \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{P}\,\boldsymbol{\sigma_n}}{ \sqrt{\boldsymbol{\sigma_n}^T\,\mathbf{P}\,\boldsymbol{\sigma_n}} } - \tfrac{1}{3} \frac{dq}{dp}(\boldsymbol{\sigma}_n)\,\mathbf{I} \right)$$

and

$$\mathbf{C} = \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} \quad \mathbf{P} = \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 6 \end{bmatrix} \quad \boldsymbol{\sigma}_n = \begin{bmatrix} \sigma^n_{11} \\ \sigma^n_{22} \\ \sigma^n_{12} \end{bmatrix} \quad \mathbf{I} = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}$$

The discrete consistency condition is

$$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})$$

In Part 3, we saw that

$$\mathbf{C} = \mathbf{Q} \mathbf{L}_C \mathbf{Q}^T \quad \mathbf{P} = \mathbf{Q} \mathbf{L}_P \mathbf{Q}^T$$

where

$$\mathbf{Q} = \begin{bmatrix} 0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ 1 & 0 & 0 \end{bmatrix} \quad \mathbf{L}_P = \begin{bmatrix} 6 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 3 \end{bmatrix} \quad \mathbf{L}_C = \begin{bmatrix} \mu & 0 & 0 \\ 0 & \frac{E}{1-\nu} & 0 \\ 0 & 0 & 2\mu \end{bmatrix}$$

Let us now try to find $\Delta\lambda$ and $\boldsymbol{\sigma}_{n+1}$ and check whether we can find simple expressions for those.

##### Finding $\Delta\lambda$

To find $\Delta\lambda$, we use the consistency condition and express all quantities in that expression in terms of the trial stress. Then,

$$\tfrac{1}{\sqrt{2}} \sqrt{ \left[\boldsymbol{\sigma}^\text{trial}_{n+1} -\Delta\lambda \mathbf{C}\boldsymbol{n}_n\right]^T \mathbf{P} \left[\boldsymbol{\sigma}^\text{trial}_{n+1} -\Delta\lambda \mathbf{C}\boldsymbol{n}_n\right]} = q(\text{tr}[\boldsymbol{\sigma}_{n+1}^\text{trial} - \Delta\lambda\mathbf{C}\boldsymbol{n}_n])$$

or

$$\left[(\boldsymbol{\sigma}^\text{trial}_{n+1})^T -\Delta\lambda \boldsymbol{n}_n^T\mathbf{C}\right] \left[\mathbf{P}\boldsymbol{\sigma}^\text{trial}_{n+1} -\Delta\lambda \mathbf{P}\mathbf{C}\boldsymbol{n}_n\right] = 2q^2(\text{tr}[\boldsymbol{\sigma}_{n+1}^\text{trial} - \Delta\lambda\mathbf{C}\boldsymbol{n}_n])$$

or

$$(\boldsymbol{\sigma}^\text{trial}_{n+1})^T\mathbf{P}\boldsymbol{\sigma}^\text{trial}_{n+1} - 2\Delta\lambda (\boldsymbol{\sigma}^\text{trial}_{n+1})^T\mathbf{P}\mathbf{C}\boldsymbol{n}_n + (\Delta\lambda)^2 (\mathbf{C}\boldsymbol{n}_n)^T\mathbf{P}\mathbf{C}\boldsymbol{n}_n \\ = 2q^2(\text{tr}[\boldsymbol{\sigma}_{n+1}^\text{trial} - \Delta\lambda\mathbf{C}\boldsymbol{n}_n])$$

Because of the potential dependence on $\Delta\lambda$ in the expression for $q^2$, this equation is best solved using a root-finding method. The quadratic equation in $\Delta\lambda$ has two solutions, only one of which should be greater than 0 for the solution to be unique. It is not straightforward to show that this is indeed the case.

##### An attempt at simplification

Let us try to find a simpler expression for $\mathbf{C}\boldsymbol{n}_n$ and related quantities. We have

$$\mathbf{C}\boldsymbol{n}_n = \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{C}\,\mathbf{P}\,\boldsymbol{\sigma_n}}{ \sqrt{\boldsymbol{\sigma_n}^T\,\mathbf{P}\,\boldsymbol{\sigma_n}} } - \tfrac{1}{3} \frac{dq}{dp}(\boldsymbol{\sigma}_n)\,\mathbf{C}\,\mathbf{I} \right)$$

Therefore,

$$\mathbf{P}\mathbf{C}\boldsymbol{n}_n = \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{P}\mathbf{C}\,\mathbf{P}\,\boldsymbol{\sigma_n}}{ \sqrt{\boldsymbol{\sigma_n}^T\,\mathbf{P}\,\boldsymbol{\sigma_n}} } - \tfrac{1}{3} \frac{dq}{dp}(\boldsymbol{\sigma}_n)\,\mathbf{P}\mathbf{C}\,\mathbf{I} \right)$$

and

$$(\mathbf{C}\boldsymbol{n}_n)^T\mathbf{P}\mathbf{C}\boldsymbol{n}_n = \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{C}\,\mathbf{P}\,\boldsymbol{\sigma_n}}{ \sqrt{\boldsymbol{\sigma_n}^T\,\mathbf{P}\,\boldsymbol{\sigma_n}} } - \tfrac{1}{3} \frac{dq}{dp}(\boldsymbol{\sigma}_n)\,\mathbf{C}\,\mathbf{I} \right)^T \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{P}\mathbf{C}\,\mathbf{P}\,\boldsymbol{\sigma_n}}{ \sqrt{\boldsymbol{\sigma_n}^T\,\mathbf{P}\,\boldsymbol{\sigma_n}} } - \tfrac{1}{3} \frac{dq}{dp}(\boldsymbol{\sigma}_n)\,\mathbf{P}\mathbf{C}\,\mathbf{I} \right)$$

Now,

$$\mathbf{C}\mathbf{P}\boldsymbol{\sigma}_n = \mathbf{Q}\mathbf{L}_C\mathbf{L}_P\mathbf{Q}^T\boldsymbol{\sigma}_n \quad \text{and} \quad \boldsymbol{\sigma}_n^T\mathbf{P}\boldsymbol{\sigma}_n = \boldsymbol{\sigma}_n^T\mathbf{Q}\mathbf{L}_P\mathbf{Q}^T\boldsymbol{\sigma}_n$$

Define $\mathbf{R}_n := \mathbf{Q}^T\boldsymbol{\sigma}_n$ and $\mathbf{R}_{n+1}^{\text{trial}} := \mathbf{Q}^T\boldsymbol{\sigma}_{n+1}^{\text{trial}}$. Then,

$$\mathbf{C}\mathbf{P}\boldsymbol{\sigma}_n = \mathbf{Q}\mathbf{L}_C\mathbf{L}_P\mathbf{R}_n \quad , \quad \mathbf{P}\mathbf{C}\mathbf{P}\boldsymbol{\sigma}_n = \mathbf{Q}\mathbf{L}_C\mathbf{L}_P^2\mathbf{R}_n \quad \text{and} \quad \boldsymbol{\sigma}_n^T\mathbf{P}\boldsymbol{\sigma}_n = \mathbf{R}_n^T\mathbf{L}_P\mathbf{R}_n$$

Then we get the following expressions for the terms in the nonlinear equation for $\Delta\lambda$ (after defining $\mathbf{Q}_I := \mathbf{Q}^T \mathbf{I}$):

\begin{align} & (\boldsymbol{\sigma}^\text{trial}_{n+1})^T\mathbf{P}\boldsymbol{\sigma}^\text{trial}_{n+1} = (\mathbf{R}_{n+1}^\text{trial})^T\mathbf{L}_P\mathbf{R}_{n+1}^\text{trial} \\ & (\boldsymbol{\sigma}^\text{trial}_{n+1})^T\mathbf{P}\mathbf{C}\boldsymbol{n}_n = (\mathbf{R}_{n+1}^{\text{trial}})^T\mathbf{L}_C\mathbf{L}_P \left( \tfrac{1}{\sqrt{2}} \frac{\mathbf{L}_P\mathbf{R}_n}{\sqrt{\mathbf{R}_n^T\mathbf{L}_P\mathbf{R}_n}} - \tfrac{1}{3} \frac{dq}{dp}\,\mathbf{Q}_I \right)\\ & (\mathbf{C}\boldsymbol{n}_n)^T\mathbf{P}\mathbf{C}\boldsymbol{n}_n = \tfrac{1}{2} \frac{\mathbf{R}_n^T\mathbf{L}_C^2\mathbf{L}_P^3\mathbf{R}_n}{\sqrt{\mathbf{R}_n^T\mathbf{L}_P\mathbf{R}_n}} - \tfrac{2}{3\sqrt{2}} \frac{dq}{dp}\frac{\mathbf{R}_n^T\mathbf{L}_C^2\mathbf{L}_P^2}{\sqrt{\mathbf{R}_n^T\mathbf{L}_P\mathbf{R}_n}}\mathbf{Q}_{I} + \tfrac{1}{9}\left(\frac{dq}{dp}\right)^2 \mathbf{Q}_I^T\mathbf{L}_C^2\mathbf{L}_P\mathbf{Q}_I \end{align}

For $J_2$ (von Mises) plasticity, the $dq/dp$ terms are zero and further simplification is possible. However, for a general Drucker-Prager material no such benefit is found from the use of the spectral decomposition, other than efficiency of computation.

In our case, the two column matrices that can help with computational efficiency are $\mathbf{R}$ and $\mathbf{Q}_I$:

$$\mathbf{R} = \tfrac{1}{\sqrt{2}}\begin{bmatrix} \sqrt{2}\sigma_{12} \\ \sigma_{11}+\sigma_{22} \\ -\sigma_{11}+\sigma_{22} \end{bmatrix} \quad \text{and} \quad \mathbf{Q}_I = \sqrt{2} \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \,.$$

#### Remarks

These results show that the “radial return” approach is applicable only to a very specific class of plasticity models. The next parts of this series will discuss the backward Euler approach and geometric approaches for return algorithms.