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 and and check whether we can find simple expressions for those.

Finding

To find , 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 in the expression for , this equation is best solved using a root-finding method. The quadratic equation in 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 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 and . 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 (after defining ):

$$ \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 (von Mises) plasticity, the 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 and :

$$ \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.

If you have questions/comments/corrections, please contact banerjee at parresianz dot com dot zen (without the dot zen).

 