Engineering at Alberta Courses » Plasticity in 3D Stress State

Mathematical Modelling of Plasticity: Plasticity in 3D Stress State

Strain Decomposition and Incremental Behaviour in Multi-axial Stress States with Isotropic Hardening (No Bauschinger Effect)

The three functions that were presented in the previous section can be extended to cover a multi-axial stress state.

Multi-axial Von Mises Yield Function

The introduction of the multi-axial yield function dates back to the early years of the twentieth century with various researchers conducting experiments on the yielding of metals. As it was observed that metals are, in general, linear elastic isotropic materials and that yielding occurs independently of the hydrostatic stress, it was proposed that the material would start yielding once the deviatoric component of the elastic strain energy stored inside the material reaches a critical state. The expression for the deviatoric strain energy is equal to:

$\overline{U}_{deviatoric} = \frac{(vonMisesstress)^2}{6G}$

Where, the von Mises stress, denoted by $\sigma_{vM}$ , is an invariant of the deviatoric stress tensor. The “critical” value at which the material starts to yield can be obtained from a uniaxial state of stress experiment. In the uniaxial state of stress, $\sigma_{vM}$ is equal to $\sigma_{11}$ and yielding will occur once $\sigma_{11}=\sigma_Y$ . Thus, the critical strain energy at the onset of yield is equal to:

$\overline{U}_{deviatoriccritical} = \frac{\sigma_{Y}^2}{6G}$

Thus, in a state of stress described by a Cauchy stress matrix $\sigma$ , the yield function can be defined as follows:

$f(\sigma) = \sigma_{vM} - \sigma_Y$

In this relationship, $\sigma_{vM}$ describes the state of the stress inside the material while $\sigma_Y$ gives the maximum value that $\sigma_{vM}$ can attain before yielding. $\sigma_Y$ is presented as $\sigma_{max}$ in the stress based failure criteria section. Similar to the uniaxial case, if $f(\sigma)<0$ , then the material is exhibiting elastic deformations while if $f(\sigma)=0$ then, the material is about to exhibit plastic strain. The three dimensional version of the consistency condition presented above has the following form:

$\dot f (\sigma) = \sum_{i,j=1}^3 \Big(\frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} \dot{\sigma_{ij}}\Big) - \frac{d\sigma_Y}{d\overline{\varepsilon^p}} \dot{\overline{\varepsilon^p}} = 0$

In other words, the rates of change of the stress components $\sigma_{ij}$ are related to the rate of change of the equivalent plastic strain $\overline{\varepsilon^p}$ . The incremental form of the consistency condition is such that $f(\sigma+\Delta\sigma)=0$ . This implies that as the stress changes, the yield function should not change and should also be equal to 0. Thus:

$\Delta f = \Delta \sigma_{vM} - \Delta \sigma_Y = \sum_{i,j=1}^3 \Big(\frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} \Delta \sigma_{ij}\Big) - \frac{d\sigma_Y}{d \overline{\varepsilon^p}} \Delta \overline{\varepsilon^p} = 0$

The incremental form of the consistency condition in a multi-axial state of stress relates the increment in the increase in the applied stress represented by $\sigma_{vM}$ to the incremental change of the yield stress of the material, or the increment in the applied stress components $\Delta\sigma_{ij}$ to the increment in the accumulated plastic strain parameter $\Delta\overline{\varepsilon^p}$ . The backward Euler integration scheme is recommended for better accuracy of the results.

As shown in the stress based failure criteria section, if a coordinate system is chosen such that stress matrix is diagonal and the only possible nonzero components are $\sigma_1$ , $\sigma_2$ , and $\sigma_3$ , then the von Mises yield function is represented in the stress space by a cylinder whose axis is the line representing the hydrostatic stress states ( $\sigma_1=\sigma_2=\sigma_3$ ). To show this, we will investigate the the relationship:

$(\sigma_1 - \sigma_2)^2 + (\sigma_2-\sigma_3)^2 + (\sigma_1 - \sigma_3)^2 = 2\sigma_Y^2$

after a coordinate transformation into the new coordinate system $e'_1=\frac{1}{\sqrt{3}}\{1,1,1\}$ , $e'_2=\frac{1}{\sqrt{2}}\{1,-1,0\}$ , and $e'_3=\frac{1}{\sqrt{6}}\{1,1,-2\}$ . Obviously, $e'_1$ represents the line of the hydrostatic stress while $e'_2$ and $e'_3$ are perpendicular to that line. The coordinate transformation matrix into the new coordinate system has the form:

$Q = \begin{pmatrix} \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & 0 \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} & \frac{-2}{\sqrt{6}} \\ \end{pmatrix}$

A point with coordinates $\sigma=\{\sigma_1,\sigma_2,\sigma_3\}$ will have coordinates $a=\{a_1,a_2,a_3\}$ in the new coordinate system such that:

$\sigma = \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \end{pmatrix} = Q^T \begin{pmatrix} a_1 \\ a_2 \\ a_3 \\ \end{pmatrix}$

By replacing $\sigma_1$ , $\sigma_2$ , and $\sigma_3$ in the yield surface relationship with the coordinates $a_1$ , $a_2$ , and $a_3$ :

$a_2^2+a_3^2=\frac{2}{3}\sigma_Y^2$

Suggesting that in this new coordinate system, the relationship represents a cylinder whose axis is the vector $e'_1$ and whose radius is $\sqrt{\frac{2}{3}}\sigma_Y$ .

View Mathematica Code

Q = {{1/Sqrt[3], 1/Sqrt[3], 1/Sqrt[3]}, {1/Sqrt[2], -1/Sqrt[2], 0}, {1/Sqrt[6], 1/Sqrt[6], -2/Sqrt[6]}};
a = {a1, a2, a3};
s = Transpose[Q].a;
s // MatrixForm
Print["2*Sigma_y="]
FullSimplify[(s[[1]] - s[[2]])^2 + (s[[2]] - s[[3]])^2 + (s[[1]] - s[[3]])^2]

View Python Code

from sympy import sqrt,Matrix,simplify
import sympy as sp
a1,a2,a3 = sp.symbols("a_1 a_2 a_3")
Q = Matrix([[1/sqrt(3),1/sqrt(3),1/sqrt(3)],
   [1/sqrt(2),-1/sqrt(2),0],
   [1/sqrt(6),1/sqrt(6),-2/sqrt(6)]])
a = Matrix([a1, a2, a3])
s = Q.T*a;
display("\u03C3 =",s)
display("2*\u03C3_y=")
simplify((s[0]-s[1])**2+(s[1]-s[2])**2+(s[0]-s[2])**2)

Another method to find the radius of the cylinder is by realizing that the stress state represented by the vector $oa$ in Figure 7 with the components $\sigma_1=\sigma_Y$ , $\sigma_2=\sigma_3=0$ lies on the surface of the cylinder. The unit vector representing the hydrostatic direction can be obtained by normalizing the vector $ob$ as follows:

$n = \frac{ob}{||ob||} = \frac{1}{\sqrt{3}} \begin{pmatrix} 1 \\ 1 \\ 1 \\ \end{pmatrix}$

The radius of the cylinder is given by the norm of the vector $ab$ as follows:

$||ab|| = ||oa - (oa \cdot n)n|| = \Big|\Big|\begin{pmatrix} \sigma_Y \\ 0 \\ 0 \\ \end{pmatrix} - \frac{\sigma_Y}{3} \begin{pmatrix} 1 \\ 1 \\ 1 \\ \end{pmatrix} \Big|\Big| = \sqrt{\frac{2}{3}} \sigma_Y$

Figure 7. von Mises yield function in a multi-axial stress state represented in the stress space of the principal stresses

Multi-axial TRESCA Yield Function

Another yield function that predicts similar results is that given by the maximum shear stress as described in the stress based failure criteria section. The Tresca Yield function predicts that a material will yield when the maximum shear stress reaches a critical value $\tau_{max}$ :

$\frac{\max(|\sigma_1-\sigma_2|,|\sigma_1 - \sigma_3|, |\sigma_2-\sigma_3|)}{2} = \tau_{max}$

The “critical” value $\tau_{max}$ can be obtained from a uniaxial state of stress experiment. In the uniaxial state of stress, $\tau_{max}$ is equal to $\frac{\sigma_Y}{2}$ . Thus, in a state of stress described by a Cauchy stress matrix $\sigma$ , the Tresca yield function can be defined as follows:

This yield function is represented graphically in the stress space of the principal stresses by a hexagon. Computational algorithms that use the von Mises yield criterion and those that use the Tresca yield criterion predict almost the same behaviour. While the von Mises yield criterion has a smoother curvature, the Tresca yield criterion may be easier for hand calculations. In general, when viewed on the plane perpendicular to the hydrostatic stress, the differences are practically negligible (Figure 8).

Figure 8. View perpendicular to the hydrostatic stress showing the difference between the von Mises (circle) and the Tresca (hexagon) yield surfaces.

Isotropic Hardening Rule

In a multi-axial state of loading, the isotropic hardening rule predicts that the radius of the von Mises circle will change according to the relationship given between $\overline{\varepsilon^p}$ and $\sigma_Y$ . Usually, additional plastic strain (or plastic energy dissipation) applied to the material results in an increase in the radius of the von Mises circle (Figure 9).

Figure 9. Isotropic hardening represented by a change (usually an increase) in the radius of the circle with the increase of the plastic strain parameter $\overline{\varepsilon^p}$

Multi-axial Flow Rule

The last ingredient of incremental plasticity is to find the rate of change of the plastic strain components $\dot{\varepsilon}_{ij}^p$ associated with the rate of change of the stress components that is causing plastic flow. Experimental evidence suggest that the rates of the change of the plastic strain components are directly proportional to the applied deviatoric stress components as follows:

$\dot{\varepsilon^p} = \dot{\lambda}S \Rightarrow \dot{\varepsilon}_{ij}^p = \dot{\lambda}S_{ij}$

i.e.,

$\begin{pmatrix} \dot{\varepsilon}_{11}^p & \dot{\varepsilon}_{12}^p & \dot{\varepsilon}_{13}^p \\ \dot{\varepsilon}_{21}^p & \dot{\varepsilon}_{22}^p & \dot{\varepsilon}_{23}^p \\ \dot{\varepsilon}_{31}^p & \dot{\varepsilon}_{32}^p & \dot{\varepsilon}_{33}^p \\ \end{pmatrix} = \dot{\lambda} \begin{pmatrix} S_{11} & S_{12} & S_{13} \\ S_{21} & S_{22} & S_{23} \\ S_{31} & S_{32} & S_{33} \\ \end{pmatrix}$

The incremental form predicts that the increments in the plastic strain components are directly proportional to the deviatoric stress components in a material as follows:

$\Delta \varepsilon^p = \Delta \lambda S \Rightarrow \Delta \varepsilon_{ij}^p = \Delta \lambda S_{ij}$

i.e.,

$\begin{pmatrix} \Delta \varepsilon_{11}^p & \Delta \varepsilon_{12}^p & \Delta \varepsilon_{13}^p \\ \Delta \varepsilon_{21}^p & \Delta \varepsilon_{22}^p & \Delta \varepsilon_{23}^p \\ \Delta \varepsilon_{31}^p & \Delta \varepsilon_{32}^p & \Delta \varepsilon_{33}^p \end{pmatrix} = \Delta \lambda \begin{pmatrix} S_{11} & S_{12} & S_{13} \\ S_{21} & S_{22} & S_{23} \\ S_{31} & S_{32} & S_{33} \\ \end{pmatrix}$

This observation conforms with the isochoric nature of the plastic deformation (why?). Notice that there is an inherent difference between the calculations of the plastic strain components and the elastic strain components. The increments of the elastic strain components are directly proportional to the increments in the stress components, while the plastic strain components, as observed from experiments, are directly proportional to the components of the deviatoric stress tensor.

In order to be able to compute the plastic strain components associated with a plastic flow instant, the constant of proportionality $\dot{\lambda}$ or $\Delta \lambda$ needs to be computed and then multiplied by the respective deviatoric stress component. This constant of proportionality can be obtained from energy considerations by assuming that the energy that goes into plastic deformation in a general state of stress is equal to the energy that goes into plastic deformation in a uniaxial state of stress. In a uniaxial state of the stress, the rate of the energy that is dissipated during the plastification of the material is equal to $\sigma_Y\dot{\overline{\varepsilon^p}}$ which is equal to $\sigma_{vM}\dot{\overline{\varepsilon^p}}$ . The energy dissipated during the plasticification of the material in a general state of stress is equal to $\sum_{i,j=1}^3 \left(\sigma_{ij}\dot{\varepsilon}_{ij}^p\right)$ . Assuming that the energy dissipated in the plastification process in a general state of stress is equal to that in a uniaxial state of stress results in the following:

$\sigma_{vM} \dot{\overline{\varepsilon^p}} = \sigma_{11} \dot{\varepsilon}_{11}^p + \sigma_{22} \dot{\varepsilon}_{22}^p + \sigma_{33} \dot{\varepsilon}_{33}^p + \sigma_{12} \dot{\varepsilon}_{12}^p + \sigma_{21} \dot{\varepsilon}_{21}^p + \sigma_{13} \dot{\varepsilon}_{13}^p + \sigma_{31} \dot{\varepsilon}_{31}^p + \sigma_{23} \dot{\varepsilon}_{23}^p + \sigma_{32} \dot{\varepsilon}_{32}^p$

Replacing the plastic strain components with the deviatoric stress components and the constant of proportionality $\dot{\lambda}$ results in the following expression:

$\sigma_{vM} \dot{\overline{\varepsilon^p}} = \dot{\lambda} \sum_{i,j=1}^3 \sigma_{ij}S_{ij}$

However, since $\sigma_{ij}=S_{ij}+p\delta_{ij}$ and since $\sum_{i=1}^3S_{ii}=0$ , the stress components $\sigma_{ij}$ can be replaced with the deviatoric components as follows:

$\sigma_{vM} \dot{\overline{\varepsilon^p}} = \dot{\lambda} \sum_{i,j=1}^3 S_{ij}S_{ij}$

Note that by definition, the von Mises stress is equal to:

$\sigma_{vM} = \sqrt{\frac{3}{2} \sum_{i,j=1}^3 S_{ij}S_{ij}}$

Thus, the relationship between $\dot{\lambda}$ and $\dot{\overline{\varepsilon^p}}$ is as follows:

$\dot{\lambda} = \frac{\frac{3}{2} \dot{\overline{\varepsilon^p}}}{\sigma_{vM}}$

Therefore, the plastic strain components can be calculated by integrating the following differential equations functions in $\dot{\overline{\varepsilon^p}}$ :

$\dot{\varepsilon}_{ij}^p = \dot{\overline{\varepsilon^p}} \frac{3S_{ij}}{2\sigma_{vM}}$

The equivalent plastic strain is always positive and is a measure of the history of loading of the material. The term “equivalent plastic strain” is similar to the term “equivalent stress” that is used to denote $\sigma_{vM}$ . In addition, the rate of the “equivalent plastic strain” can be related to the rates of the plastic strain components with a relationship similar in form to that between $\sigma_{vM}$ and the deviatoric stress components:

$\sum_{i,j=1}^3 \dot{\varepsilon}_{ij}^p \dot{\varepsilon}_{ij}^p = \sum_{i,j=1}^3 \Big(\frac{3}{2} \big(\dot{\overline{\varepsilon^p}}\big)^2 \frac{\frac{3}{2}S_{ij}S_{ij}}{\sigma_{vM}^2}\Big)$

Therefore:

$\dot{\overline{\varepsilon^p}} = \sqrt{\frac{2}{3} \sum_{i,j=1}^3 \dot{\varepsilon}_{ij}^p \dot{\varepsilon}_{ij}^p}$

The above relationship can be written in the following incremental form:

$\Delta \overline{\varepsilon^p} = \sqrt{\frac{2}{3} \sum_{i,j=1}^3 \Delta \varepsilon_{ij}^p \Delta \varepsilon_{ij}^p}$

Associated Flow Rule

The assumption that the incremental plastic strain components are directly proportional to the deviatoric stress components resulted in the final expression (as shown above):

$\dot{\varepsilon}_{ij}^p = \dot{\overline{\varepsilon^p}} \frac{3S_{ij}}{2\sigma_{vM}}$

It has been suggested in the early twentieth century that, similar to the existence of an “Elastic Potential Energy Function” from which the stress components and the elastic strain components are derived, the plastic strain components can also be derived from a “Plastic Potential Function” $\phi$ where:

$\dot{\varepsilon}_{ij}^p = \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{ij}}$

Careful investigation of the above expression in comparison with the previous expression for the incremental plastic strain components suggests that $\phi$ can be taken to be the yield function $\phi=f(\sigma)=\sigma_{vM}-\sigma_Y$ . In this case:

$\frac{\partial \phi}{\partial \sigma_{ij}} = \frac{\partial(\sigma_{vM}-\sigma_Y)}{\partial \sigma_{ij}} = \frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} = \frac{3S_{ij}}{2\sigma_{vM}}$

Bonus Exercise: Show the following relation:

$\frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} = \frac{3S_{ij}}{2\sigma_{vM}}$

The term “associated flow rule” is used to signify that the plastic strain components are derived from the yield function and thus, the flow rule is associated with the yield function. Associated flow rules are usually used for metal plasticity, while for plasticity of other materials, the incremental plastic strain components could be derived from different functions that are not associated with the yield function. In the latter case, the flow rule is termed: “non-associated flow rule”. Associated flow rules in metals that are derived from the von Mises or the Tresca yield surfaces predict zero volumetric plastic strain, however, in soils and many other materials, plasticity is associated with a change in volume and different “plastic potential” functions are utilized.

The Normality Rule

Consider a smooth function $\phi:\mathbb{M}^3\rightarrow\mathbb{R}$ such that $\phi:\sigma\mapsto\mathbb{R}$ . With 9 independent components, $\sigma$ can also be viewed as an element of a 9 dimensional vector space. The function $\phi$ can then be considered as a smooth scalar valued function $\phi:\mathbb{R}^9\rightarrow\mathbb{R}$ . The gradient of $\phi$ with respect to its variables is denoted by $\nabla\phi$ and is defined as a vector field in $\mathbb{R}^9$ with components $\frac{\partial \phi}{\partial \sigma_{ij}}$ (See the section on the gradient of scalar fields). The direction pointing towards constant $\phi$ in $\mathbb{R}^9$ can be obtained as the unit vector $n$ along which the directional derivative is equal to zero:

$\nabla \phi \cdot n = 0$

Thus, the vector $\nabla\phi$ is normal to the vectors pointing in the direction of constant $\phi$ . If $\phi$ is taken to be the yield function defined on the nine dimensional Euclidean vector space of the nine stress components (termed the stress space), then the vector $\nabla\phi$ is perpendicular to the direction of constant $\phi$ in the stress space. Since $\phi$ is constant ( $\phi=0$ ) on the yield surface, then $\nabla\phi$ is normal to the yield surface. The incremental plastic strain components can also be considered to be elements of a nine dimensional Euclidean vector space and the associated flow rule asserts that the nine dimensional vector $\Delta\varepsilon^p$ with components $\Delta\varepsilon_{ij}^p$ is linearly dependent on (parallel to) $\nabla\phi$ as:

$\Delta \varepsilon_{ij}^p = \Delta \overline{\varepsilon^p} \frac{\partial \phi}{\partial \sigma_{ij}} \Rightarrow \Delta \varepsilon^p = \Delta \overline{\varepsilon^p} \nabla \phi$

Thus, the nine dimensional vector $\Delta\varepsilon^p$ can also be viewed as normal (perpendicular) to the yield surface and thus, the associated flow rule leads to the normality rule. In a three dimensional vector space of the principal stresses of $\sigma$ , the components of $\Delta\varepsilon^p$ are perpendicular to the yield surface (Figure 10).

See Example 1 in the examples and exercises page for a detailed example of an isotropic hardening material model.

Open Educational Resources