Plasticity#

Overview#

Plasticity models describe the (typically irreversible and dissipative) history-dependent deformation of solids. NEML2’s macroscale plasticity catalog under neml2/models/solid_mechanics/plasticity/ provides small composable ingredients that you wire together with the Syntax catalog primitives in common/ to assemble a return-map step:

Stress measures — IsotropicMandelStress, MandelStress, and SR2Invariant (from common/) project the Cauchy/Mandel stress onto the effective stress used by the yield surface.
Yield surfaces — YieldFunction for classical \(J_2\)-style envelopes, GTNYieldFunction and GursonCavitation for porous metal plasticity.
Hardening models — isotropic (LinearIsotropicHardening, VoceIsotropicHardening, SlopeSaturationVoceIsotropicHardening, PowerLawIsotropicHardeningStaticRecovery), kinematic (LinearKinematicHardening, ChabochePlasticHardening, FredrickArmstrongPlasticHardening, PowerLawKinematicHardeningStaticRecovery), and rate-temperature couplings built from the Kocks-Mecking family (KocksMeckingActivationEnergy, KocksMeckingFlowSwitch, KocksMeckingFlowViscosity, KocksMeckingIntercept, KocksMeckingRateSensitivity, KocksMeckingYieldStress).
Flow directions — Normality takes the gradient of the yield function with respect to its variational arguments, or AssociativeJ2FlowDirection provides a closed-form \(J_2\) direction.
Flow rules and rates — AssociativePlasticFlow and the hardening-variable rates AssociativeIsotropicPlasticHardening / AssociativeKinematicPlasticHardening multiply the consistency parameter by the corresponding associative direction. For rate-dependent flow, PerzynaPlasticFlowRate gives the consistency parameter itself as an explicit function of the yield function.

These pieces are typically glued together by ComposedModel and closed by a complementarity condition (consistent plasticity) or a Perzyna-style rate equation (viscoplasticity), then wrapped in ImplicitUpdate for the implicit return map. For consistent plasticity the recommended closure is the hard-switch MinMapComplementarity; FBComplementarity is a smooth alternative (see below). A closed-form radial-return path is also available for the J2-elastic-perfectly-plastic case via LinearIsotropicElasticJ2TrialStressUpdate; see the rate_independent_plasticity/radial_return/ regression scenario.

Math#

Plastic flow is governed by the Karush-Kuhn-Tucker (KKT) loading/unloading conditions

\[ f^p \leq 0, \quad \dot{\gamma} \geq 0, \quad \dot{\gamma}\, f^p = 0, \]

where \(f^p\) is the yield function and \(\gamma\) is the consistency parameter.

Consistent (rate-independent) plasticity#

NEML2 enforces the KKT conditions by recasting them as a single complementarity residual that the nonlinear solver drives to its convergence tolerance. The recommended form is the minimum map, a hard switch that reduces to whichever condition is active:

\[ r = \min(\dot{\gamma},\, -f^p), \]

implemented by MinMapComplementarity with a_inequality = 'LE'. Being piecewise-linear, it gives semismooth-Newton / active-set behavior that converges robustly in return mapping.

An alternative is the smooth Fischer-Burmeister residual,

\[ r = \dot{\gamma} - f^p - \sqrt{\dot{\gamma}^{\,2} + {f^p}^{\,2}}, \]

implemented by FBComplementarity with the same a_inequality = 'LE'. It is differentiable everywhere but strongly nonlinear near the origin, which can slow Newton convergence; prefer the minimum map unless a globally smooth residual is specifically required.

Note

“Consistent” plasticity is often called rate-independent. Rate sensitivity can still enter through the yield function or hardening laws expressed in terms of internal-variable rates, so the “consistent” label is the more precise one.

Viscoplasticity#

The Perzyna regularization replaces the complementarity with an explicit power-law overstress relation,

\[ \dot{\gamma} = \left( \frac{\langle f^p \rangle}{\eta} \right)^{n}, \]

where \(\eta\) is the reference stress, \(n\) is the rate exponent, and \(\langle \cdot \rangle\) denotes the Macaulay bracket. This is the role of PerzynaPlasticFlowRate.

Effective stress and yield function#

For classical \(J_2\) plasticity the effective stress is the von Mises norm of the deviatoric (over-)stress,

(16)#\[\begin{align} \bar{\sigma} &= \sqrt{3 J_2}, \\ J_2 &= \tfrac{1}{2}\,\mathrm{dev}\,\boldsymbol{\Xi}\,:\, \mathrm{dev}\,\boldsymbol{\Xi}, \\ \boldsymbol{\Xi} &= \boldsymbol{\sigma} - \sum_i \boldsymbol{X}_i. \end{align}\]

The yield function combines the effective stress with the yield strength \(\sigma_y\), isotropic hardening \(k(\bar{\varepsilon}^p)\), and (implicitly, through \(\boldsymbol{\Xi}\)) the kinematic back stresses \(\boldsymbol{X}_i\):

\[ f^p = \bar{\sigma} - \sigma_y - k(\bar{\varepsilon}^p). \]

Associative flow#

Associative flow rules derive the rates of the plastic strain and internal variables from the maximum-dissipation principle,

(17)#\[\begin{align} \dot{\boldsymbol{\varepsilon}}^p &= \dot{\gamma}\, \frac{\partial f^p}{\partial \boldsymbol{\sigma}}, \\ \dot{\bar{\varepsilon}}^p &= -\dot{\gamma}\, \frac{\partial f^p}{\partial k}, \\ \dot{\boldsymbol{K}}^p &= \dot{\gamma}\, \frac{\partial f^p}{\partial \boldsymbol{X}}. \end{align}\]

Normality computes these gradients by symbolically differentiating a sub-model (the yield function) with respect to a list of input variables; AssociativePlasticFlow, AssociativeIsotropicPlasticHardening, and AssociativeKinematicPlasticHardening then scale the directions by \(\dot{\gamma}\).

Example: rate-independent \(J_2\) plasticity with mixed hardening#

The input file below assembles a fully implicit return-map for elastic-plastic behaviour with both isotropic and kinematic hardening: linear-elastic stress-strain, VoceIsotropicHardening isotropic hardening, LinearKinematicHardening back stress, associative \(J_2\) flow, and the minimum-map consistency condition.

Listing 2 tests/regression/solid_mechanics/rate_independent_plasticity/isokinharden/model.i#

# neml2
[Tensors]
  [end_time]
    type = Python
    expr = 'Scalar(torch.logspace(-1.0, 5.0, 20, dtype=torch.float64))'
  []
  [times]
    type = Python
    expr = 'Scalar(end_time.data.unsqueeze(0) * torch.linspace(0.0, 1.0, 100, dtype=torch.float64).unsqueeze(-1))'
  []
  [max_strain]
    type = Python
    expr = 'SR2.fill(0.1, -0.05, -0.05, 0.0, 0.0, 0.0).dynamic_batch.expand(20)'
  []
  [strains]
    type = Python
    expr = 'SR2(max_strain.data.unsqueeze(0) * torch.linspace(0.0, 1.0, 100, dtype=torch.float64).reshape(100, 1, 1))'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'times'
    prescribed_SR2_names = 'E'
    prescribed_SR2_values = 'strains'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Models]
  [isoharden]
    type = VoceIsotropicHardening
    saturated_hardening = 100
    saturation_rate = 1.2
  []
  [kinharden]
    type = LinearKinematicHardening
    hardening_modulus = 1000
    back_stress = 'X'
  []
  [elastic_strain]
    type = SR2LinearCombination
    from = 'E plastic_strain'
    to = 'elastic_strain'
    weights = '1 -1'
  []
  [elasticity]
    type = LinearIsotropicElasticity
    coefficients = '1e5 0.3'
    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
    strain = 'elastic_strain'
  []
  [mandel_stress]
    type = IsotropicMandelStress
    cauchy_stress = 'stress'
  []
  [overstress]
    type = SR2LinearCombination
    to = 'O'
    from = 'mandel_stress X'
    weights = '1 -1'
  []
  [vonmises]
    type = SR2Invariant
    invariant_type = 'VONMISES'
    tensor = 'O'
    invariant = 'effective_stress'
  []
  [yield_surface]
    type = YieldFunction
    yield_stress = 1000
    isotropic_hardening = 'isotropic_hardening'
  []
  [flow]
    type = ComposedModel
    models = 'overstress vonmises yield_surface'
  []
  [normality]
    type = Normality
    model = 'flow'
    function = 'yield_function'
    from = 'mandel_stress X isotropic_hardening'
    to = 'flow_direction kinematic_hardening_direction isotropic_hardening_direction'
  []
  [eprate]
    type = AssociativeIsotropicPlasticHardening
  []
  [Kprate]
    type = AssociativeKinematicPlasticHardening
  []
  [Eprate]
    type = AssociativePlasticFlow
  []
  [integrate_ep]
    type = ScalarBackwardEulerTimeIntegration
    variable = 'equivalent_plastic_strain'
  []
  [integrate_Kp]
    type = SR2BackwardEulerTimeIntegration
    variable = 'kinematic_plastic_strain'
  []
  [integrate_Ep]
    type = SR2BackwardEulerTimeIntegration
    variable = 'plastic_strain'
  []
  [consistency]
    type = MinMapComplementarity
    a = 'yield_function'
    a_inequality = 'LE'
    b = 'flow_rate'
  []
  [surface]
    type = ComposedModel
    models = 'isoharden kinharden elastic_strain elasticity
              mandel_stress overstress vonmises
              yield_surface normality eprate Kprate Eprate
              consistency integrate_ep integrate_Kp integrate_Ep'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'surface'
    unknowns = 'plastic_strain kinematic_plastic_strain equivalent_plastic_strain flow_rate'
    residuals = 'plastic_strain_residual kinematic_plastic_strain_residual equivalent_plastic_strain_residual complementarity'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_SR2 = 'plastic_strain kinematic_plastic_strain'
    unknowns_Scalar = 'equivalent_plastic_strain flow_rate'
  []
  [return_map]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [model]
    type = ComposedModel
    models = 'return_map elastic_strain elasticity'
    additional_outputs = 'plastic_strain kinematic_plastic_strain equivalent_plastic_strain'
  []
[]

How the pieces wire together#

Hardening laws — [isoharden] (VoceIsotropicHardening) maps equivalent_plastic_strain to isotropic_hardening; [kinharden] (LinearKinematicHardening) maps kinematic_plastic_strain to the back-stress contribution X.
Stress side — [elastic_strain] (SR2LinearCombination) subtracts the plastic strain from the total strain, [elasticity] (LinearIsotropicElasticity) gives the Cauchy stress, and [mandel_stress] (IsotropicMandelStress) lifts it to the work-conjugate Mandel stress consumed by the yield surface.
Yield surface — [overstress] builds \(\boldsymbol{\Xi} = \boldsymbol{\sigma} - \boldsymbol{X}\) as O, [vonmises] (SR2Invariant) produces the effective stress, and [yield_surface] (YieldFunction) assembles \(f^p\) from the effective stress, yield strength, and isotropic hardening.
Flow — [flow] is a ComposedModel glueing the overstress, invariant, and yield steps so that Normality can differentiate the resulting yield_function with respect to mandel_stress, X, and isotropic_hardening, producing the three associative directions. [Eprate], [Kprate], [eprate] scale each direction by the consistency parameter flow_rate.
Time integration — three time-integration blocks (SR2BackwardEulerTimeIntegration and ScalarBackwardEulerTimeIntegration) turn the rates into the residuals plastic_strain_residual, kinematic_plastic_strain_residual, equivalent_plastic_strain_residual.
Consistency — [consistency] (MinMapComplementarity) closes the system with the minimum-map residual between yield_function and flow_rate.
Return map — [surface] composes every residual-producing piece into one model. NonlinearSystem solves for plastic_strain, kinematic_plastic_strain, equivalent_plastic_strain, flow_rate using Newton + DenseLU. [return_map] (ImplicitUpdate) wraps the equation system, warm-started by [predictor] (ConstantExtrapolationPredictor). The top-level [model] re-evaluates elastic_strain and elasticity with the converged plastic strain to expose the final stress.

Switching to viscoplasticity#

Swapping the FB closure for a direct Perzyna flow-rate evaluation is the basic move from rate-independent to viscoplastic; depending on the variation you may also integrate stress in rate form rather than plastic strain. The reference recipe lives at tests/regression/solid_mechanics/viscoplasticity/perfect/model.i.

Variations shipped in the test tree#

The tests/regression/solid_mechanics/ tree carries ready-to-read compositions for the most common dialects:

rate_independent_plasticity/perfect/ — perfectly plastic, no hardening.
rate_independent_plasticity/isoharden/ — isotropic hardening only.
rate_independent_plasticity/kinharden/ — kinematic hardening only.
rate_independent_plasticity/isokinharden/ — combined (used above).
rate_independent_plasticity/radial_return/ — closed-form \(J_2\) radial-return via LinearIsotropicElasticJ2TrialStressUpdate.
rate_independent_plasticity/gurson/ — porous metal plasticity using GTNYieldFunction and GursonCavitation.
viscoplasticity/perfect/ — Perzyna viscoplastic counterpart.
recovery/ — adds isotropic/kinematic static recovery using PowerLawIsotropicHardeningStaticRecovery and PowerLawKinematicHardeningStaticRecovery.
km_flow/ — Kocks-Mecking temperature- and rate-dependent flow.

Plasticity#

Overview#

Math#

Consistent (rate-independent) plasticity#

Viscoplasticity#

Effective stress and yield function#

Associative flow#

Example: rate-independent \(J_2\) plasticity with mixed hardening#

How the pieces wire together#

Switching to viscoplasticity#

Variations shipped in the test tree#

See also#