Anisotropic, Elastic-Plastic Material

Constitutive Law

This MPM material is identical to an orthotropic material in the elastic regime, but can plastically deform according to a built-in, anistropic Hill yielding criterion. The two Hill plastic yield criteria available are HillStyle=1:

      [math]\displaystyle{ \sqrt{F(\sigma_{yy}-\sigma_{zz})^2 + G(\sigma_{xx}-\sigma_{zz})^2 + H(\sigma_{yy}-\sigma_{xx})^2 + 2L\tau_{yz}^2 + 2M\tau_{xz}^2 + 2N\tau_{xy}^2} = H(\alpha) }[/math]

and HillStyle=2:

      [math]\displaystyle{ F(\sigma_{yy}-\sigma_{zz})^2 + G(\sigma_{xx}-\sigma_{zz})^2 + H(\sigma_{yy}-\sigma_{xx})^2 + 2L\tau_{yz}^2 + 2M\tau_{xz}^2 + 2N\tau_{xy}^2 = \bigl(H(\alpha)\bigr)^2 }[/math]

where σ and τ are normal and shear stresses in the material axis system after rotation from the analysis coordinates, [math]\displaystyle{ H(\alpha) }[/math] is a hardening law, and [math]\displaystyle{ \alpha }[/math] is a plastic hardening variable. The reason for two styles is discussed below. The remaining constants are determined by the yield stresses:

      [math]\displaystyle{ F = {1\over 2}\left({1\over \sigma_{Y,yy}^2} + {1\over \sigma_{Y,zz}^2} - {1\over \sigma_{Y,xx}^2}\right) \qquad\qquad L = {1\over 2\tau_{Y,yz}^2} }[/math]

      [math]\displaystyle{ G= {1\over 2}\left({1\over \sigma_{Y,xx}^2} + {1\over \sigma_{Y,zz}^2} - {1\over \sigma_{Y,yy}^2}\right) \qquad\qquad M = {1\over 2\tau_{Y,xz}^2} }[/math]

      [math]\displaystyle{ H = {1\over 2}\left({1\over \sigma_{Y,xx}^2} + {1\over \sigma_{Y,yy}^2} - {1\over \sigma_{Y,zz}^2}\right) \qquad\qquad N = {1\over 2\tau_{Y,xy}^2} }[/math]

where σ_Y and τ_Y are yield stresses for loading in the indicated direction. These yields stresses are given by:

      [math]\displaystyle{ {1\over \sigma_{Y,xx}^2} = G+H \qquad {1\over \sigma_{Y,yy}^2} = F+H \qquad {1\over \sigma_{Y,zz}^2} = F+G }[/math]

Clearly only one of F, G, and H can be negative (otherwise a square yield stress would be negative). The yield surface also has to describe a contained elastic region. This condition is satisfied by rewriting the yield surfaces as

      [math]\displaystyle{ \sqrt{\vec \sigma\mathbf{A}\vec \sigma} = H(\alpha) \quad({\rm style=1}) \qquad{\rm and}\qquad \vec \sigma\mathbf{A}\vec \sigma = H(\alpha)^2 \quad({\rm style=2}) }[/math]

where [math]\displaystyle{ \mathbf{A} }[/math] (when using Voight notation for stresses) is

      [math]\displaystyle{ \mathbf{A} = \left( \begin{array}{cccccc} G+H & -H & -G & 0 & 0 & 0 \\ -H & F+H & -F & 0 & 0 & 0\\ -G & -F & F+G & 0 & 0 & 0\\ 0 & 0 & 0 & 2L & 0 & 0 \\ 0 & 0 & 0 & 0 & 2M & 0 \\ 0 & 0 & 0 & 0 & 0 & 2N \end{array} \right) }[/math]

For a valid yield surface, the Eigenvalues of matrix [math]\displaystyle{ \mathbf{A} }[/math] must be non-negative. In a more practical form, given any two yield stresses related by some ratio:

      [math]\displaystyle{ R = {\sigma_{Y,ii}\over \sigma_{Y,jj}} }[/math]

then the third yield stress is bracketed by:

      [math]\displaystyle{ {\sigma_{Y,ii}\over {1+R}} \le \sigma_{Y,kk} \le {\sigma_{Y,ii}\over {|1-R|}} }[/math]

When using this material, the code automatically verifies the normal yield stresses define a physically acceptable yield surface. The shear yield stresses are allowed to have any positive values (although some might not make sense due to symmetry considerations).

Example Yield Strengths

Some examples in various types of materials follow. If one direction is prevented from yielding by setting its yield strength to ∞, then

      [math]\displaystyle{ \sigma_{Y,jj} = \infty \quad\implies\quad R = 0 \quad\implies\quad \sigma_{Y,kk} = \sigma_{Y,ii} }[/math]

In other words, the other two directions must have the same yield stress. The elastic region inside the yield surface is bounded by two parallel planes with normal (1,-1,0)/√2 separated by [math]\displaystyle{ \sqrt{2}\sigma_{Y,kk} }[/math] on either side of the [math]\displaystyle{ \sigma_{yy}=\sigma_{xx} }[/math] plane through the origin. If two directions have the same yield strength, then the third direction must be greater than half that value:

      [math]\displaystyle{ \sigma_{Y,jj} = \sigma_{Y,ii}\quad\implies\quad R = 1 \quad\implies\quad {\sigma_{Y,ii}\over 2} \le \sigma_{Y,kk} \le \infty }[/math]

In the lower limit of [math]\displaystyle{ \sigma_{Y,kk}={\sigma_{Y,ii}/2} }[/math], the elastic region inside the yield surface is bounded by two parallel planes with normal (1,1,-2)/√6 separated by [math]\displaystyle{ \sqrt{2/3}\sigma_{Y,ii} }[/math] on either side of the [math]\displaystyle{ \sigma_{yy}+\sigma_{xx}-2\sigma_{zz}=0 }[/math] plane through the origin. For [math]\displaystyle{ \sigma_{Y,kk}\gt \sigma_{Y,ii}/2 }[/math], the elastic region is an "ellipsoidal cylinder" along the [math]\displaystyle{ \sigma_{xx}=\sigma_{yy}=\sigma_{zz} }[/math] hydrostatic stress axis, which becomes a circular cylinder for the isotropic case of three equal yield stresses.

Similarly for [math]\displaystyle{ 0\lt R\lt 1 }[/math], the limits in [math]\displaystyle{ \sigma_{Y,kk} }[/math] define two limits that degenerate into parallel planes while yeild stresses within those limits are ellipsoidal surfaces. Picking values outside the limits gives opposing hyperbolic surfaces that allow a non-physical expanding elastic domain.

This limitation of yield stress limits ability to model materials that might have high yield stresses in one plane, but a low yield stress in the thickness directions. Some examples are in-plane, random fiber composites such as medium density fiber board, particle board, fiberglass composites, or paperboard. After picking reasonable, and similar, in-plane yield stresses, the thickness direction yield stress must be:

      [math]\displaystyle{ \sigma_{Y,jj} \sim \sigma_{Y,ii} \quad\implies\quad R \sim 1 \quad\implies\quad \sigma_{Y,zz} \gtrapprox \frac{1}{2}\min(\sigma_{Y,xx},\sigma_{Y,yy}) }[/math]

If that limit does not include experimental results for thickness direction yielding, Hill plasticity may need be a good model. Overall, Hill plasticity was developed for modeling sheet metals where the three anisotropic yields stresses are expected to be a similar magnitude.

Hardening Laws

There currently available hardening laws:

      [math]\displaystyle{ {\rm Nonlinear\ 1:} \qquad H(\alpha) = 1 + K \alpha^n }[/math]

      [math]\displaystyle{ {\rm Nonlinear\ 2:} \qquad H(\alpha) = (1 + K \alpha)^n }[/math]

      [math]\displaystyle{ {\rm Exponential:} \qquad H(\alpha) = \left\{ \begin{array}{ll} 1 + \frac{K}{\beta}(1-\exp(-\beta\alpha)) & \alpha \le \alpha_{max} \\ 1 + K^*\alpha & \alpha \gt \alpha_{max}\end{array} \right. }[/math]

Note, that because [math]\displaystyle{ H'(\alpha) }[/math] in the first law is infinite for [math]\displaystyle{ \alpha=0 }[/math] when n is less than one, this law becomes unstable for n less than about 0.7.

The parameter [math]\displaystyle{ \alpha_{max} }[/math] is to prevent exponential hardening from devolving into elastic-plastic deformation. When this parameter is used, the hardening is converted to linear hardening once [math]\displaystyle{ \alpha }[/math] exceeds [math]\displaystyle{ \alpha_{max} }[/math] and the slope of the linear hardening, [math]\displaystyle{ K^* }[/math], is set equal to the slope of the exponential hardening law at [math]\displaystyle{ \alpha_{max} }[/math]:

      [math]\displaystyle{ K^* = K \exp(-\beta\alpha_{max}) }[/math]

Thus, both the hardening law and its slope are continuous at [math]\displaystyle{ \alpha_{max} }[/math]. One way to pick [math]\displaystyle{ \alpha_{max} }[/math] is to choose a minimum relative slope compared to initial slope and then set:

      [math]\displaystyle{ \alpha_{max} = -\frac{\ln s_{rel}}{\beta} \qquad{\rm where}\qquad s_{rel} = \frac{K^*}{K} }[/math]

Note that these are slopes of hardening law and not slope of a uniaxial stress strain curves. In 2D, the value of [math]\displaystyle{ \alpha_{max} }[/math] needed to match experimental final slopes will depend on if the analysis is plane stress or plane strain.

Why Two Hill Criteria

The two Hill criteria described above are:

      [math]\displaystyle{ \sqrt{\sigma\cdot \mathbf{A} \sigma } = H(\alpha) \quad {\rm and} \quad \sigma\cdot \mathbf{A} \sigma = \bigl(H(\alpha)\bigr)^2 }[/math]

where elements of [math]\displaystyle{ \mathbf{A} }[/math] are derived from material properties F through N. A discussion with implementation details for HillStyle=1 (the first one with a square root) can be found in Simo and Hughes. Unfortunately, that derivation has a step that multiplies by [math]\displaystyle{ \mathbf{A}^{-1}\mathbf{A} }[/math] but the tensor for Hill plasticity is singular. Despite this inconsistency, the overall implementation seems to work well.

To avoid this problem, the criteria can use HillStyle=2 instead (with squared terms). Implementation by this approach does not need to use [math]\displaystyle{ \mathbf{A}^{-1} }[/math]. It is now the default style for this material type. Furthermore, HillStyle=2 can be used in both 2D plane stress and plane strain analysis while HillStyle=1 can only be used in plane strain.

Material Properties

Property	Description	Units	Default
yldxx	Yield stress for axial loading in the x direction	pressure units	∞
yldyy	Yield stress for axial loading in the y direction	pressure units	∞
yldzz	Yield stress for axial loading in the z direction	pressure units	∞
yldxy	Yield stress for shear loading in the x-y plane	pressure units	∞
yldxz	Yield stress for shear loading in the x-z plane	pressure units	∞
yldyz	Yield stress for shear loading in the y-z plane	pressure units	∞
Khard	Hardening law K parameter	dimensionless	0
nhard	Hardening law n parameter. If n>0, the Nonlinear 1 law is used; if n<0, the Nonlinear 2 law is used with abs(n).	dimensionless	1
exphard	Hardening law β parameter in the exponential hardening law. Use of this parameter also switches to using exponential hardening	dimensionless	0
alphaMax	Optional [math]\displaystyle{ \alpha_{max} }[/math] where exponential hardening converts to linear hardening	dimensionless	∞
HillStyle	Use 1 for square-root criterion or 2 for squared criterion. Style 2 is required when using 2D plane stress anlaysis.	dimensionless	2
largeRotation	If used, this setting is ignored and material always uses 1	dimensionless	1 (fixed)
(other)	All other properties are identical to the properties for an othotropic material.	varies	varies

Although default tensile yield stresses are all infinite, one of them must be finite to use this material. Similarly, the combination of properties must satisfy conditions for positive definiteness described above. If not, an error message will appear and simulation will not run.

If only nhard is provided, the modeling will use Nonlinear 1 or Nonlinear 2 hardening law. If only exphard is provided, it will use the exponential hardening law. If both are given, the last one provided determines the hardening law.

History Variables

The one history variable is the value of plastic hardening variable [math]\displaystyle{ \alpha }[/math]. This variable can be archived as history variable 1.