Hardening Plasticity

Overview

After the elastic limit of a material is exceeded and the material yields, further increases in load are usually required for plastic flow to continue. This phenomenon is known as hardening and is broadly classified as either work or strain-hardening. Hardening theories of plasticity describe the evolution of the strength through time. This notebook

• presents an introduction to hardening plasticity,
• implements a hardening "$J_2$" plasticity model in Matmodlab for the case of linear isotropic strain hardening,
• verifies the $J_2$ model against analytic solutions.

See Also

• User Defined Materials
• Linear Elastic Material
• Nonhardening Plasticity

Contents

1. Hardening Plasticity
2. $J_2$ Plasticity
3. Standard Model Implementation

Hardening Plasticity

Overview

Similar to nonhardening theories of plasticity, the hardening theory of plasticity presumes the existence of an elastic limit defined by a yield surface

$$ f\left(\pmb{\sigma}, Y\right) = 0 $$

where $f$ is the yield function. Unlike the nonhardening theory, the yield function depends not only on the mechanical stress $\pmb{\sigma}$ but also the (nonconstant) yield strength $Y$.

The rate of mechanical stress is given by

$$ \dot{\pmb{\sigma}} = \mathbb{C}{:}\dot{\pmb{\epsilon}}^{\rm e} $$

where $\mathbb{C}$ is the elastic stiffness and $\dot{\pmb{\epsilon}}^{\rm e}$ the rate of elastic strain. Presuming that the rate of strain is the sum of elastic and plastic parts, the mechanical response is defined by

$$ \dot{\pmb{\sigma}} = \mathbb{C}{:}\left(\dot{\pmb{\epsilon}} - \dot{\pmb{\epsilon}}^{\rm p}\right) $$

where $\dot{\pmb{\epsilon}}^{\rm p}$ is the rate of plastic strain. Replacing $\dot{\pmb{\epsilon}}^{\rm p}$ with $\dot{\lambda}\pmb{m}$, $\dot{\lambda}$ being the magnitude of $\dot{\pmb{\epsilon}}^{\rm p}$ and $\pmb{m}$ its direction, the mechanical response of the material is

$$ \dot{\pmb{\sigma}} = \mathbb{C}{:}\left(\dot{\pmb{\epsilon}} - \dot{\lambda}\pmb{m}\right) $$

The solution to the plasticity problem is reduced to determining $\dot{\lambda}$ such that $f\left(\pmb{\sigma}(t), Y(t)\right)\leq 0 \ \forall t>0$

Solution Process

Given the current state of stress $\pmb{\sigma}_n$, the solution to the plasticity problem begins with the hypothesis that the entire strain increment is elastic:

$$ \pmb{\sigma}_{\rm new} \stackrel{?}{=} \pmb{\sigma}_{\rm old} + \mathbb{C}{:}\dot{\pmb{\epsilon}}dt = \pmb{\sigma}^{\rm test} $$

where the descriptor "test" is used to signal the fact that at this point $\pmb{\sigma}^{\rm test}$ is merely a hypothesized solution. The hypothesis is validated if $\pmb{\sigma}^{\rm test}$ satisfies the yield condition

$$f\left(\pmb{\sigma}^{\rm test}, Y^{\rm test}\right)\leq 0$$

so that $\pmb{\sigma}_{\rm new}=\pmb{\sigma}^{\rm test}$.

If instead the hypothesis is falsefied, i.e., the predicted test stress falls outside of the yield surface defined by $f=0$, the plasticity problem,

$$\begin{align} \pmb{\sigma}_{\rm new} = \pmb{\sigma}_{\rm old} + \mathbb{C}{:}\left(\dot{\pmb{\epsilon}} - \dot{\lambda}\pmb{m}\right)dt &= \pmb{\sigma}^{\rm trial} - \dot{\lambda}\pmb{A}dt\\ f\left(\pmb{\sigma}^{\rm trial} - \dot{\lambda}\pmb{A}dt, Y\right) &= 0 \end{align}$$

where $\pmb{A}=\mathbb{C}{:}\pmb{m}$, is solved. $\dot{\pmb{\sigma}}^{\rm trial}=\mathbb{C}{:}\dot{\pmb{\epsilon}}$ is distinguished from $\dot{\pmb{\sigma}}^{\rm test}$ in that for stress driven problems $\dot{\pmb{\sigma}}^{\rm trial}$ is not necessarily known because the strain rates $\dot{\epsilon}$ are not known.

The unknown scalar $\dot{\lambda}$ is determined by noting the following observation: if $f\left(\pmb{\sigma}_{\rm old}, Y_{\rm old}\right)=0$ and, after continued loading, $f\left(\pmb{\sigma}_{\rm new}, Y_{\rm new}\right)=0$, the rate of change of $f$ itself must be zero. Thus, by the chain rule,

$$ \dot{f}{\left(\pmb{\sigma}, Y\right)} =\frac{df}{d\pmb{\sigma}}{:}\dot{\pmb{\sigma}} +\frac{df}{dY}\dot{Y} = 0 $$

During elastic loading, the yield strength does not change. Therefore, we presume that the rate of $Y$ is of the form

$$ \dot{Y} = \dot{\lambda}h_{Y} $$

where $\dot{\lambda}$ is the rate of plastic strain and $h_{Y}$ the modulus of $Y$. Since $Y$ is tied to the internal state of the material, it is regarded as a "solution-dependent variable" (SDV) and $h_Y$ a SDV modulus. SDV moduli must be determined from experiment or microphysical considerations.

Substituting the expression for $\dot{Y}$ in to the consistency condition gives

$$ \frac{df}{d\pmb{\sigma}}{:}\left(\mathbb{C}{:}\dot{\epsilon}-\dot{\lambda}\pmb{A}\right) + \frac{df}{dY}\dot{\lambda}h_Y = 0 $$

Letting

$$ \pmb{n} = \frac{df}{d\pmb{\sigma}}\Big/\Big\lVert\frac{df}{d\pmb{\sigma}}\Big\rVert $$

the preceding equation can be solved $\dot{\lambda}$, giving

$$ \dot{\lambda} = \frac{\pmb{n}{:}\mathbb{C}{:}\dot{\epsilon}}{\pmb{n}{:}\pmb{A}+H} $$

where

$$ H = -\frac{df}{dY}\Big/\Big\lVert\frac{df}{d\pmb{\sigma}}\Big\rVert h_Y $$

Substituting $\dot{\lambda}$ in to the expression for stress rate gives

$$\begin{align} \dot{\pmb{\sigma}} &= \mathbb{C}{:}\dot{\pmb{\epsilon}} - \frac{\pmb{n}{:}\mathbb{C}{:}\dot{\epsilon}}{\pmb{n}{:}\pmb{A}+H}\pmb{A}\\ &= \mathbb{C}{:}\dot{\pmb{\epsilon}} - \frac{1}{\pmb{n}{:}\pmb{A}+H}\pmb{Q}\pmb{A}{:}\dot{\pmb{\epsilon}}\\ &= \mathbb{D}{:}\dot{\pmb{\epsilon}} \end{align}$$

where

$$ \pmb{Q} = \mathbb{C}{:}\pmb{n} $$

and $$ \mathbb{D} = \mathbb{C} - \frac{1}{\pmb{n}{:}\pmb{A}+H}\pmb{Q}\pmb{A} $$

The stress rate is then integrated through time to determine $\pmb{\sigma}$

Integration Procedure

Unlike the nonhardening case, the expression for the stress rate

$$ \dot{\pmb{\sigma}} = \dot{\pmb{\sigma}}^{\rm trial} - \frac{\pmb{n}{:}\dot{\pmb{\sigma}}^{\rm trial}}{\pmb{n}{:}\pmb{A}+H}\pmb{A} $$

is not a projection of the trial stress rate, but it is still true that

$$\pmb{\sigma}_{\rm new} = \pmb{\sigma}^{\rm trial} - \Gamma\pmb{A}$$

$\Gamma$ is determined by satisfying the now evolving yield condition. In other words, $\Gamma$ is the solution to

$$f\left(\pmb{\sigma}^{\rm trial} - \Gamma\pmb{A}, Y(\Gamma)\right)=0$$

The unknown $\Gamma$ is found such that $f\left(\pmb{\sigma}(\Gamma), Y(\Gamma)\right)=0$. The solution can be found by solving the preceding equation iteratively by applying Newton's method so that

$$ \Gamma^{i+1} = \Gamma^i - \frac{f\left(\pmb{\sigma}(\Gamma^{n}), Y(\Gamma)\right)} {\frac{df}{d\pmb{\sigma}}{:}\frac{d\pmb{\sigma}}{d\Gamma} + \frac{df}{dY}\frac{dY}{d\Gamma}} = \Gamma^i + \frac{g\left(\pmb{\sigma}(\Gamma^{n}), Y(\Gamma)\right)} {\pmb{n}{:}\pmb{A}-\frac{df}{dY}\frac{dY}{d\Gamma}} $$

where $g = f\Big/\lVert df/d\pmb{\sigma}\rVert$.

When $\Gamma^{i+1}-\Gamma^i<\epsilon$, where $\epsilon$ is a small number, the iterations are complete and the updated stress can be determined.

Note that the scalar $\Gamma$ is also equal to $\Gamma=\dot{\lambda}dt$, but since $\dot{\lambda}=0$ for elastic loading, $\Gamma=\dot{\lambda}dt^p$, where $dt^p$ is the plastic part of the time step. This gives $\Gamma$ the following physical interpretation: it is the magnigute of the total plastic strain increment.

$J_2$ Plasticity

The equations developed thus far are general in the sense that they apply to any material that can be modeled by hardening plasticity. The equations are now specialized to the case of isotropic hypoelasticity and $J_2$ plasticity by defining

$$\begin{align} \dot{\pmb{\sigma}} &= 3K\,\mathrm{iso}\dot{\pmb{\epsilon}}^{\rm e} + 2G\,\mathrm{dev}\dot{\pmb{\epsilon}}^{\rm e} \\ f\left(\pmb{\sigma}\right) &= \sqrt{J_2} - \frac{1}{\sqrt{3}}Y\left(\epsilon^p_{eq}\right) \end{align} $$

where $J_2$ is the second invariant of the stress deviator, defined as

$$J_2 = \frac{1}{2}\pmb{s}{:}\pmb{s}, \quad \pmb{s} = \pmb{\sigma} - \frac{1}{3}\mathrm{tr}\left(\pmb{\sigma}\right)\pmb{I}$$

and $Y\left(\epsilon^p_{eq}\right)$ is the plastic strain dependent yield strength in tension.

Additionally, we adopt the assumption of an associative flow rule wherein $\pmb{m}=\pmb{n}$. Accordingly,

$$\begin{align} \frac{df}{d\pmb{\sigma}}&=\frac{1}{2\sqrt{J_2}}\pmb{s}, &\pmb{n}=\frac{1}{\sqrt{2J_2}}\pmb{s} \\ \pmb{A}&=\frac{2G}{\sqrt{2J_2}}\pmb{s}, &\pmb{Q}=\frac{2G}{\sqrt{2J_2}}\pmb{s} \end{align}$$

Linear Hardening

Returning now to the definition of $Y$. Let

$$Y = Y_0 + Y_1 \epsilon^p_{eq}$$

where $Y_0$ is the initial strength and $Y_1$ is a fitting parameter. Taking the rate of $Y$ allows determining the SDV modulus $h_Y$:

$$\dot{Y} = Y_1\dot{\epsilon}^p_{eq} = h_Y\dot{\lambda}, \quad h_Y = \sqrt{\frac{2}{3}}Y_1$$

Required Parameters

The model as described above requires at minimum 4 parameters: 2 independent elastic moduli and the yield strength parameters $Y_0$ and $Y_1$.

Solution Dependent Variables

Evaluation of $Y$ requires storing and tracking the equivalent plastic strain $\epsilon_{eq}^p$

Power-Law Hardening

A phenomenological power-law model presumes that

$$Y = Y_0 + Y_1 \left(\epsilon^p_{eq}\right)^m$$

where $Y_0$ is the initial strength and $K_1$ is a fitting parameter. Taking the rate of $Y$ allows determining the SDV modulus $h_Y$:

$$\dot{Y} = h_Y\dot{\lambda}, \quad h_Y = \sqrt{\frac{2}{3}}mY_1\left[\frac{Y-Y_0}{Y_1}\right]^{\frac{m-1}{m}}$$

Required Parameters

The model as described above requires at minimum 4 parameters: 2 independent elastic moduli and the yield strength parameters $Y_0$, $Y_1$, and $m$.

Solution Dependent Variables

Evaluation of $Y$ requires storing and tracking the equivalent plastic strain $\epsilon_{eq}^p$

Work Hardening

A phenomenological linear work hardening model presumes that

$$Y = Y_0 + Y_1 W_P, \quad W_P = \int_0^t \pmb{\sigma}{:}\dot{\pmb{\epsilon}}^p dt$$

where $Y_0$ is the initial strength and $K_1$ is a fitting parameter. Taking the rate of $Y$ allows determining the SDV modulus $h_Y$:

$$\dot{Y} = h_Y\dot{\lambda}, \quad h_Y = mY_1\left[\frac{Y-Y_0}{Y_1}\right]^{\frac{m-1}{m}}\sqrt{\pmb{m}{:}\pmb{m}}$$

Required Parameters

The model as described above requires at minimum 4 parameters: 2 independent elastic moduli and the yield strength parameters $Y_0$, $Y_1$, and $m$.

Solution Dependent Variables

Evaluation of $Y$ requires storing and tracking the equivalent plastic strain $\epsilon_{eq}^p$

Model Implementation

The plastic material described above is implemented as HardeningPlasticMaterial in matmodlab2/materials/plastic3.py. HardeningPlasticMaterial is implemented as a subclass of the matmodlab2.core.material.Material class. HardeningPlasticMaterial defines

• name: class attribute

  Used for referencing the material model in the MaterialPointSimulator.

• eval: instance method

  Updates the material stress, stiffness (optional), and state dependent variables to the end of the time increment.

In the model, in addition to some standard functions imported from Numpy, several helper functions are imported from various locations in Matmodlab:

• matmodlab.core.tensor

  • isotropic_part, deviatoric_part: computes the isotropic and deviatoric parts of a second-order symmetric tensor stored as an array of length 6
  • magnitude: computes the magnitude $\left(\lVert x_{ij} \rVert=\sqrt{x_{ij}x_{ij}}\right)$ of a second-order symmetric tensor stored as an array of length 6
  • double_dot: computes the double dot product of second-order symmetric tensors stored as an array of length 6
  • VOIGT: mulitplier for converting tensor strain components to engineering strain components

Model Source

The source of the file can be viewed by executing the following cell.

Verification

Exercising the elastic model through a path of uniaxial stress should result in the slope of axial stress vs. axial strain being equal to the input parameter E for the elastic portion. The maximum stress should be equal to the input parameter Y.

Note: the input parameters to a standard material are given as a dictionary of name:value pairs for each paramter. Parameters not specified are initialized to a value of zero.

Linear Hardening

Nonlinear Hardening

Validation

Let us now fit experimental data

Determine elastic response

Power law model fit