Crystal plasticity¶

Overview¶

Crystal plasticity in NEML2 is built around an incremental, rate-form view of single-crystal kinematics: a forward operator that, given the current spatial velocity gradient, integrates the elastic stretch, the lattice orientation, and any internal hardening variables forward in time. The crystal-plasticity catalog carves this physics into reusable primitives — kinematic operators (ElasticStrainRate, OrientationRate, PlasticDeformationRate, PlasticVorticity, PlasticSpatialVelocityGradient), constitutive maps from stress to resolved shear (ResolvedShear) and from internal state to slip strength (SingleSlipStrengthMap, DislocationObstacleStrengthMap), slip kinetics (PowerLawSlipRule), and hardening laws (LinearSingleSlipHardeningRule, VoceSingleSlipHardeningRule, PerSlipForestDislocationEvolution). The crystal geometry itself (slip systems, twin systems) is generated by an ancillary CubicCrystal data object, so users specify Miller indices, not enumerated slip systems.

Math¶

The fundamental kinematics derive from the multiplicative split of the deformation gradient:

\[ F = F^e F^p \]

so that the spatial velocity gradient is

\[ l = \dot{F} F^{-1} = \dot{F}^e F^{e-1} + F^e \dot{F}^p F^{p-1} F^{e-1}. \]

The plastic deformation \(\bar{l}^p = \dot{F}^p F^{p-1}\) defines the crystal-plasticity kinematics. NEML2 assumes the elastic stretch is small, \(F^e = (I + \varepsilon)\,R^e\), so

\[ l = \dot{\varepsilon} + \Omega^e - \Omega^e \varepsilon + \varepsilon \Omega^e + l^p + \varepsilon l^p - l^p \varepsilon, \]

with \(l^p = R^e \bar{l}^p R^{eT}\) the plastic velocity gradient pushed forward into the current configuration and \(\Omega^e = \dot{R}^e R^{eT}\) the elastic spin. Two further simplifications are made:

terms quadratic in the elastic stretch \(\varepsilon\) are neglected;
terms quadratic in the rate of elastic stretch \(\dot{\varepsilon}\) are also neglected.

The first holds well for metal plasticity; the second can break down at very high strain rates.

Defining the current orientation as the composition of the initial lattice rotation \(Q_0\) with the elastic rotation, \(Q = R^e Q_0\), the elastic spin reduces to \(\Omega^e = \dot{Q} Q^T\). Splitting the spatial velocity gradient into symmetric and skew parts (\(l = d + w\), \(l^p = d^p + w^p\)) and rearranging yields the two governing rate equations integrated by this module:

(18)¶\[\begin{align} \dot{\varepsilon} &= d - d^p - \varepsilon w + w \varepsilon \\ \dot{Q} &= \left(w - w^p - \varepsilon d^p + d^p \varepsilon\right) Q. \end{align}\]

The Cauchy stress closes the system via an (in general anisotropic) elastic relation rotated into the current configuration,

\[ \sigma = C : \varepsilon. \]

To close the system, the model also needs an expression for \(l^p\). The catalog adopts Asaro’s slip-system decomposition,

\[ l^p = \sum_{i=1}^{n_\mathrm{slip}} \dot{\gamma}_i \, Q \left(d_i \otimes n_i\right) Q^T, \]

where the slip rate \(\dot{\gamma}_i\) on each system is the user’s constitutive choice, driven by the resolved shear stress

\[ \tau_i = \sigma : Q \operatorname{sym}\left(d_i \otimes n_i\right) Q^T. \]

Internally, orientations are stored as modified Rodrigues parameters; conversion to a rotation matrix is available via RotationMatrix.

Note

The catalog supports two integration strategies for the orientation equation: fully coupled implicit (integrate elastic strain, hardening, and orientation as one nonlinear system, via WR2ImplicitExponentialTimeIntegration) and decoupled (integrate strain and hardening first, then advance orientation with an explicit exponential update). The example below uses the coupled form.

Example model composition¶

The input file below assembles a single-crystal model with isotropic elasticity, a power-law slip rule, Voce hardening on a single (shared) slip strength, and fully coupled implicit integration of elastic strain, hardening, and orientation. The driver sweeps 20 batched initial orientations under a constant deformation-rate / vorticity history.

Listing 24 tests/regression/solid_mechanics/crystal_plasticity/single_crystal_coupled/model.i¶

# neml2
# Single-crystal coupled crystal plasticity, native port of
# tests/regression/solid_mechanics/crystal_plasticity/single_crystal_coupled/model.i.
# Dynamic batch axis = (20,) (varying end-time / orientation); time axis = (100,).
[Tensors]
  # end_time = LinspaceScalar(1, 10, 20) -> shape (20,)
  [end_time]
    type = Python
    expr = 'Scalar.linspace(1.0, 10.0, 20)'
  []
  # times = LinspaceScalar(0, end_time, 100) -> shape (100, 20)
  [times]
    type = Python
    expr = 'Scalar(end_time.data.unsqueeze(0) * torch.linspace(0.0, 1.0, 100, dtype=torch.float64).unsqueeze(-1))'
  []
  # deformation_rate single = FillSR2(dxx=0.1, dyy=-0.05, dzz=-0.05) batched (20,)
  [deformation_rate_single]
    type = Python
    expr = 'SR2.fill(0.1, -0.05, -0.05, 0.0, 0.0, 0.0).dynamic_batch.expand(20)'
  []
  # deformation_rate = LinspaceSR2(d_single, d_single, 100) -> shape (100, 20, 6)
  [deformation_rate]
    type = Python
    expr = 'SR2(deformation_rate_single.data.unsqueeze(0).expand(100, 20, 6).contiguous())'
  []
  # vorticity single = FillWR2(w1=0.1, w2=-0.05, w3=-0.05) batched (20,)
  [vorticity_single]
    type = Python
    expr = 'WR2(torch.tensor([0.1, -0.05, -0.05], dtype=torch.float64).unsqueeze(0).expand(20, 3).contiguous())'
  []
  # vorticity = LinspaceWR2(w_single, w_single, 100) -> shape (100, 20, 3)
  [vorticity]
    type = Python
    expr = 'WR2(vorticity_single.data.unsqueeze(0).expand(100, 20, 3).contiguous())'
  []

  # Crystal geometry inputs: lattice parameter + slip direction + slip plane
  [a]
    type = Python
    expr = 'Scalar(1.0)'
  []
  [sdirs]
    type = Python
    expr = 'MillerIndex(torch.tensor([1, 1, 0], dtype=torch.int64))'
  []
  [splanes]
    type = Python
    expr = 'MillerIndex(torch.tensor([1, 1, 1], dtype=torch.int64))'
  []

  # Initial orientation = FillRot(R1, R2, R3, method='standard'):
  # convert standard Rodrigues r_std to modified-Rodrigues parameters via
  # r = r_std / (sqrt(|r_std|^2 + 1) + 1). Shape (20, 3).
  [initial_orientation]
    type = Python
    expr = 'Rot((lambda r: r / (torch.sqrt((r * r).sum(-1, keepdim=True) + 1.0) + 1.0))(torch.stack([torch.linspace(0.0, 0.75, 20, dtype=torch.float64), torch.linspace(0.0, -0.25, 20, dtype=torch.float64), torch.linspace(-0.1, 0.1, 20, dtype=torch.float64)], dim=-1)))'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model_with_stress'
    prescribed_time = 'times'
    force_SR2_names = 'deformation_rate'
    force_SR2_values = 'deformation_rate'
    force_WR2_names = 'vorticity'
    force_WR2_values = 'vorticity'
    ic_Rot_names = 'orientation'
    ic_Rot_values = 'initial_orientation'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Data]
  [crystal_geometry]
    type = CubicCrystal
    lattice_parameter = 'a'
    slip_directions = 'sdirs'
    slip_planes = 'splanes'
  []
[]

[Models]
  [euler_rodrigues]
    type = RotationMatrix
    from = 'orientation'
    to = 'orientation_matrix'
  []
  [elasticity]
    type = LinearIsotropicElasticity
    coefficients = '1e5 0.25'
    coefficient_types = 'YOUNGS_MODULUS POISSONS_RATIO'
    strain = 'elastic_strain'
    stress = 'cauchy_stress'
  []
  [resolved_shear]
    type = ResolvedShear
    stress = 'cauchy_stress'
  []
  [elastic_stretch]
    type = ElasticStrainRate
  []
  [plastic_spin]
    type = PlasticVorticity
  []
  [plastic_deformation_rate]
    type = PlasticDeformationRate
  []
  [orientation_rate]
    type = OrientationRate
  []
  [sum_slip_rates]
    type = SumSlipRates
  []
  [slip_rule]
    type = PowerLawSlipRule
    n = 8.0
    gamma0 = 2.0e-1
  []
  [slip_strength]
    type = SingleSlipStrengthMap
    constant_strength = 50.0
  []
  [voce_hardening]
    type = VoceSingleSlipHardeningRule
    initial_slope = 500.0
    saturated_hardening = 50.0
  []
  [integrate_slip_hardening]
    type = ScalarBackwardEulerTimeIntegration
    variable = 'slip_hardening'
  []
  [integrate_elastic_strain]
    type = SR2BackwardEulerTimeIntegration
    variable = 'elastic_strain'
  []
  [integrate_orientation]
    type = WR2ImplicitExponentialTimeIntegration
    variable = 'orientation'
  []
  [implicit_rate]
    type = ComposedModel
    models = 'euler_rodrigues elasticity orientation_rate resolved_shear
              elastic_stretch plastic_deformation_rate plastic_spin
              sum_slip_rates slip_rule slip_strength voce_hardening
              integrate_slip_hardening integrate_elastic_strain integrate_orientation'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'implicit_rate'
    unknowns = 'elastic_strain slip_hardening orientation'
    residuals = 'elastic_strain_residual slip_hardening_residual orientation_residual'
  []
[]

[Solvers]
  [newton]
    type = NewtonWithLineSearch
    max_linesearch_iterations = 5
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [cp_warmup_1]
    type = CrystalPlasticityStrainPredictor
    scale = 0.1
  []
  [cp_warmup_2]
    type = ConstantExtrapolationPredictor
    unknowns_Rot = 'orientation'
    unknowns_Scalar = 'slip_hardening'
  []
  [predictor]
    type = ComposedModel
    models = 'cp_warmup_1 cp_warmup_2'
  []
  [model]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [model_with_stress]
    type = ComposedModel
    models = 'model elasticity'
    additional_outputs = 'elastic_strain'
  []
[]

Explanation¶

The input is organized into the standard NEML2 blocks: [Tensors] builds the load history and crystal-geometry inputs, [Data] materializes the slip systems, [Models] declares the constitutive pieces, and [EquationSystems] / [Solvers] wrap them in an implicit update.

Crystal geometry. CubicCrystal enumerates the cubic slip systems given a lattice parameter, a Miller-index slip direction (1 1 0), and a Miller-index slip plane (1 1 1). Downstream models pick this up by referencing the standard name (crystal_geometry) without the user listing individual systems.

Kinematic chain. Working from the lattice frame to stress:

RotationMatrix converts the modified-Rodrigues orientation to a matrix orientation_matrix consumed by the other operators.
LinearIsotropicElasticity maps elastic_strain to cauchy_stress via a Young’s modulus / Poisson’s ratio pair (anisotropic crystal elasticity is supported by swapping in an anisotropic elasticity model from the Elasticity catalog).
ResolvedShear projects cauchy_stress onto each slip system to produce the resolved shears \(\tau_i\).
ElasticStrainRate and OrientationRate realize the right-hand sides of the two governing rate equations above.
PlasticDeformationRate and PlasticVorticity assemble \(d^p\) and \(w^p\) from the slip rates \(\dot{\gamma}_i\) (Asaro’s sum) — or use PlasticSpatialVelocityGradient directly if you need \(l^p\) as a single tensor.

Slip kinetics and hardening. Three pieces close the constitutive loop:

SingleSlipStrengthMap supplies the slip strength on every system — here a single constant 50.0 shared across systems, exposed as the variable slip_strength.
PowerLawSlipRule computes the slip rate on each system from the resolved shear and the slip strength, \(\dot{\gamma}_i = \dot{\gamma}_0 \,|\tau_i/\bar{\tau}|^{n} \operatorname{sgn}(\tau_i)\), with n = 8.0 and reference rate gamma0 = 2e-1.
SumSlipRates reduces the per-system rates to a scalar accumulated slip used as the hardening driver.
VoceSingleSlipHardeningRule evolves a single slip_hardening scalar Voce-style toward a saturation value (initial slope 500, saturated hardening 50). For density-based forms swap in PerSlipForestDislocationEvolution and DislocationObstacleStrengthMap; for linear hardening swap in LinearSingleSlipHardeningRule.

Integration and solve. Three backward-Euler-flavored integrators turn the rate equations into algebraic residuals:

SR2BackwardEulerTimeIntegration on elastic_strain,
ScalarBackwardEulerTimeIntegration on slip_hardening,
WR2ImplicitExponentialTimeIntegration on orientation — an exponential map appropriate for rotations, used here in its implicit (coupled) form.

The [EquationSystems] block names the three unknowns (elastic_strain slip_hardening orientation) and their residuals, and NewtonWithLineSearch drives them to zero with a DenseLU inner solver.

Predictor. A fresh-from-elastic guess often struggles to converge for crystal plasticity, so a two-stage predictor warms the unknowns: CrystalPlasticityStrainPredictor seeds the elastic strain from a fraction (scale = 0.1) of the imposed strain increment, and ConstantExtrapolationPredictor copies the previous-step values of orientation and slip_hardening forward. Wrapped in ImplicitUpdate, this becomes the per-step solve.

The final model_with_stress composes the implicit update with the elasticity model and an additional_outputs = 'elastic_strain' so the driver records both stress and strain over the load path.

Tip

If you need to post-process an orientation (e.g. wrap the modified-Rodrigues parameters back into the fundamental zone before writing output), append a FixOrientation model to the chain.