SlipSystemElasticInteraction

Source: models/solid_mechanics/crystal_plasticity/SlipSystemElasticInteraction.py

The coupling matrix \(A_{ij} = \Delta t\, M_i : \mathbb{C} : M_j\).

Slip on system \(j\) relaxes the resolved shear on system \(i\) through the elastic response: with \(\varepsilon^e = \varepsilon^{e,\rm trial} - \Delta t \sum_j \dot\gamma_j M_j\), the resolved shears follow \(\tau_i = \tau^{\rm trial}_i - \sum_j A_{ij}\dot\gamma_j\). That \(A\) is what turns the slip rule into the condensed system \(\varphi(\dot\gamma) + A\dot\gamma = b\) that CoordinateDescentPredictor solves, with \(b = \tau^{\rm trial}\).

\(A\) is symmetric positive semi-definite – it is \(\Delta t\,B^{\mathsf T}\mathbb{C}B\) with \(B\) the Schmid map and \(\mathbb{C}\) the elastic stiffness – and singular, because the Schmid map sends \(\mathbb{R}^{n_{\rm slip}}\) into a lower-dimensional space of symmetric tensors (the classical Taylor ambiguity). Both properties are relied on downstream: positive semi-definiteness is what brackets each coordinate solve, and singularity is harmless because the slip rule’s own monotone nonlinearity is what makes the condensed system strictly convex.

It is emitted as a Scalar carried on two sub-batch axes, (n_slip, n_slip) – a matrix of scalars is exactly what the sub-batch machinery represents.

Frames. Everything is contracted in the lab frame, matching the rest of the crystal-plasticity library: ResolvedShear rotates the Schmid tensors and contracts them with the lab-frame stress, and GeneralElasticity rotates the stiffness (T_rot = rotate(T, R)) before applying it to the lab-frame strain. So both Schmid operands and the stiffness are rotated here. Getting this wrong yields a plausible but silently wrong matrix, so it is pinned by a test: for an isotropic stiffness, \(A\) must not depend on the orientation at all.

The spin convection is deliberately omitted. ElasticStrainRate carries \(\dot\varepsilon^e = d - d^p + \Omega[\varepsilon^e]\) with \(\Omega[V] = [W, V]\) the material spin, so the elastic-strain block of the Jacobian is \(I - \Delta t\,\Omega\) rather than \(I\) and the exact condensed coupling is \(\Delta t\,M^{\mathsf T}\mathbb{C}(I - \Delta t\,\Omega)^{-1}M\). Because \(\Omega\) generates a rotation it is skew-adjoint, which splits that exactly: the symmetric part is what is computed here to within \(O((\Delta t\lVert\Omega\rVert)^2)\), and the whole deviation is a skew contribution of size \(\Delta t\lVert\Omega\rVert = 2\Delta t\lvert w\rvert\) – the incremental rotation over the step.

That omission is required rather than tolerated. Coordinate descent converges because \(\varphi(\dot\gamma) + A\dot\gamma - b\) is the gradient of a strictly convex potential, which needs \(A\) symmetric; the exact coupling is not symmetric and so is not the gradient of anything. This model computes the variational part and drops the rest. One corollary worth knowing: a skew matrix has zero diagonal, so the \(A_{ii}\ge 0\) precondition that brackets each coordinate solve holds exactly, not to within the spin error.

The \(O((\Delta t\lVert\Omega\rVert)^2)\) bound on the symmetric part assumes \([\mathbb{C}, \Omega] = 0\), i.e. isotropic elasticity. With an anisotropic stiffness the symmetric part also picks up \(\tfrac{1}{2}\Delta t^2 M^{\mathsf T}[\mathbb{C},\Omega]M\), one order lower. Still small, still \(O(\Delta t)\) in the step.

Inputs

orientation — input · MRP · required

Active-convention orientation

time — input · Scalar · default t

Time

Outputs

coupling — output · Scalar · required

The coupling matrix, sub-batched over (n_slip, n_slip)

Parameters

elastic_stiffness_tensor — parameter · SSR4 · required

Elastic stiffness tensor, in the crystal frame

Other options

crystal_geometry — str · default crystal_geometry

The data object carrying the crystallographic information