Crystal-plasticity FE simulation#

Overview#

Crystal-plasticity finite-element (CPFE) solution of a reconstructed microstructure under prescribed boundary conditions, producing per-element field histories (stress, strain, Nye tensor, orientation). The stage is driven by CPFESimulation, which writes the MOOSE input decks, bakes the material parameters into a NEML2 v3 model, AOTI-compiles it with neml2-compile, and launches the MOOSE/PUMA solver (puma-opt).

Method#

The material response at each quadrature point is a NEML2 v3 single-crystal plasticity model. Plastic flow is carried entirely by crystallographic slip: the lattice itself only stretches elastically and rotates, and every bit of permanent shape change comes from dislocations gliding on specific planes in specific directions. For a cubic metal that is the \(\{111\}\langle 110 \rangle\) family, supplied by the CubicCrystal geometry.

Symbols follow Notation and conventions.

Resolved shear stress#

Elasticity is a cubic (or isotropic) linear law on the elastic strain, \(\boldsymbol{\sigma} = \mathbb{C} : \boldsymbol{\varepsilon}^{e}\), with \(\mathbb{C}\) built from elastic_E, elastic_nu and elastic_G. Three constants rather than two is what makes it cubic: an isotropic material satisfies \(G = E / 2(1+\nu)\), and the departure from that identity is the elastic anisotropy ratio, which for most metals is far from 1 and matters at grain boundaries.

Each slip system \(\alpha\) has a unit slip direction \(\mathbf{d}^{\alpha}\) and a unit plane normal \(\mathbf{n}^{\alpha}\), both fixed in the crystal frame. The stress the system feels is the Cauchy stress projected onto the Schmid tensor, after that tensor has been rotated into the sample frame by the orientation matrix \(Q\) (ResolvedShear):

(23)#\[\tau^{\alpha} = \boldsymbol{\sigma} : \Bigl[\, Q \operatorname{sym}\!\bigl(\mathbf{d}^{\alpha} \otimes \mathbf{n}^{\alpha}\bigr) Q^{\mathsf T} \Bigr].\]

The \(Q \,\cdot\, Q^{\mathsf T}\) sandwich is the whole coupling between microstructure and mechanics. Two grains under identical macroscopic stress see different \(\tau^{\alpha}\) purely because \(Q\) differs, so which systems yield, how much they slip, and how the grain hardens are all decided by orientation. Everything the reconstruction contributes to the simulation enters through Eq. 23, and every orientation convention error in the pipeline shows up here as a grain deforming the wrong way.

Flow rule#

The shear rate on system \(\alpha\) follows a power law (PowerLawSlipRule):

(24)#\[ \dot{\gamma}^{\alpha} = \dot{\gamma}_0 \, \left| \frac{\tau^{\alpha}}{\hat{\tau}^{\alpha}} \right|^{\,n-1} \frac{\tau^{\alpha}}{\hat{\tau}^{\alpha}},\]

with reference rate \(\dot{\gamma}_0\) (power_slip_g0), rate exponent \(n\) (power_slip_n), and slip resistance \(\hat{\tau}^{\alpha}\). The \(|\cdot|^{n-1}(\cdot)\) form is how NEML2 writes it and is what the code evaluates; for positive \(\hat{\tau}^{\alpha}\) it equals the more familiar \(\dot{\gamma}_0 |\tau^{\alpha}/\hat{\tau}^{\alpha}|^{n} \operatorname{sgn}(\tau^{\alpha})\), but it is written this way so the derivative stays continuous through \(\tau^{\alpha} = 0\) instead of needing a sign function.

Eq. 24 is a viscoplastic regularisation of a yield criterion, and that is the point of it. Rate-independent crystal plasticity has to decide which subset of slip systems is active, a combinatorial problem with a non-unique answer for a cubic crystal (the Taylor ambiguity). The power law sidesteps it: every system always slips a little, and as \(n \to \infty\) the response approaches the rate-independent limit. So \(n\) is not a fitted material property so much as a numerical choice.

Important

\(n\) controls the stiffness of the integration directly. The Jacobian of Eq. 24 scales as \(n/\hat{\tau}\), so raising \(n\) for a sharper yield tightens the Newton solve at every quadrature point. power_slip_n = 20 in the shipped deck and 25 in the documented recipe are already firm; going much above that trades convergence for a yield corner you probably cannot measure.

Hardening#

The deck uses SingleSlipStrengthMap, so all systems share one scalar resistance:

(25)#\[\hat{\tau}^{\alpha} = \bar{\tau} + \tau_0 \quad \text{for every } \alpha ,\]

with \(\tau_0\) = slip_constant_strength a fixed offset and \(\bar{\tau}\) the single evolving hardening variable. There is no latent hardening and no system-to-system interaction matrix; slip on one system hardens all of them equally. That is a real modelling simplification, and it is why this model reproduces a macroscopic stress-strain curve well while under-predicting the texture-dependent spread between grains.

\(\bar{\tau}\) follows a Voce law (VoceSingleSlipHardeningRule) driven by the total slip rate \(\dot{\gamma}_{\mathrm{tot}} = \sum_{\alpha} |\dot{\gamma}^{\alpha}|\):

(26)#\[\dot{\bar{\tau}} = \theta_0 \left( 1 - \frac{\bar{\tau}}{\tau_f} \right) \dot{\gamma}_{\mathrm{tot}},\]

with initial slope \(\theta_0\) = voce_hardening_initial_slope and saturation \(\tau_f\) = voce_hardening_saturation.

Warning

\(\tau_f\) is the saturated value of \(\bar{\tau}\) alone, and \(\tau_0\) is added on top of it by Eq. 25. The slip resistance therefore saturates at \(\tau_0 + \tau_f\), not at \(\tau_f\). With the documented recipe’s \(\tau_0 = 100\) and \(\tau_f = 220\) MPa the resistance runs from 100 to 320 MPa. Setting voce_hardening_saturation to a measured saturation stress overshoots by \(\tau_0\). (NEML2’s own parameter description for saturated_hardening calls it “the final, saturated value of the slip system strength”, which is loose: it is the saturated value of the evolving variable, before the constant offset.)

Kinematics and reorientation#

Slip produces both a stretching and a spin, obtained from the symmetric and skew parts of the same rotated Schmid tensor (PlasticDeformationRate and PlasticVorticity):

(27)#\[d^{p} = \sum_{\alpha} \dot{\gamma}^{\alpha}\, Q \operatorname{sym}\!\bigl(\mathbf{d}^{\alpha} \otimes \mathbf{n}^{\alpha}\bigr) Q^{\mathsf T}, \qquad w^{p} = \sum_{\alpha} \dot{\gamma}^{\alpha}\, Q \operatorname{skew}\!\bigl(\mathbf{d}^{\alpha} \otimes \mathbf{n}^{\alpha}\bigr) Q^{\mathsf T}.\]

The elastic strain and the lattice orientation then evolve as (ElasticStrainRate, OrientationRate)

(28)#\[\dot{\varepsilon} = d - d^{p} - \varepsilon w + w \varepsilon, \qquad \Omega^{e} \equiv \dot{Q} Q^{\mathsf T} = w - w^{p} - \varepsilon d^{p} + d^{p} \varepsilon,\]

where \(d\) and \(w\) are the total deformation rate and vorticity and \(\varepsilon\) is the elastic strain. The leading \(w - w^{p}\) is the familiar statement that the lattice takes up whatever spin the slip systems do not; the \(\varepsilon d^{p}\) pair is the finite-elastic-strain correction, small for metals but not dropped.

That first term is the mechanism behind texture evolution and intragranular fragmentation. A grain whose active systems produce a plastic spin that does not match the imposed one must rotate its lattice to compensate, and neighbouring regions of the same grain with different active sets rotate differently. Grain segmentation and the reorientation REI criterion both measure the result.

The state variables are integrated implicitly: backward Euler for the elastic strain and the slip hardening, and an implicit exponential map for the orientation. The exponential update is not a refinement – integrating \(\dot{Q} = \Omega^{e} Q\) additively would take \(Q\) off \(SO(3)\) and accumulate a spurious stretch over a long load history.

Note

The orientation state is stored and output as NEML2 v3 MRP, \(\tan(\theta/4)\,\hat{\mathbf{e}}\). The output columns are named ori_rodrigues_{x,y,z}, which is a misnomer: they do not hold classical Rodrigues (Gibbs) vectors, \(\tan(\theta/2)\,\hat{\mathbf{e}}\). Every in-repository consumer treats them correctly as MRP; a third-party tool expecting Rodrigues will silently misinterpret them, and the error grows with rotation angle. See ori_rodrigues_* holds MRP, not Rodrigues.

MOOSE/PUMA assembles these point-wise responses into the global finite-element equilibrium solve over the mesh, advancing in time as the load ramps between initialize_time and total_time. graintrace sets up the model, parameters, mesh, and boundary conditions; NEML2 evaluates the constitutive model and MOOSE/PUMA performs the FE solve.

Algorithm#

  1. Load the mesh and per-block MRP orientations (and optional initial elastic-strain field).

  2. Write the MOOSE input decks (main run, initial conditions, transfer, grid output) from the templates in cpfe_base/.

  3. Bake the material parameters into the NEML2 model and AOTI-compile it with neml2-compile for the target device(s) (recompile rebuilds the .pt2 when parameters change).

  4. Launch puma-opt under the configured launcher (mpiexec/srun); multi-GPU runs map to MPI ranks over a device list.

  5. Collect outputs: the native-mesh Exodus, per-element CSVs (mesh_out/), and, if requested, regular-grid CSVs (grid_out/) at the configured sync/step frequency.

Parameters that matter#

See Configuration reference for the full list.

  • material: slip_constant_strength, voce_hardening_initial_slope, voce_hardening_saturation, power_slip_n, power_slip_g0, elastic_E / elastic_nu / elastic_G, burger_scale.

  • simulation_parameters: dt, total_time, initialize_time, sync_times, device / device_batch, and the output-frequency knobs grid_transfer / exodus_output / mesh_csv.

  • boundary: bounding_box and the bc dict (displacement or stress_free per face).

  • grid_properties: number_of_elements and the inset bounding_box for regular-grid output.

Further details#

For the full constitutive-model formulation, see the NEML2 documentation: https://applied-material-modeling.github.io/neml2/ . For the finite-element framework and the nonlinear/transient solve, see the MOOSE documentation: https://mooseframework.inl.gov/ .

See also#