Constitutive Models and Plasticity#

This chapter develops the mathematical foundations for plasticity models in femlabpy.materials. We derive the von Mises and Drucker-Prager yield criteria, establish the return mapping algorithms, and connect them to the implementation.

1. Stress Invariants and Decomposition#

1.1 Stress tensor decomposition#

Any stress tensor \(\boldsymbol{\sigma}\) can be decomposed into volumetric (hydrostatic) and deviatoric parts:

\[ \boldsymbol{\sigma} = \underbrace{\sigma_m \mathbf{I}}_{\text{volumetric}} + \underbrace{\mathbf{s}}_{\text{deviatoric}} \]

where the mean (hydrostatic) stress is:

\[ \sigma_m = \frac{1}{3} \text{tr}(\boldsymbol{\sigma}) = \frac{1}{3}(\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) \]

and the deviatoric stress tensor is:

\[ \mathbf{s} = \boldsymbol{\sigma} - \sigma_m \mathbf{I} \]

In Voigt notation for 3D: \(\mathbf{S} = [\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{zx}]^T\), the deviatoric components are:

\[ s_{ii} = \sigma_{ii} - \sigma_m, \quad s_{ij} = \sigma_{ij} \text{ for } i \neq j \]

1.2 Principal stress invariants#

The three invariants of the stress tensor are:

\[ I_1 = \text{tr}(\boldsymbol{\sigma}) = \sigma_{xx} + \sigma_{yy} + \sigma_{zz} \]

\[ I_2 = \frac{1}{2}\left[ (\text{tr}\boldsymbol{\sigma})^2 - \text{tr}(\boldsymbol{\sigma}^2) \right] \]

\[ I_3 = \det(\boldsymbol{\sigma}) \]

The deviatoric invariants \(J_2\) and \(J_3\) are defined as:

\[ J_2 = \frac{1}{2} \mathbf{s} : \mathbf{s} = \frac{1}{2} s_{ij} s_{ij} \]

\[ J_3 = \det(\mathbf{s}) \]

Expanding \(J_2\) in terms of Voigt components:

\[ J_2 = \frac{1}{2}\left( s_{xx}^2 + s_{yy}^2 + s_{zz}^2 \right) + s_{xy}^2 + s_{yz}^2 + s_{zx}^2 \]

1.3 Von Mises equivalent stress#

The von Mises equivalent stress (or effective stress) is defined as:

\[ \sigma_{eq} = \sqrt{3 J_2} = \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} \]

This can be written explicitly as:

\[ \sigma_{eq} = \sqrt{\frac{1}{2}\left[ (\sigma_{xx}-\sigma_{yy})^2 + (\sigma_{yy}-\sigma_{zz})^2 + (\sigma_{zz}-\sigma_{xx})^2 \right] + 3(\sigma_{xy}^2 + \sigma_{yz}^2 + \sigma_{zx}^2)} \]

For plane stress (\(\sigma_{zz} = \sigma_{yz} = \sigma_{zx} = 0\)):

\[ \sigma_{eq} = \sqrt{\sigma_{xx}^2 - \sigma_{xx}\sigma_{yy} + \sigma_{yy}^2 + 3\sigma_{xy}^2} \]

This is implemented in eqstress(S).

2. Von Mises Yield Criterion#

2.1 Yield function#

The von Mises yield criterion states that yielding occurs when the equivalent stress equals the yield stress:

\[ f(\boldsymbol{\sigma}) = \sigma_{eq} - \sigma_y = \sqrt{3 J_2} - \sigma_y = 0 \]

This defines a cylindrical yield surface in principal stress space, with axis along the hydrostatic line \(\sigma_1 = \sigma_2 = \sigma_3\).

2.2 Isotropic hardening#

With linear isotropic hardening, the yield stress evolves with accumulated plastic strain:

\[ \sigma_y = \sigma_{y0} + H \bar{\varepsilon}^p \]

where:

\(\sigma_{y0}\) is the initial yield stress
\(H\) is the hardening modulus
\(\bar{\varepsilon}^p\) is the equivalent plastic strain

The increment of equivalent plastic strain is related to the plastic multiplier \(\Delta\lambda\) by:

\[ \Delta \bar{\varepsilon}^p = \Delta\lambda \]

2.3 Associative flow rule#

The plastic strain rate follows the normality rule (associated flow):

\[ \dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} = \dot{\lambda} \frac{\partial \sigma_{eq}}{\partial \boldsymbol{\sigma}} = \dot{\lambda} \frac{3 \mathbf{s}}{2 \sigma_{eq}} \]

The direction tensor (normal to the yield surface) is:

\[ \mathbf{n} = \frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{3 \mathbf{s}}{2 \sigma_{eq}} \]

2.4 Elastic-plastic tangent modulus#

The consistent (algorithmic) tangent for 3D is:

\[ \mathbf{D}^{ep} = \mathbf{D} - \frac{(\mathbf{D} : \mathbf{n}) \otimes (\mathbf{n} : \mathbf{D})}{\mathbf{n} : \mathbf{D} : \mathbf{n} + H} \]

For radial return in 3D von Mises plasticity, this simplifies because the elastic stiffness \(\mathbf{D}\) can be split into deviatoric and volumetric parts.

3. Return Mapping Algorithm#

3.1 Elastic predictor#

Given the strain increment \(\Delta\boldsymbol{\varepsilon}\), compute the trial (elastic) stress:

\[ \boldsymbol{\sigma}^{trial} = \boldsymbol{\sigma}_n + \mathbf{D} : \Delta\boldsymbol{\varepsilon} \]

where \(\boldsymbol{\sigma}_n\) is the stress at the end of the previous increment.

3.2 Yield check#

Evaluate the yield function at the trial state:

\[ f^{trial} = \sigma_{eq}^{trial} - \sigma_{y,n} \]

If \(f^{trial} \leq 0\), the step is elastic: \(\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}^{trial}\).

If \(f^{trial} > 0\), plastic correction is needed.

3.3 Plastic corrector (3D radial return)#

For J2 plasticity in 3D, the return mapping has a closed-form solution. The deviatoric trial stress is:

\[ \mathbf{s}^{trial} = \boldsymbol{\sigma}^{trial} - \sigma_m^{trial} \mathbf{I} \]

The plastic multiplier is found by enforcing the consistency condition:

\[ f(\boldsymbol{\sigma}_{n+1}) = 0 \implies \sigma_{eq,n+1} = \sigma_{y,n+1} \]

This gives:

\[ \Delta\lambda = \frac{f^{trial}}{3G + H} \]

where \(G = \frac{E}{2(1+\nu)}\) is the shear modulus.

The corrected deviatoric stress is:

\[ \mathbf{s}_{n+1} = \left(1 - \frac{3G \Delta\lambda}{\sigma_{eq}^{trial}}\right) \mathbf{s}^{trial} \]

The updated stress is:

\[ \boldsymbol{\sigma}_{n+1} = \mathbf{s}_{n+1} + \sigma_m^{trial} \mathbf{I} \]

3.4 Plane stress return mapping#

Under plane stress constraints (\(\sigma_{zz} = 0\)), the radial return is no longer closed-form. The constraint must be enforced iteratively.

The residual system involves finding \(\Delta\lambda\) such that:

\[\begin{split} \begin{cases} f = \sigma_{eq} - \sigma_y = 0 \\ \sigma_{zz} = 0 \end{cases} \end{split}\]

The implementation uses Newton-Raphson iteration on the scalar consistency equation. Functions yieldvm(S, G, dL, Sy) and dyieldvm(S, G, dL, Sy) evaluate the residual and its derivative.

The iteration proceeds:

\[ \Delta\lambda^{(k+1)} = \Delta\lambda^{(k)} - \frac{f(\Delta\lambda^{(k)})}{\frac{\partial f}{\partial \Delta\lambda}\big|_{\Delta\lambda^{(k)}}} \]

until \(|f| < \text{tol}\).

4. Drucker-Prager Plasticity#

The Drucker-Prager criterion models pressure-dependent yielding, relevant for soils, rocks, and concrete.

4.1 Yield function#

\[ f(\boldsymbol{\sigma}) = \sigma_{eq} + \phi \sigma_m - \sigma_y = \sqrt{3J_2} + \phi \cdot \frac{I_1}{3} - \sigma_y = 0 \]

where \(\phi\) is the friction parameter controlling pressure sensitivity.

When \(\phi = 0\), this reduces to von Mises plasticity.

4.2 Geometric interpretation#

In principal stress space, the Drucker-Prager surface is a cone with:

Apex on the hydrostatic axis
Cone angle determined by \(\phi\)
Von Mises cylinder as the limiting case (\(\phi \to 0\))

4.3 Flow rule and derivatives#

The gradient of the yield function:

\[ \frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{3 \mathbf{s}}{2 \sigma_{eq}} + \frac{\phi}{3} \mathbf{I} \]

This non-deviatoric term means plastic dilation (volume change) accompanies yielding.

4.4 Drucker-Prager return mapping#

The stressdp function implements a coupled Newton iteration solving for both stress correction \(\Delta\mathbf{S}\) and plastic multiplier \(\Delta\lambda\) simultaneously.

The residual vector is:

\[ \mathbf{R} = \Delta\boldsymbol{\varepsilon} - \mathbf{C} : (\Delta\boldsymbol{\sigma} + \delta\boldsymbol{\sigma}) - \Delta\lambda \frac{\partial f}{\partial \boldsymbol{\sigma}} \]

The \(4 \times 4\) tangent system (for 3 stress components + multiplier) is solved at each iteration:

\[\begin{split} \begin{bmatrix} \mathbf{C} + \Delta\lambda \frac{\partial^2 f}{\partial \boldsymbol{\sigma}^2} & \frac{\partial f}{\partial \boldsymbol{\sigma}} \\ \left(\frac{\partial f}{\partial \boldsymbol{\sigma}}\right)^T & -H \end{bmatrix} \begin{bmatrix} \delta(\Delta\boldsymbol{\sigma}) \\ \delta(\Delta\lambda) \end{bmatrix} = -\begin{bmatrix} \mathbf{R} \\ f \end{bmatrix} \end{split}\]

5. Material Row Layout#

The plasticity routines expect material parameters in a flat array:

Index	Parameter	Description
0	\(E\)	Young’s modulus
1	\(\nu\)	Poisson’s ratio
2	\(\sigma_{y0}\)	Initial yield stress
3	\(H\)	Hardening modulus
4	\(\phi\)	Friction parameter (Drucker-Prager only)

6. Implementation Summary#

Module structure#

Function	Description
`devstress(S)`	Splits stress into deviatoric and mean parts
`eqstress(S)`	Computes von Mises equivalent stress
`yieldvm(S, G, dL, Sy)`	Evaluates plane-stress consistency residual
`dyieldvm(S, G, dL, Sy)`	Derivative of consistency residual
`stressvm(S, G, Sy)`	Plane-stress von Mises return mapping
`stressdp(S, G, Sy0, dE, dS)`	Drucker-Prager return mapping

Design principles#

Pure functions: Invariant helpers do not mutate inputs
Explicit state: Trial stress, plastic multiplier, and corrected stress are separate variables
Local solvers: Each material call is self-contained; no global state caching
Modular tangent: Consistent tangent computed at solver level, not material level

Source file locations#

femlabpy/
├── materials/
│   ├── __init__.py      # Public exports, devstres alias
│   ├── invariants.py    # devstress, eqstress
│   └── plasticity.py    # yieldvm, dyieldvm, stressvm, stressdp