Biswajit Banerjee

### Geometric closest point return algorithm

Part 8 of the series on plasticity return algorithms

##### Introduction

In Part 7, we saw that for isotropic elastic materials and perfect associated plasticity, the trial stress and the actual stress are at the shortest distance from each other in a transformed stress space. We also saw that the transformed stress can be expressed as

$$\boldsymbol{\sigma}^\star = \frac{z}{\sqrt{3\kappa}}\,\mathbf{E}_z + \frac{r}{\sqrt{2\mu}}\,\mathbf{E}_r$$

where

$$z = \tfrac{1}{\sqrt{3}}\,\text{tr}(\boldsymbol{\sigma}) ~,~~ r = \lVert \boldsymbol{s} \rVert$$

and

$$\mathbf{E}_z = \tfrac{1}{\sqrt{3}}\,\mathbf{I} ~,~~ \mathbf{E}_r = \frac{\boldsymbol{s}}{\lVert\boldsymbol{s}\rVert} \,.$$

We can show that the transformed stress vector remains geometrically unchanged (in the sense that angles are unchanged) if we express it as

$$\boldsymbol{\sigma}^\star = z\,\mathbf{E}_z + \sqrt{\frac{3\kappa}{2\mu}}\,r\,\mathbf{E}_r =: z\,\mathbf{E}_z + r'\,\mathbf{E}_r$$

This observation can be used to perform a purely geometric stress update that can be quite efficient for nonlinear Drucker-Prager plasticity and many other phenomenological plasticity models. Let us see how this concept has been applied to the Arena plasticity model.

##### The Arena yield function

The Arena model is an extension to partially saturated soils of the Arenisca model developed by R. M. Brannon and co-workers. The yield function used by this model is a nonlinear Drucker-Prager model with a compression cap. Further details can be found in the Arena theory manual.

If the volumetric and deviatoric components of the total stress are

$$\bar{p} := -\tfrac{1}{3} \text{tr}(\boldsymbol{\sigma}) \quad \text{and} \quad \boldsymbol{s} := \boldsymbol{\sigma} + \bar{p} \mathbf{I}$$

we can define

\begin{align} \bar{p}_\text{eff} & := -\tfrac{1}{3} \text{tr}(\boldsymbol{\sigma}_\text{eff}) = \bar{p} - B \bar{p}^w \\ \boldsymbol{s}_\text{eff} & := \boldsymbol{\sigma}_\text{eff} + \bar{p}_\text{eff} \mathbf{I} = \boldsymbol{\sigma} + \bar{p} \mathbf{I} = \boldsymbol{s} \\ J_2^\text{eff} & := \tfrac{1}{2} \boldsymbol{s}_\text{eff}:\boldsymbol{s}_\text{eff} \,. \end{align}

Then the Arena yield function is

$$f(\boldsymbol{\sigma}, B, \bar{p}^w, \bar{X}) = \sqrt{J_2^\text{eff}} - F_f(\bar{p}_\text{eff}) \, F_c(\bar{p}_\text{eff}, \bar{X})$$

where

$$F_f(\bar{p}_\text{eff}) = a_1 - a_3 \exp[- 3 a_2 \bar{p}_\text{eff}] + 3 a_4 \bar{p}_\text{eff}$$

and

$$F_c(\bar{p}_\text{eff}, \bar{X}) = \begin{cases} 1 & \quad \text{for}\quad 3\bar{p}_\text{eff} \le \bar{\kappa} \\ \sqrt{1 - \left(\cfrac{3\bar{p}_\text{eff} - \bar{\kappa}}{\bar{X}_\text{eff} - \bar{\kappa}}\right)^2} & \quad \text{for}\quad 3\bar{p}_\text{eff} > \bar{\kappa} \,. \end{cases}$$

Here $a_i$ are material parameters, $\bar{X}_\text{eff}(\boldsymbol{\varepsilon}^p, B, \bar{p}^w) = \bar{X} - 3B\bar{p}^w$ is the shifted form of the apparent hydrostatic compressive strength ($\bar{X}/3$) of the partially saturated material, and $\bar{\kappa}$ is the branch point at which the cap function $F_c$ starts decreasing until it reaches the hydrostatic strength point ($\bar{X}$):

$$\bar{\kappa} = 3\bar{p}_\text{eff}^\text{peak} - (3\bar{p}_\text{eff}^\text{peak} - \bar{X}_\text{eff}) R_c$$

where $\bar{p}_\text{eff}^\text{peak}$ is the maximum hydrostatic tensile stress that the material can support and $R_c$ is a cap ratio.

##### The non-hardening return algorithm

We use the closest-point return approach discussed in the previous articles in this series to compute a return to the yield surface while keeping all internal variables fixed. The pseudocode of the algorithm is listed below:

\begin{align} &\text{Require:} && \boldsymbol{\sigma}^k, \delta\boldsymbol{\varepsilon}, X^k, K^k, G^k, (\bar{p}^w)^k, \boldsymbol{s}^\text{trial}, \sqrt{J_2^\text{trial}}, r^\text{trial}, z_\text{eff}^\text{trial}, \\ & && a_1, a_2, a_3, a_4, I_1^\text{peak}, R_c \\ &\text{Procedure:} && \text{nonHardeningReturn} \\ & 1.\quad && r'_\text{trial} \leftarrow r^\text{trial} \sqrt{\cfrac{3K^k}{2G^k}} \qquad \mbox{Transform the trial r coordinate} \\ & 2.\quad && X_\text{eff}^k \leftarrow X^k + 3 (\bar{p}^w)^k \\ & 3.\quad && z_\text{eff}^\text{closest}, r'_\text{closest} \leftarrow \text{getClosestPoint}(K^k, G^k, X_\text{eff}^k, a_1, a_2, a_3, a_4,\\ & && \qquad \qquad\qquad I_1^\text{peak}, R_c, z_\text{eff}^\text{trial}, r'_\text{trial}) \\ & 4.\quad && I_1^{\text{closest}} \leftarrow \sqrt{3} z_\text{eff}^\text{closest} - 3 (\bar{p}^w)^k, \sqrt{J_2^\text{closest}} \leftarrow \sqrt{\frac{G^k}{3K^k}}\,r'_\text{closest} \\ & 5.\quad && \text{If} {\sqrt{J_2^\text{trial}} > 0} \\ & 6. && \qquad \boldsymbol{\sigma}^\text{fixed} = \tfrac{1}{3} I_1^{\text{closest}} \mathbf{I} + \frac{\sqrt{J_2^\text{closest}}}{\sqrt{J_2^\text{trial}}}\,\boldsymbol{s}^\text{trial} \\ & && \qquad \mbox{Compute updated total stress} \\ & 7.\quad && \text{Else} \\ & 8. && \qquad \boldsymbol{\sigma}^\text{fixed} = \tfrac{1}{3} I_1^{\text{closest}} \mathbf{I} + \boldsymbol{s}^\text{trial}\\ & && \qquad \mbox{Compute updated total stress from hydrostatic trial stress} \\ & 9.\quad && \text{EndIf} \\ & 10.\quad && \delta\boldsymbol{\sigma}_\text{fixed} \leftarrow \boldsymbol{\sigma}^\text{fixed} - \boldsymbol{\sigma}^k \qquad \mbox{Compute stress increment} \\ & 11.\quad && \delta\boldsymbol{\sigma}_\text{fixed}^\text{iso} \leftarrow \tfrac{1}{3} \text{tr}(\delta\boldsymbol{\sigma}_\text{fixed}) \mathbf{I}, \quad \delta\boldsymbol{\sigma}_\text{fixed}^\text{dev} \leftarrow \delta\boldsymbol{\sigma}_\text{fixed} - \delta\boldsymbol{\sigma}_\text{fixed}^\text{iso} \\ & 12.\quad && \delta\boldsymbol{\varepsilon}^{p,\text{fixed}} = \delta\boldsymbol{\varepsilon} - \frac{1}{3K^k}\,\delta\boldsymbol{\sigma}_\text{fixed}^\text{iso} - \frac{1}{2G^k}\,\delta\boldsymbol{\sigma}_\text{fixed}^\text{dev} \\ & && \quad\mbox{Compute plastic strain increment} \\ & 13.\quad && \text{Return} \qquad \boldsymbol{\sigma}^\text{fixed}, \delta\boldsymbol{\varepsilon}^{p,\text{fixed}} \\ \end{align}
##### The closest point algorithm

Let us look at the actual implementation to get a feel for how the closest point can be found geometrically.

The bisection algorithm for the closest point can be implemented as shown below.

The yield surface points are computed using

And, finally, the actual geometric closest point algorithm is

#### Remarks

This geometric algorithm is remarkably accurate and avoids complications associated with computing gradients in the transformed space. In the next article we will discuss how the animation in this was created.