Introduction
In Part 6, I
explained why a backward Euler stress update and a closest point return from the trial stress
to the yield surface are closely related. More specifically, the correct updated stress
is at the shortest distance from the trial stress to the yield surface
in a 9-dimensional space that has the Euclidean distance measure
$$
\lVert \boldsymbol{\sigma} \rVert_{\mathsf{C}^{-1}} = \sqrt{\boldsymbol{\sigma}:\mathsf{C}^{-1}:\boldsymbol{\sigma}}
$$
where \(\mathsf{C}\) is the stiffness tensor. We will explore some of the implications of this idea
in this article.
Note that this particular closest-point interpretation applies only for perfect plasticity and
only associative flow rules. For hardening plasticity, the space in which the actual stress is
closest to the trial stress is different. For non-associative plasticity, it is unclear whether any
closest-point approach can be rigorously justified.
Eigendecompositions in linear elasticity
The stiffness tensor for an isotropic elastic material is
$$
\mathsf{C} = \lambda\,\mathbf{I}\otimes\mathbf{I} + 2\mu\,\mathsf{I}^s
$$
where \(\lambda,\mu\) are the Lamé elastic constants, \(\mathbf{I}\) is the rank-2 identity tensor, and
\(\mathsf{I}^s\) is the symmetric rank-4 identity tensor. The inverse of \(\mathsf{C}\) is the
compliance tensor
$$
\mathsf{C}^{-1} = \mathsf{S} =
-\frac{\lambda}{2\mu(3\lambda+2\mu)}\,\mathbf{I}\otimes\mathbf{I} + \frac{1}{2\mu}\,\mathsf{I}^s
$$
Eigendecompositions of the stiffness and compliance tensors are defined via
$$
\mathsf{C}:\mathbf{V} = \lambda \mathbf{V} ~,~~
\mathsf{S}:\mathbf{V} = \frac{1}{\lambda} \mathbf{V}
$$
where \(\lambda\) are the eigenvalues (not to be confused with the Lamé modulus)
and \(\mathbf{V}\) are rank-2 tensors that form the eigenbasis. Because of the symmetries
of the stiffness matrix, there are six or less unique eigenvalues and the corresponding
eigentensors are orthogonal, i.e.,
$$
\mathbf{V}_i : \mathbf{V}_i = 1 \quad \text{and} \quad
\mathbf{V}_i : \mathbf{V}_j = 0 \,.
$$
The stiffness and compliance tensors may then be represented as:
$$
\mathsf{C} = \sum_{i=1}^m \lambda_i \mathbf{V}_i\otimes\mathbf{V}_i ~,~~
\mathsf{S} = \sum_{i=1}^m \frac{1}{\lambda_i} \mathbf{V}_i\otimes\mathbf{V}_i
$$
where \(m\) is the number of non-zero and distinct eigenvalues. Note also that, in this
eigenbasis, the symmetric rank-4 identity tensor is
$$
\mathsf{I}^s = \sum_{i=1}^m \mathbf{V}_i\otimes\mathbf{V}_i \,.
$$
Eigenprojectors are defined as rank-4 tensors that have the property (for \(i \ne j\))
$$
\mathsf{P}_i:\mathbf{V}_i = \mathbf{V}_i ~,~~
\mathsf{P}_i:\mathbf{V}_j = \mathbf{0}
$$
If we apply the eigenprojector to the rank-4 identity tensor, we get
$$
\mathsf{P}_k = \mathsf{P}_k : \mathsf{I}^s =
\sum_{i=1}^m \mathsf{P}_k: (\mathbf{V}_i\otimes\mathbf{V}_i)
= \sum_{i=1}^m (\mathsf{P}_k: \mathbf{V}_i)\otimes\mathbf{V}_i
= (\mathsf{P}_k : \mathbf{V}_k)\otimes\mathbf{V}_k
= \mathbf{V}_k\otimes\mathbf{V}_k
$$
Therefore we may also write the eigendecomposition in terms of the eigenprojectors
$$
\mathsf{C} = \sum_{i=1}^m \lambda_i \mathsf{P}_i ~,~~
\mathsf{S} = \sum_{i=1}^m \frac{1}{\lambda_i} \mathsf{P}_i ~,~~
\mathsf{I}^s = \sum_{i=1}^m \mathsf{P}_i \,.
$$
For isotropic materials, a small amount of algebra shows that there are two unique eigenvectors
which lead to the decomposition
$$
\mathsf{C} = \lambda_1 \mathsf{P}_1 + \lambda_2 \mathsf{P}_2
\quad \text{where} \quad
\mathsf{S} = \frac{1}{\lambda_1} \mathsf{P}_1 + \frac{1}{\lambda_2} \mathsf{P}_2
$$
and
$$
\mathsf{P}_1 = \tfrac{1}{3} \mathbf{I}\otimes\mathbf{I} ~,~~
\mathsf{P}_2 = \mathsf{I}^s - \mathsf{P}_1 \,.
$$
We can now express the stiffness and compliance tensors in terms of these eigenprojections:
$$
\begin{align}
\mathsf{C} & = \left(\kappa-\tfrac{2}{3}\mu\right)\,\mathbf{I}\otimes\mathbf{I} +
2\mu\,\mathsf{I}^s \\
& = 3\kappa\left(\tfrac{1}{3}\mathbf{I}\otimes\mathbf{I}\right) +
2\mu\,\left(\mathsf{I}^s - \tfrac{1}{3}\mathbf{I}\otimes\mathbf{I}\right)
\end{align}
$$
where \(\kappa\) is the bulk modulus and \(\mu\) is the shear modulus. Also,
$$
\begin{align}
\mathsf{S}
& = \frac{1}{3}\left(\frac{1}{3\kappa} - \frac{1}{2\mu}\right)\,\mathbf{I}\otimes\mathbf{I} + \frac{1}{2\mu}\,\mathsf{I}^s \\
& = \frac{1}{3\kappa}\left(\tfrac{1}{3}\mathbf{I}\otimes\mathbf{I}\right) +
\frac{1}{2\mu}\left(\mathsf{I}^s - \tfrac{1}{3}\mathbf{I}\otimes\mathbf{I}\right)
\end{align}
$$
Therefore, we can write
$$
\mathsf{C} = 3\kappa\,\mathsf{P}^{\text{iso}} + 2\mu\,\mathsf{P}^{\text{symdev}}
\quad \text{and} \quad
\mathsf{S} = \frac{1}{3\kappa}\,\mathsf{P}^{\text{iso}} +
\frac{1}{2\mu}\,\mathsf{P}^{\text{symdev}}
$$
where \(\mathsf{P}^{\text{iso}} = \mathsf{P}_1\) and
\(\mathsf{P}^{\text{symdev}} = \mathsf{P}_2\).
It is also worth noting that if
$$
\mathsf{C}^{1/2} : \mathsf{C}^{1/2} := \mathsf{C} \quad \text{and} \quad
\mathsf{S}^{1/2} : \mathsf{S}^{1/2} := \mathsf{S} \\
$$
then, using the property that \(\mathsf{P}_1 : \mathsf{P}_2 = \mathsf{0}\),
$$
\mathsf{C}^{1/2} = \sqrt{\lambda_1} \mathsf{P}_1 + \sqrt{\lambda_2} \mathsf{P}_2
\quad \text{where} \quad
\mathsf{S}^{1/2} = \frac{1}{\sqrt{\lambda_1}} \mathsf{P}_1 + \frac{1}{\sqrt{\lambda_2}} \mathsf{P}_2
$$
In that case, we have
$$
\mathsf{C}^{1/2} = \sqrt{3\kappa}\,\mathsf{P}^{\text{iso}} + \sqrt{2\mu}\,\mathsf{P}^{\text{symdev}}
\quad \text{and} \quad
\mathsf{S}^{1/2} = \frac{1}{\sqrt{3\kappa}}\,\mathsf{P}^{\text{iso}} +
\frac{1}{\sqrt{2\mu}}\,\mathsf{P}^{\text{symdev}}
$$
Details of the transformed space for isotropic linear elasticity were worked out by M. Homel in
his 2014 PhD dissertation. We will follow his approach in this section.
The distance measure
$$
\lVert \boldsymbol{\sigma} \rVert_{\mathsf{S}} = \sqrt{\boldsymbol{\sigma}:\mathsf{S}:\boldsymbol{\sigma}}
$$
can be interpreted as a standard Euclidean distance measure in a transformed stress space by
observing that
$$
\begin{align}
\lVert \boldsymbol{\sigma} \rVert_{\mathsf{S}} & =
\sqrt{(\boldsymbol{\sigma}:\mathsf{S}^{1/2}):(\mathsf{S}^{1/2}:\boldsymbol{\sigma})}
\quad \text{where} \quad
\mathsf{S}^{1/2} : \mathsf{S}^{1/2} := \mathsf{S} \\
& = \sqrt{(\mathsf{S}^{1/2}:\boldsymbol{\sigma}):(\mathsf{S}^{1/2}:\boldsymbol{\sigma})}
\quad \text{using the major symmetry of} \quad \mathsf{S} \\
& = \sqrt{\boldsymbol{\sigma}^\star:\boldsymbol{\sigma}^\star}
= \lVert \boldsymbol{\sigma}^\star \rVert
\end{align}
$$
We would like to calculate the transformed stress tensor.
The Lode invariants and the Lode basis
The Lode basis (described by R. M. Brannon in 2009) is an alternative basis that can
be used to decompose the stress tensor. Let us define the following deviatoric quantities:
$$
\boldsymbol{s} = \text{dev}(\boldsymbol{\sigma})
= \boldsymbol{\sigma} - \tfrac{1}{3}\text{tr}(\boldsymbol{\sigma})\mathbf{I}
\quad \text{and} \quad
\boldsymbol{t} = \text{dev}(\boldsymbol{s}\cdot\boldsymbol{s})
= \boldsymbol{s}\cdot\boldsymbol{s} - \tfrac{1}{3}\text{tr}(\boldsymbol{s}\cdot\boldsymbol{s})\mathbf{I}
$$
The quantity \(\boldsymbol{t}\) is also called the Hill tensor.
The Lode invariants of a stress tensor are
$$
z = \tfrac{1}{\sqrt{3}}\,\text{tr}(\boldsymbol{\sigma}) ~,~~
r = \lVert \boldsymbol{s} \rVert ~,~~
\sin 3\theta = 3\sqrt{6}\det\left(\frac{\boldsymbol{s}}{\lVert\boldsymbol{s}\rVert}\right)
$$
These invariants are associated with an orthonormal set of unit tensors
$$
\mathbf{E}_z = \tfrac{1}{\sqrt{3}}\,\mathbf{I} ~,~~
\mathbf{E}_r = \frac{\boldsymbol{s}}{\lVert\boldsymbol{s}\rVert} ~,~~
\mathbf{E}_\theta =
\frac{\frac{\boldsymbol{t}}{\lVert\boldsymbol{t}\rVert} -
\sin 3\theta\frac{\boldsymbol{s}}{\lVert\boldsymbol{s}\rVert}}{\cos 3\theta}
$$
The stress can be expressed in terms of the Lode basis as
$$
\boldsymbol{\sigma} = z\,\mathbf{E}_z + r\,\mathbf{E}_r \,.
$$
We can now compute the transformed stress tensor:
$$
\boldsymbol{\sigma}^\star = \mathsf{S}^{1/2}:\boldsymbol{\sigma}
= \left[\frac{1}{\sqrt{3\kappa}}\,\mathsf{P}^{\text{iso}} +
\frac{1}{\sqrt{2\mu}}\,\mathsf{P}^{\text{symdev}} \right]:
( z\,\mathbf{E}_z + r\,\mathbf{E}_r) \,.
$$
We can show that
$$
\mathsf{P}^{\text{iso}}:\mathbf{E}_z = \mathbf{E}_z ~,~~
\mathsf{P}^{\text{iso}}:\mathbf{E}_r = \mathbf{0} ~,~~
\mathsf{P}^{\text{symdev}}:\mathbf{E}_z = \mathbf{0} ~,~~
\mathsf{P}^{\text{symdev}}:\mathbf{E}_r = \mathbf{E}_r
$$
Therefore,
$$
\boldsymbol{\sigma}^\star
= \frac{z}{\sqrt{3\kappa}}\,\mathbf{E}_z +
\frac{r}{\sqrt{2\mu}}\,\mathbf{E}_r
$$
So we have a straightforward way of computing stresses in the transformed space
and will use this idea in the geometrical closest point return algorithm that we
will discuss in the next part of this series.