A dual-length-scale gradient crystal plasticity framework: Finite-element implementation and validation against single-crystal copper micro-bending experiments

Abstract

This paper develops a gradient-enhanced finite-strain crystal plasticity framework with two distinct intrinsic length scales associated with lattice incompatibility and higher-order orientation gradients. The formulation is built on a curvature-based description of elastic rotation and its associated incompatibility measures, and it provides a clear separation between dislocation-driven mechanisms and curvature-gradient (disclination-type) mechanisms within a unified kinematic setting. A differential-geometric interpretation is used to connect these measures to torsion- and curvature-type incompatibilities, thereby clarifying the physical meaning of the additional fields and the origin of size effects. The constitutive structure is posed within a thermodynamically consistent internal-variable framework, yielding coupled relations for classical stress and higher-order stress measures, along with the corresponding balance laws at multiple scales. Analytical validation is established through a set of canonical boundary-value problems. A benchmark single-crystal strip in simple shear demonstrates size-dependent hardening under constrained deformation and reveals boundary-layer formation. Complementary analytical studies treat torsion of a cylindrical bar, bending of a crystal plate, and a bicrystal tilt boundary, providing closed-form reductions that highlight how the two length scales govern the localization of orientation gradients and the spatial organization of defect-related fields. A finite-element implementation is then presented, including an integration-point constitutive update and a mixed variational formulation tailored to the higher-order structure of the theory. Numerical simulations of cantilever micro-bending quantify thickness-dependent strengthening and distinguish the roles of curvature-amplitude strengthening versus curvature-gradient regularization. Finally, the computed micro-bending trends are compared with reported single-crystal copper micro-beam experiments, demonstrating that the proposed two-length framework captures the primary thickness dependence while offering a unified continuum route to incorporate defect geometry into predictive crystal plasticity simulations.

Introduction

Crystal plasticity modeling is a cornerstone of modern solid mechanics and materials science, offering critical insights into the mechanical behavior of crystalline materials under complex loading conditions. The primary challenge lies in developing constitutive frameworks that accurately capture the anisotropic response driven by crystallographic slip while also accounting for microstructural phenomena such as dislocation interactions, lattice curvature, and subgrain formation. Classical crystal plasticity models, rooted in the multiplicative decomposition of the deformation gradient $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p$ , ¹³ have been highly effective in describing macroscopic plastic flow and slip system evolution. However, these models are inherently local and lack intrinsic length scales, which limits their ability to predict size-dependent phenomena observed in experiments such as nanoindentation, grain boundary strengthening, and micro-bending.^19,25

The absence of length scales in classical models stems from their focus on homogeneous deformation fields, neglecting the spatial gradients associated with microstructural defects like geometrically necessary dislocations (GNDs). GNDs arise to accommodate lattice incompatibilities caused by non-uniform plastic deformation, significantly influencing mechanical responses at small scales.^3,35 For instance, phenomena such as the Hall-Petch effect, where yield strength increases with decreasing grain size, or the enhanced hardening observed in micro-scale torsion, underscore the need for models that incorporate gradient effects.^10,11 To address these limitations, strain gradient plasticity (SGP) theories have emerged, introducing higher-order stresses and energetic length scales to capture the effects of strain gradients and defect density.^19,27

The concept of strain-gradient plasticity (SGP) was first introduced in the pioneering works of Aifantis,^1,2 who emphasized the role of internal length scales in regularizing plasticity models. Later developments, such as Fleck and Hutchinson, ¹⁹ extended these ideas into higher-order continuum frameworks, demonstrating that strain gradients introduce size-dependent hardening through additional energetic contributions. Gurtin and Anand ²⁷ and Gurtin ²⁶ extended these concepts to single-crystal plasticity, incorporating GNDs within a microforce balance framework that accounts for the work conjugate to plastic strain gradients. Similarly, Acharya and Bassani ³ proposed a theory where the incompatibility of the plastic distortion field directly informs the dislocation content, linking lattice curvature to defect density. From a computational perspective, prior work has demonstrated that complex constitutive systems can be implemented in robust finite-element settings once the internal variables and their evolution equations are cast in a thermodynamically consistent form; see, for example, the implementation-oriented development of evolving internal state variables in He et al. ³¹ These advancements highlight the importance of embedding microstructural length scales into constitutive models to capture phenomena driven by non-local effects.

Among gradient-enhanced models, Bammann’s contributions^8,9 are particularly notable for introducing a natural length scale through a curvature tensor defined as ${\kappa }_e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T$ , where ${\mathbf{R}}^e$ is the elastic rotation tensor. This curvature tensor, inspired by Nye’s dislocation density tensor, ³⁵ quantifies lattice incompatibility and provides a regularized internal backstress that enhances predictive accuracy in scenarios involving rotational gradients. ⁹ However, Bammann’s formulation primarily addresses first-order incompatibilities associated with GNDs, neglecting higher-order effects such as disclinations—topological defects arising from gradients of the curvature tensor. Disclinations are critical in materials undergoing significant rotations or twinning, where they contribute to rotational non-uniformity and complex defect interactions.^17,24

To overcome these shortcomings, micromorphic and higher-order theories have been proposed. Eringen and Claus ¹⁷ introduced a micromorphic approach that accounts for additional degrees of freedom associated with micro-rotations, providing a framework to model disclination effects. Forest and Sievert ²³ and Forest ²⁴ further developed micromorphic crystal plasticity models, incorporating lattice curvature and micro-deformation fields to capture size effects and defect evolution. Ohno and Okumura ³⁶ and Steinmann ³⁸ explored higher-order stress contributions, demonstrating their role in regularizing strain localization and enhancing hardening predictions. These efforts underscore the need for a unified framework that integrates both curvature and its gradients to fully describe the geometric complexity of defect fields.

The present work proposes a novel crystal plasticity model that systematically incorporates two internal length scales associated with lattice curvature and its gradient. Building on Bammann’s curvature-based approach, ⁹ the framework is extended to include higher-order effects by introducing the gradient of the curvature tensor, $\mathrm{Curl}\left({\kappa }_e\right)$ , which captures disclination-driven incompatibilities. ¹⁷ The model employs a generalized Helmholtz free energy that includes quadratic contributions from elastic strain, curvature, and curvature gradients, ensuring thermodynamic consistency. ⁶ This approach aligns with the geometric perspective of Acharya and Bassani ³ and the micromorphic formulations of Forest, ²⁴ embedding GNDs and disclinations within a Riemann–Cartan framework where dislocations correspond to torsion and disclinations to curvature. Crystal plasticity has also been coupled to microstructure-evolution models; see, for example, the CP-phase-field integration for process–structure–property prediction in additively manufactured metals in Liu ³⁴ The present contribution is not the existence of such couplings, but a two-length-scale, defect-geometric CP core that energetically distinguishes torsion- and curvature-driven incompatibilities, enabling more interpretable control of strengthening and localization regularization in coupled settings.

By incorporating both ${\kappa }_e$ and $\mathrm{Curl}\left({\kappa }_e\right)$ as internal variables, the model introduces distinct length scales that penalize rotational incompatibility and non-uniformity, respectively. These length scales regularize strain localization and enable the prediction of size-dependent hardening, addressing limitations of classical models in micro-scale applications. ¹⁸ The resulting higher-order balance laws, inspired by Anand et al., ⁶ include classical stresses and higher-order moment stresses conjugate to curvature gradients, providing a robust framework for analytical validation. ⁴ This work bridges mesoscopic defect mechanics with macroscopic continuum theory, offering a comprehensive tool for studying complex deformation phenomena in crystalline materials.

Recent developments in strain-gradient and higher-order plasticity theories have highlighted the importance of internal length scales in capturing size-dependent phenomena such as strain localization, strengthening in small-scale crystals, and dislocation patterning.^1,19,24,26 While conventional formulations often rely on phenomenological gradients of plastic strain, our approach embeds these effects into a rigorous geometric framework grounded in differential geometry.

Specifically, the model is formulated on a Riemann–Cartan manifold, where dislocations and disclinations arise naturally as torsion and curvature of the material connection. The elastic rotation ${\mathbf{R}}^e\in \mathrm{SO}\left(3\right)$ defines a local crystallographic frame, and its incompatibilities are quantified through the dislocation density tensor ${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^{\top }$ and the disclination density tensor $\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right)$ . These measures provide a direct link between continuum kinematics and microscopic defect structures, extending the classical Kröner–Bilby theory.^12,16

The dual internal length scales of the present model acquire a clear physical interpretation: ${\ell }_1$ penalizes torsion-driven incompatibilities associated with dislocations, while ${\ell }_2$ penalizes curvature-driven incompatibilities associated with disclinations. This separation is a central added value of the formulation: rather than collapsing distinct size-effect mechanisms into a single gradient parameter, it enables one to disentangle dislocation-mediated strengthening from curvature-mediated regularization, thereby providing direct control over both size-dependent hardening and the thickness of boundary layers.

Beyond interpretation, the two-length-scale structure has concrete computational benefits. First, the energetic embedding of torsion and curvature terms yields well-defined generalized stresses and higher-order boundary conditions, which clarifies what must be prescribed (or left natural) in boundary-value problems and supports a systematic weak formulation suitable for finite-element discretization. Second, the distinct numerical roles of ${\ell }_1$ and ${\ell }_2$ improve robustness and identifiability: ${\ell }_1$ primarily governs strengthening through torsion-related defect storage, while ${\ell }_2$ primarily governs field smoothness and localization control through curvature-related incompatibilities. As a result, the model can regularize mesh-sensitive localization in a physically interpretable way and enables targeted calibration strategies (e.g. using bending/torsion size-effect data to identify ${\ell }_1$ and boundary-layer or gradient-dominated responses to identify ${\ell }_2$ ). Third, the framework is implementation-ready for multiscale coupling because it provides a consistent defect-geometric “interface”: curvature and its gradient appear as energetically conjugate measures that can be exchanged with phase-field microstructure variables or discrete dislocation dynamics inputs without introducing ad hoc smoothing.

In this way, the model provides not only a higher-order energetic extension of multiplicative plasticity, but also a genuine differential-geometric foundation. Section 3 develops this framework in detail, explicitly introducing connections on $\mathrm{SO}\left(3\right)$ , torsion, curvature, and their energetic embedding into the free energy, while the subsequent sections detail the associated weak form, finite-element implementation, and illustrative numerical benchmarks.

The remainder of this paper is structured as follows:

• Kinematics: This section establishes the kinematic framework, detailing the multiplicative decomposition of the deformation gradient, polar decomposition, velocity gradient, and definitions of the curvature tensor ${\kappa }_e$ and disclination density tensor $\theta$ , introducing two internal length scales.

• Geometric interpretation: A differential geometric framework is developed, linking curvature and its gradient to torsion and curvature in a Riemann–Cartan manifold, providing a rigorous description of defect-induced incompatibilities.

• Thermodynamics: A thermodynamically consistent model is formulated, incorporating internal variables into a generalized Helmholtz free energy and deriving constitutive relations for stresses and higher-order moment stresses.

• Crystal plasticity model: The model couples backstress and moment stress to crystallographic slip evolution, using a power-law flow rule and specifying the evolution of statistically stored dislocations.

• Balance laws at multiple scales: Mesoscopic and macroscopic balance laws are derived to account for internal stress and couple stress fields, ensuring regularization of strain localization.

• Analytical validation: Beyond simple shear of a single-crystal strip—where the framework is compared against classical single-slip formulations and shown to reproduce size-dependent hardening consistent with established theories—analytical reductions are presented for bending of crystal plates, torsion of cylindrical bars, and bicrystal tilt boundaries. These examples highlight the distinct roles of the two internal length scales: ${\ell }_1$ governs the backstress associated with GNDs, while ${\ell }_2$ controls boundary-layer formation and disclination localization.

• Finite element simulations and comparison with experiments: This section details the finite-element formulation and the local constitutive update used for the computational implementation. The same section presents some numerical simulations, with emphasis on cantilever micro-bending size effects and comparisons with representative experimental trends.

• Conclusion: The model’s contributions are summarized, emphasizing its thermodynamic and geometric rigor and outlining potential for future additional analytical or numerical investigations.

Kinematics

The formulation of finite-strain crystal plasticity relies on a precise kinematic framework to capture large deformations, rotations, and their gradients, which are essential for constitutive modeling and numerical implementation. ¹³ Unlike infinitesimal plasticity theories limited to small displacements, finite-strain models accommodate arbitrary deformations and microstructural incompatibilities, such as those induced by dislocations and disclinations. ²⁶ This section establishes the kinematic foundation for a gradient-enhanced crystal plasticity model incorporating two internal length scales, associated with lattice curvature and its gradient, to address size-dependent phenomena.

Multiplicative decomposition of the deformation gradient

The kinematic description begins with the multiplicative decomposition of the deformation gradient $\mathbf{F}$ , which maps material vectors from the reference configuration ${\mathcal{B}}_0$ to the current configuration ${\mathcal{B}}_t$ :

$$\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p,$$

(1)

where ${\mathbf{F}}^p$ represents the plastic deformation, preserving the lattice structure and volume, and ${\mathbf{F}}^e$ accounts for elastic distortions and lattice rotations. ¹³ This decomposition defines an intermediate configuration ${\mathcal{B}}_i$ , which is stress-free but generally incompatible, as ${\mathbf{F}}^p$ does not satisfy integrability conditions due to plastic distortion. ³

Before introducing the curvature and disclination measures, it is useful to recall the role of the intermediate configuration in the multiplicative decomposition and to specify where Curl operations are consistently evaluated.

The intermediate configuration ${\mathcal{B}}_i$ is defined through the multiplicative split $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p$ , where ${\mathbf{F}}^p$ maps the reference configuration ${\mathcal{B}}_0$ to ${\mathcal{B}}_i$ , and ${\mathbf{F}}^e$ maps ${\mathcal{B}}_i$ to the current configuration ${\mathcal{B}}_t$ . This configuration is stress-free but generally incompatible, since ${\mathbf{F}}^p$ need not satisfy integrability conditions. Following Mandel’s isoclinic framework, all Curl operations that define the curvature tensor:

$${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}_e\right){\mathbf{R}}_e^{\top }$$

and the disclination density:

$$\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right)$$

are consistently evaluated in ${\mathcal{B}}_i$ . In particular, slip-system vectors ${\mathit{m}}^i$ and ${\mathit{s}}^i$ used in equation (32) are defined in this intermediate configuration. This ensures frame invariance and coherence of the higher-order measures throughout the formulation.

To avoid any ambiguity between different classes of variables, the present framework explicitly distinguishes between the degrees of freedom, state variables, and internal variables. The degrees of freedom are the primary kinematic fields, such as the deformation gradient and the elastic rotation, which evolve from balance laws. The state variables are those that enter the Helmholtz free energy, namely the elastic strain, curvature tensor, and disclination density, which define reversible energetic contributions. Finally, the internal variables represent history-dependent quantities, such as the statistically stored dislocation density, that evolve according to phenomenological laws and capture dissipative microstructural mechanisms. Table 1 summarizes this taxonomy, ensuring clarity in the logical structure of the model.

Table 1.

Classification of variables in the dual-length-scale crystal plasticity model.

Category	Variables	Role in the model
DOFs	• Deformation gradient $\mathbf{F}$ • Elastic rotation ${\mathbf{R}}_e$ (from ${\mathbf{F}}^e$ in the multiplicative split)	Primary kinematic fields governed by balance laws (Section 6). They evolve from equilibrium and compatibility conditions
State variables (entering free energy)	• Elastic strain ${\mathbf{E}}^e$ • Curvature tensor ${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}_e\right){\mathbf{R}}_e^T$ • Disclination density $\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right)$	Energetic quantities entering the Helmholtz free energy (Section 4.4). They represent reversible storage of elastic, curvature, and curvature-gradient energies
Internal variables (history-dependent)	• SSD density ${\mathbf{\rho }}_{\mathrm{SS}}$ • Equivalent measure ${\epsilon }_{\mathrm{SS}}=b\sqrt{{\mathbf{\rho }}_{\mathrm{SS}}}$	History-dependent variables describing irreversible microstructural evolution and hardening. They obey phenomenological evolution laws (Section 5)

DOFs: degrees of freedom; SSD: statistically stored dislocation.

Polar decomposition and velocity gradient

To incorporate lattice rotations, the elastic deformation gradient ${\mathbf{F}}^e$ is expressed via polar decomposition:

$${\mathbf{F}}^e={\mathbf{V}}^e{\mathbf{R}}^e,$$

(2)

where ${\mathbf{V}}^e$ is the symmetric, positive-definite elastic left stretch tensor, and ${\mathbf{R}}^e\in \mathrm{SO}\left(3\right)$ is the orthogonal rotation tensor. ²⁷ Combining equations (2) with (1), the total deformation gradient is:

$$\mathbf{F}={\mathbf{V}}^e{\mathbf{R}}^e{\mathbf{F}}^p.$$

(3)

The velocity gradient $\mathbf{L}=\stackrel{·}{\mathbf{F}}{\mathbf{F}}^{-1}$ describes the rate of deformation. Differentiating equation (3) and using the inverse ${\mathbf{F}}^{-1}={\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^{\top }{\left({\mathbf{V}}^e\right)}^{-1}$ , yields:

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{L}=\stackrel{·}{\mathbf{F}}{\mathbf{F}}^{-1}}\\ \multicolumn{1}{c}{={\stackrel{·}{\mathbf{V}}}^e{\left({\mathbf{V}}^e\right)}^{-1}+{\mathbf{V}}^e{\stackrel{·}{\mathbf{R}}}^e{\left({\mathbf{R}}^e\right)}^{\top }{\left({\mathbf{V}}^e\right)}^{-1}}\\ \multicolumn{1}{c}{+{\mathbf{V}}^e{\mathbf{R}}^e{\stackrel{·}{\mathbf{F}}}^p{\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^{\top }{\left({\mathbf{V}}^e\right)}^{-1}}\\ \multicolumn{1}{c}{={\mathbf{L}}^e+{\mathbf{\omega }}^e+{\mathbf{L}}^p,}\end{array}$$

(4)

with:

• ${\mathbf{L}}^e={\stackrel{·}{\mathbf{V}}}^e{\left({\mathbf{V}}^e\right)}^{-1}$ , the elastic stretching rate,

• ${\mathbf{\omega }}^e={\mathbf{V}}^e{\stackrel{·}{\mathbf{R}}}^e{\left({\mathbf{R}}^e\right)}^{\top }{\left({\mathbf{V}}^e\right)}^{-1}$ , the elastic spin,

• ${\mathbf{L}}^p={\mathbf{V}}^e{\mathbf{R}}^e{\stackrel{·}{\mathbf{F}}}^p{\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^{\top }{\left({\mathbf{V}}^e\right)}^{-1}$ , the plastic flow contribution.

Equation (4) decomposes the deformation rate into elastic, rotational, and plastic components, consistent with finite-strain kinematics. ⁸

To analyze these quantities in the intermediate configuration ${\mathcal{B}}_i$ , the intermediate velocity gradient is defined as:

$$\tilde{\mathbf{L}}={\left({\mathbf{V}}^e\right)}^{-1}\mathbf{L}{\mathbf{V}}^e.$$

A transformation of equation (4) yields:

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\mathbf{L}}={\left({\mathbf{V}}^e\right)}^{-1}{\stackrel{·}{\mathbf{V}}}^e+{\stackrel{·}{\mathbf{R}}}^e{\left({\mathbf{R}}^e\right)}^{\top }+{\mathbf{R}}^e{\stackrel{·}{\mathbf{F}}}^p{\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^{\top }}\\ \multicolumn{1}{c}{={\tilde{\mathbf{L}}}^e+{\tilde{\mathbf{\Omega }}}^e+{\tilde{\mathbf{L}}}^p,}\end{array}$$

(5)

where:

• ${\tilde{\mathbf{L}}}^e={\left({\mathbf{V}}^e\right)}^{-1}{\stackrel{·}{\mathbf{V}}}^e$ , the elastic velocity gradient,

• ${\tilde{\mathbf{\Omega }}}^e={\stackrel{·}{\mathbf{R}}}^e{\left({\mathbf{R}}^e\right)}^{\top }$ , the lattice spin,

• ${\tilde{\mathbf{L}}}^p={\mathbf{R}}^e{\stackrel{·}{\mathbf{F}}}^p{\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^{\top }$ , the plastic velocity gradient.

The components are split into symmetric and skew-symmetric parts:

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathbf{L}}}^e={\tilde{\mathbf{D}}}^e+{\tilde{\mathbf{W}}}^e,{\tilde{\mathbf{D}}}^e=\frac{1}{2}({\left({\mathbf{V}}^e\right)}^{-1}{\stackrel{·}{\mathbf{V}}}^e+{\stackrel{·}{\mathbf{V}}}^e{\left({\mathbf{V}}^e\right)}^{-1}),}\\ \multicolumn{1}{c}{{\tilde{\mathbf{W}}}^e=\frac{1}{2}({\left({\mathbf{V}}^e\right)}^{-1}{\stackrel{·}{\mathbf{V}}}^e-{\stackrel{·}{\mathbf{V}}}^e{\left({\mathbf{V}}^e\right)}^{-1})+{\tilde{\mathbf{\Omega }}}^e,}\\ \multicolumn{1}{c}{{\tilde{\mathbf{L}}}^p={\tilde{\mathbf{D}}}^p+{\tilde{\mathbf{W}}}^p,{\tilde{\mathbf{D}}}^p=\frac{1}{2}({\tilde{\mathbf{L}}}^p+{\left({\tilde{\mathbf{L}}}^p\right)}^{\top }),}\\ \multicolumn{1}{c}{{\tilde{\mathbf{W}}}^p=\frac{1}{2}({\tilde{\mathbf{L}}}^p-{\left({\tilde{\mathbf{L}}}^p\right)}^{\top }),}\end{array}$$

(6)

where ${\tilde{\mathbf{D}}}^e$ and ${\tilde{\mathbf{D}}}^p$ are the elastic and plastic rate-of-deformation tensors, and ${\tilde{\mathbf{W}}}^e$ and ${\tilde{\mathbf{W}}}^p$ are the vorticity tensors. The lattice spin ${\tilde{\mathbf{\Omega }}}^e$ drives lattice rotation, while ${\tilde{\mathbf{W}}}^p$ represents material spin relative to the lattice. ²⁴ For small elastic strains, ${\tilde{\mathbf{L}}}^e\approx 0$ , simplifying computations. ⁹

Curvature tensor and geometrically necessary dislocations

Lattice incompatibility due to geometrically necessary dislocations (GNDs) is quantified by the curvature tensor:

$${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^{\top },$$

(7)

where $\mathrm{Curl}\left({\mathbf{R}}^e\right)$ is computed row-wise:

$$[\mathrm{Curl}\left({\mathbf{R}}^e\right){]}_{\mathrm{ij}}={ϵ}_{\mathrm{ikl}}{\partial }_k{\left({\mathbf{R}}^e\right)}_{\mathrm{lj}},$$

with ${ϵ}_{\mathrm{ikl}}$ as the Levi–Civita symbol. The conjugation by ${\left({\mathbf{R}}^e\right)}^{\top }$ ensures frame invariance, and ${\mathbf{\kappa }}^e$ aligns with Nye’s dislocation density tensor, capturing GNDs from rotational incompatibilities.^9,35 This tensor introduces an internal length scale, as its magnitude scales inversely with the characteristic length of lattice misorientation, influencing backstress and size effects in micro-scale deformation.^3,19

Geometrically necessary dislocations emerge to accommodate nonuniform slip, and their density is represented by the curvature tensor ${\mathbf{\kappa }}^e$ . As deformation proceeds, spatial gradients of the elastic rotation ${\mathbf{R}}_e$ lead to an increase in ${\mathbf{\kappa }}^e$ , corresponding to the buildup of GNDs. These dislocations generate a backstress $\mathbf{\chi }\propto {\ell }_1^2\mu {\mathbf{\kappa }}^e$ that resists slip and contributes to size-dependent hardening. Local concentration of ${\mathbf{\kappa }}^e$ near boundaries also influences the onset of strain localization. This interpretation is consistent with recent discussions of dislocation evolution and backstress effects in crystalline solids.^14,32

Curvature gradient and disclination effects

Higher-order incompatibilities, associated with disclinations, are captured by the disclination density tensor:

$$\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right),$$

(8)

with components:

$${\theta }_{\mathrm{ij}}={ϵ}_{\mathrm{ikl}}{\partial }_k{\left({\mathbf{\kappa }}^e\right)}_{\mathrm{lj}}.$$

This tensor measures the spatial variation of lattice curvature, introducing a second length scale tied to rotational non-uniformity. ¹⁷ Disclinations are critical in materials undergoing large rotations, twinning, or grain boundary evolution, contributing to subgrain formation and enhanced hardening. ²⁴ By including ${\mathbf{\kappa }}^e$ and $\mathit{\theta }$ , the model accounts for the energetic costs of GNDs and disclinations, enabling predictions of complex microstructural evolution.^6,27

These higher-order tensors are not only formal geometric objects but also possess clear experimental significance, as illustrated below.

The curvature tensor ${\mathbf{\kappa }}^e$ introduced here corresponds to measurable lattice rotation gradients and therefore aligns with experimentally estimated GND densities from EBSD or Laue diffraction data. Likewise, the disclination density:

$$\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right)$$

becomes nonzero in regions of sharp misorientation such as grain boundaries or twin interfaces, which can be characterized by EBSD kernel misorientation or Laue diffraction mapping. This correspondence provides a direct bridge between the geometric tensors of the model and experimentally accessible quantities.

Geometric interpretation

The present gradient crystal plasticity model is grounded in a differential-geometric framework where internal variables naturally align with the constructs of a Riemann–Cartan manifold. ³⁹ In such a setting, the elastic rotation ${\mathbf{R}}^e\in \mathrm{SO}\left(3\right)$ defines a moving orthonormal frame, the connection on $\mathrm{SO}\left(3\right)$ governs its variation, torsion corresponds to dislocations, and curvature corresponds to disclinations. This identification provides a precise and self-consistent description of defect-induced incompatibilities that drive size-dependent mechanical responses.^3,12,16,24

The crystal as a differentiable manifold

The reference configuration ${\mathcal{B}}_0$ is considered as a smooth three-dimensional differentiable manifold endowed with a metric induced by the Euclidean structure. In the absence of defects, the crystallographic frame is globally integrable, and the elastic distortion field derives from a compatible deformation map $\mathbf{u}:{\mathcal{B}}_0\to {\mathbb{R}}^3$ ,

$$\mathbf{F}=\mathbf{I}+\nabla \mathbf{u}.$$

(9)

In the presence of defects, global compatibility fails, and the manifold must be described with a nontrivial connection possessing torsion and curvature.

Within the multiplicative decomposition,

$$\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p,{\mathbf{F}}^e={\mathbf{R}}^e{\mathbf{U}}^e,$$

(10)

the plastic part ${\mathbf{F}}^p$ represents an incompatible mapping that carries a defect-free configuration into one with distributed lattice distortions. The elastic rotation ${\mathbf{R}}^e\in \mathrm{SO}\left(3\right)$ describes the local orientation of the lattice frame, while ${\mathbf{U}}^e$ captures the stretch.

Connections on $\mathrm{SO}\left(3\right)$

Let $\{{\mathbf{e}}_i\left(x\right)\}$ denote the orthonormal frame defined by the columns of ${\mathbf{R}}^e\left(x\right)$ . The variation of this frame is captured by the connection one-forms ${\omega }_{\mathrm{ij}}$ defined through:

$$\nabla {\mathbf{e}}_i={\omega }_{\mathrm{ij}}\otimes {\mathbf{e}}_j,{\omega }_{\mathrm{ij}}=-{\omega }_{\mathrm{ji}}.$$

(11)

Equivalently, for the rotation tensor itself,

$${\nabla }_j{\mathbf{R}}^e={\mathbf{\omega }}_j{\mathbf{R}}^e,$$

(12)

where ${\mathbf{\omega }}_j\in \mathfrak{so}\left(3\right)$ encodes local gradients of the lattice orientation. This construction ensures objectivity and frame invariance.

Torsion and dislocations

The torsion tensor is the skew part of the connection acting on the basis vectors,

$$\mathcal{T}(X,Y)={\nabla }_{X}Y-{\nabla }_{Y}X-[X,Y],$$

(13)

with $X,Y$ tangent vector fields and [·,·] their Lie bracket. A nonzero torsion implies closure failure of parallelograms spanned by $X$ and $Y$ , which is physically interpreted as the presence of dislocations carrying a Burgers vector.

In terms of the elastic rotation, torsion reduces to the lattice curvature tensor:

$${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^{\top }.$$

(14)

This tensor coincides with the density of geometrically necessary dislocations (GNDs) introduced in the kinematic formulation. Thus, dislocations are identified with torsion in the Riemann–Cartan framework.

To quantify dislocation content explicitly, one recovers the Kröner–Bilby relation,

$$\mathit{\alpha }=\mathrm{Curl}{\mathbf{F}}^p{\left({\mathbf{F}}^p\right)}^{-1},$$

(15)

where $\mathit{\alpha }$ is the dislocation density tensor. At finite strains, ${\mathbf{\kappa }}^e$ provides an equivalent torsion-based measure consistent with this classical construction.

Curvature and disclinations

The curvature of the connection is obtained from the commutator of covariant derivatives,

$$\mathcal{R}(X,Y)Z={\nabla }_X{\nabla }_{Y}Z-{\nabla }_Y{\nabla }_{X}Z-{\nabla }_{[X,Y]}Z.$$

(16)

A nonzero curvature indicates rotational mismatch of the lattice frame upon parallel transport around a closed loop, corresponding to disclinations with nonzero Frank vector.

In our framework, curvature is captured by the disclination density tensor:

$$\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right).$$

(17)

Thus, torsion ${\mathbf{\kappa }}^e$ and curvature $\mathit{\theta }$ are precisely the continuum measures of dislocations and disclinations, respectively.

Incompatibility in Riemann–Cartan geometry

In a conventional (Riemannian) elastic body, compatibility requires that the deformation gradient $\mathbf{F}$ derives from a smooth displacement field. This condition is equivalent to vanishing torsion and curvature of the associated connection. In linear elasticity, this reduces to the Saint–Venant compatibility condition, which requires the strain tensor $e_{\mathrm{ij}}$ to satisfy:

$${ϵ}_{\mathrm{ikl}}{ϵ}_{\mathrm{jmn}}{\partial }_k{\partial }_m{e}_{\ln }=0,$$

(18)

ensuring that the strain field is compatible with a single-valued displacement field $\mathbf{u}$ .

In a Riemann–Cartan manifold, incompatibility arises naturally from nonzero torsion and curvature:

• The torsion tensor ${\mathbf{\kappa }}^e$ measures the closure failure of infinitesimal Burgers circuits, encoding the density of dislocations.

• The curvature tensor $\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right)$ measures the rotation mismatch upon parallel transport around closed loops, encoding the density of disclinations.

Thus, dislocations and disclinations are not auxiliary constructs, but intrinsic geometric features of the connection on the intermediate configuration.

Compatibility conditions

A defect-free state is recovered when both torsion and curvature vanish:

$${\mathbf{\kappa }}^e=0,\mathit{\theta }=0.$$

(19)

In this limit, ${\mathbf{R}}^e$ becomes globally compatible and derives from a single-valued deformation mapping, consistent with equation (18) in the linearized regime. Conversely, nonzero ${\mathbf{\kappa }}^e$ or $\mathit{\theta }$ signals geometric incompatibility, directly tied to the presence of crystalline defects.

Energetic embedding with dual length scales

The Riemann–Cartan interpretation motivates a natural embedding of incompatibility into the free energy. The dual length scales ${\ell }_1$ and ${\ell }_2$ penalize torsion- and curvature-driven defects, respectively:

$$\psi =\frac{1}{2}\mu \parallel {\mathbf{E}}^e{\parallel }^2+\frac{1}{2}\mu {\ell }_1^2\parallel {\mathbf{\kappa }}^e{\parallel }^2+\frac{1}{2}\mu {\ell }_2^4\parallel \mathit{\theta }{\parallel }^2.$$

(20)

Here,

• ${\ell }_1$ controls the energetic cost of dislocations (torsion), capturing size effects linked to GNDs.

• ${\ell }_2$ controls the energetic cost of disclinations (curvature), capturing the energetic penalty of sharp orientation gradients at boundaries and interfaces.

These terms extend the conventional multiplicative decomposition $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p$ into a fully geometric framework, that is, both thermodynamically admissible and consistent with defect theory. In particular, equation (20) shows how the dual length scales arise as geometric regulators of incompatibility, generalizing Saint–Venant’s classical condition to the nonlinear Riemann–Cartan setting.

The framework of ${\mathbf{\kappa }}^e$ , $\mathit{\alpha }$ , and $\mathit{\theta }$ provides a geometrically rigorous description of microstructural incompatibilities, aligning with the Riemann–Cartan perspective where dislocations and disclinations manifest as torsion and curvature, respectively. ³⁹ This formulation supports the model’s predictive capability for size-dependent phenomena and strain localization, consistent with gradient plasticity theories.^19,26

Thermodynamics

The present section develops a thermodynamically consistent internal variable model for finite-strain crystal plasticity, extending the classical framework of Anand and Gurtin. ⁶ Following the conventional structure of continuum thermodynamics, the Helmholtz free energy is first postulated as a function of the chosen state variables; the internal power is then defined in terms of standard and higher-order stresses; and finally, the dissipation inequality is derived in general form and subsequently reduced using the corresponding energetic conjugates. This sequence ensures that the reduced dissipation inequality emerges naturally from the free-energy postulate and the internal power statement, in accordance with the classical approaches of Coleman–Noll, Germain, and Gurtin.

The formulation is inspired by continuum dislocation theory^3; however, the generalized defect energy is postulated in quadratic form in order to capture both torsion- and curvature-driven mechanisms, incorporating geometrically necessary dislocations (GNDs) and disclinations through higher-order tensorial measures. Defined in the intermediate configuration, the model accounts for lattice curvature and its gradient, thereby introducing two intrinsic length scales associated with dislocation and disclination fields. These length scales enable the prediction of size-dependent phenomena, such as enhanced hardening and strain localization, which are critical for microscale applications.^19,27

Internal variables and curvature measures

To capture the microstructural effects of GNDs and disclinations, the Helmholtz free energy $\tilde{\psi }$ is postulated to depend on the elastic Green–Lagrange strain ${\tilde{\mathbf{E}}}^e=\frac{1}{2}[{\left({\mathbf{F}}^e\right)}^T{\mathbf{F}}^e-\mathbf{I}]$ , temperature $\mathcal{T}$ , a scalar variable ${\epsilon }_{\mathrm{SS}}$ representing statistically stored dislocations (SSDs), the curvature tensor ${\tilde{\mathbf{\kappa }}}^e$ , and its Curl $\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}^e\right)$ . The curvature tensor, which quantifies lattice incompatibility due to GNDs, ³⁵ is defined as:

$${\tilde{\mathbf{\kappa }}}^e=\left(\mathrm{Curl}{\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T,$$

(21)

where ${\mathbf{R}}^e$ is the elastic rotation tensor and $\mathrm{Curl}$ operates row-wise on ${\mathbf{R}}^e$ in the intermediate configuration. ⁹ The Curl of the curvature tensor, $\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}^e\right)$ , provides a higher-order measure associated with disclinations by capturing spatial variations of lattice curvature. ¹⁷

The SSD variable is defined as:

$${\epsilon }_{\mathrm{SS}}=b\sqrt{{\rho }_{\mathrm{SS}}},$$

(22)

where $b$ is the Burgers vector magnitude and ${\rho }_{\mathrm{SS}}$ is the SSD density. This scalar quantifies the isotropic hardening contribution of randomly distributed dislocations, complementing the anisotropic effects of GNDs. ¹¹

Constitutive relations from free energy

Assuming $\tilde{\psi }$ is a smooth function of its arguments, the Coleman–Gurtin method ⁶ is applied to derive the thermodynamic conjugates. For small elastic strains, the elastic strain rate is approximated by ${\stackrel{·}{\stackrel{~}{\mathbf{E}}}}^e\approx {\tilde{\mathbf{D}}}^e$ , where ${\tilde{\mathbf{D}}}^e$ denotes the elastic rate-of-deformation tensor in the intermediate configuration. ²⁶ The entropy inequality equation (44) is satisfied by postulating:

$$\tilde{\mathit{\sigma }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{E}}}^e},$$

(23)

$$\eta =-\frac{\partial \tilde{\psi }}{\partial \mathbb{T}},$$

(24)

$$\tau =\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\epsilon }_{\mathrm{SS}}},$$

(25)

$$\tilde{\mathbf{\chi }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{\kappa }}}_e},$$

(26)

$$\tilde{\mathbf{\Gamma }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial \mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right)},$$

(27)

where $\tilde{\mathit{\sigma }}$ is the stress tensor, $\eta$ is the entropy, $\tau$ is the conjugate force to SSDs, $\tilde{\mathbf{\chi }}$ is the backstress associated with GNDs, and $\tilde{\mathbf{\Gamma }}$ is the higher-order moment stress conjugate to curvature gradients. These relations ensure non-negative mechanical dissipation, consistent with thermodynamic principles. ²⁷

These general Coleman–Noll relations provide the thermodynamic structure of the theory. To obtain explicit constitutive laws, the form of the Helmholtz free energy is specified next.

Free energy specification

Guided by the physical mechanisms of elasticity, dislocation curvature, and disclination density, a quadratic Helmholtz free energy functional is proposed that incorporates energetic contributions from elastic strain, SSDs, curvature, and curvature gradients, motivated by prior gradient plasticity models^19,24:

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\rho }\tilde{\psi }({\tilde{\mathbf{E}}}^e,\mathcal{T},{\epsilon }_{\mathrm{SS}},{\tilde{\mathbf{\kappa }}}_e,\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right))=\mu \left(\mathcal{T}\right)\parallel \mathrm{dev}\left({\tilde{\mathbf{E}}}^e\right){\parallel }^2}\\ \multicolumn{1}{c}{+\frac{K\left(\mathcal{T}\right)}{2}{(\mathrm{tr}{\tilde{\mathbf{E}}}^e)}^2+\frac{1}{2}{c}_{\tau }\mu \left(\mathcal{T}\right){\epsilon }_{\mathrm{SS}}^2+\frac{1}{2}{\ell }_1^2\mu \left(\mathcal{T}\right){\parallel {\tilde{\mathbf{\kappa }}}_e\parallel }^2}\\ +\frac{1}{2}{\ell }_2^4\mu \left(\mathcal{T}\right){\parallel \mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right)\parallel }^2+g\left(\Theta \right),\end{array}$$

(28)

where $\mu \left(\mathcal{T}\right)$ and $K\left(\mathcal{T}\right)$ are the temperature-dependent shear and bulk moduli, $c_{\tau }$ is a dimensionless hardening parameter, ${\ell }_1$ and ${\ell }_2$ are internal length scales associated with GNDs and disclinations, respectively, and $g\left(\mathcal{T}\right)$ is a thermal potential. The deviatoric and volumetric elastic strain terms ensure positive definiteness, while the curvature terms introduce size-dependent effects.^9,17 The norms are defined as $\parallel \mathbf{A}{\parallel }^2=A_{\mathrm{ij}}{A}_{\mathrm{ij}}$ , ensuring frame indifference.

Substituting this free energy into the constitutive framework of Section 4.2 yields explicit expressions for the stresses, presented next.

Resulting constitutive equations

Differentiating the free energy equation (28) with respect to its arguments yields the constitutive relations:

$$\tilde{\mathit{\sigma }}=2\mu \left(\mathcal{T}\right)\mathrm{dev}\left({\tilde{\mathbf{E}}}^e\right)+K\left(\mathcal{T}\right)\mathrm{tr}\left({\tilde{\mathbf{E}}}^e\right)\mathbf{I},$$

(29)

$$\tau =c_{\tau }\mu \left(\mathcal{T}\right){\epsilon }_{\mathrm{SS}},$$

(30)

$$\tilde{\mathbf{\chi }}={\ell }_1^2\mu \left(\mathcal{T}\right){\tilde{\mathbf{\kappa }}}_e,$$

(31)

$$\tilde{\mathbf{\Gamma }}={\ell }_2^4\mu \left(\mathcal{T}\right)\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right),$$

(32)

where $\mathbf{I}$ is the identity tensor. Substituting equation (22) into equation (30) recovers $\tau =c_{\tau }\mu \left(\mathcal{T}\right)b\sqrt{{\rho }_{\mathrm{SS}}}$ , consistent with isotropic hardening models. ¹⁰ The backstress $\tilde{\mathbf{\chi }}$ is proportional to the curvature tensor, reflecting the energetic cost of GNDs, ⁹ while $\tilde{\mathbf{\Gamma }}$ penalizes curvature gradients, capturing disclination effects. ²⁴

Connection to balance laws

The micro- and macro-balance equations presented in this section follow directly from the principle of virtual power applied to standard and higher-order stresses. The explicit derivation is provided here to clarify the role of the curvature tensor ${\mathbf{\kappa }}^e$ and to justify the form of equation (32).

Let $\mathcal{B}$ denote the body (intermediate configuration ${\mathcal{B}}_i$ ), and use the standard differential operators $\mathrm{Grad}$ , $\mathrm{Div}$ , and the row-wise $\mathrm{Curl}$ acting on second-order tensors. Let $\delta \mathit{v}$ be an admissible virtual velocity, and $\delta {\mathit{R}}_e$ a virtual elastic rotation (skew), with the induced virtual measures $\delta {\mathit{E}}^e$ , $\delta {\mathbf{\kappa }}^e$ , and $\delta \mathit{\theta }=\mathrm{Curl}\left(\delta {\mathbf{\kappa }}^e\right)$ .

Internal and external virtual power

$${\mathcal{P}}_{\mathrm{int}}={\int }_{\mathcal{B}}(\mathit{\sigma }:\delta {\mathit{E}}^e+\mathbf{\chi }:\delta {\mathbf{\kappa }}^e+\mathbf{\Gamma }:\delta \mathit{\theta })\mathrm{d}V,$$

(33)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{P}}_{\mathrm{ext}}={\int }_{\mathcal{B}}\mathit{b}·\delta \mathit{v}\mathrm{d}V+{\int }_{\partial \mathcal{B}}\overline{\mathit{t}}·\delta \mathit{v}\mathrm{d}S+{\int }_{\partial \mathcal{B}}\overline{\mathit{m}}:\delta {\mathit{R}}_e\mathrm{d}S}\\ \multicolumn{1}{c}{+{\int }_{\partial \mathcal{B}}\overline{\mathit{q}}:\delta {\mathbf{\kappa }}^e\mathrm{d}S.}\end{array}$$

(34)

Here $\mathit{\sigma }$ , $\mathbf{\chi }$ , $\mathbf{\Gamma }$ are the energetic conjugates of ${\mathit{E}}^e$ , ${\mathbf{\kappa }}^e$ , $\mathit{\theta }$ , respectively; $\mathit{b}$ is the body force, and $(\overline{\mathit{t}},\overline{\mathit{m}},\overline{\mathit{q}})$ are boundary (micro)tractions work-conjugate to $(\delta \mathit{v},\delta {\mathit{R}}_e,\delta {\mathbf{\kappa }}^e)$ .

Linearization of virtual measures

Use is made of the standard kinematic linearizations (details omitted here):

$$\delta {\mathit{E}}^e~\mathrm{sym}\mathrm{Grad}\delta \mathit{v},\delta {\mathbf{\kappa }}^e~\mathrm{Curl}\left(\delta {\mathit{R}}_e\right){\mathit{R}}_e^{\top },\delta \mathit{\theta }=\mathrm{Curl}\left(\delta {\mathbf{\kappa }}^e\right).$$

Integration by parts (adjoint identities)

Applying tensorial Green-type identities (adjoint of $\mathrm{Grad}$ and row-wise $\mathrm{Curl}$ ):

$${\int }_{\mathcal{B}}\mathit{\sigma }:\mathrm{Grad}\delta \mathit{v}\mathrm{d}V=-{\int }_{\mathcal{B}}\mathrm{Div}\mathit{\sigma }·\delta \mathit{v}\mathrm{d}V+{\int }_{\partial \mathcal{B}}\left(\mathit{\sigma }\mathit{n}\right)·\delta \mathit{v}\mathrm{d}S,$$

(35)

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathcal{B}}\mathbf{\chi }:\mathrm{Curl}\left(\delta {\mathit{R}}_e\right)\mathrm{d}V={\int }_{\mathcal{B}}\left({\mathrm{Curl}}^{*}\mathbf{\chi }\right):\delta {\mathit{R}}_e\mathrm{d}V}\\ \multicolumn{1}{c}{+{\int }_{\partial \mathcal{B}}\hat{\mathit{m}}(\mathbf{\chi },\mathit{n}):\delta {\mathit{R}}_e\mathrm{d}S,}\end{array}$$

(36)

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_{\mathcal{B}}\mathbf{\Gamma }:\mathrm{Curl}\left(\delta {\mathbf{\kappa }}^e\right)\mathrm{d}V={\int }_{\mathcal{B}}\left({\mathrm{Curl}}^{*}\mathbf{\Gamma }\right):\delta {\mathbf{\kappa }}^e\mathrm{d}V}\\ \multicolumn{1}{c}{+{\int }_{\partial \mathcal{B}}\hat{\mathit{q}}(\mathbf{\Gamma },\mathit{n}):\delta {\mathbf{\kappa }}^e\mathrm{d}S,}\end{array}$$

(37)

where ${\mathrm{Curl}}^{*}$ denotes the $L^2$ -adjoint of the row-wise $\mathrm{Curl}$ (the precise index form can be given in an Appendices A and B), $\mathit{n}$ is the outward unit normal, and $\hat{\mathit{m}},\hat{\mathit{q}}$ are boundary operators induced by the integration by parts.

Principle of virtual power

Enforcing ${\mathcal{P}}_{\mathrm{int}}={\mathcal{P}}_{\mathrm{ext}}$ for all admissible virtual fields and grouping volume and boundary terms yields:

Field (Euler–Lagrange) equations in $\mathcal{B}$ :

$$\mathrm{Div}\mathit{\sigma }+\mathit{b}=0,(\mathrm{macro}-\mathrm{balance})$$

(38)

$${\mathrm{Curl}}^{*}\mathbf{\chi }-\mathrm{Div}\mathbf{\Gamma }+\mathit{h}=0,(\mathrm{micro}-\mathrm{balance})$$

(39)

where $\mathit{h}$ collects possible micro-inertial or dissipative sources (zero in the purely energetic setting).

Natural boundary conditions on $\partial \mathcal{B}$ :

$$\mathit{\sigma }\mathit{n}=\overline{\mathit{t}},$$

(40)

$$\hat{\mathit{m}}(\mathbf{\chi },\mathit{n})=\overline{\mathit{m}},$$

(41)

$$\hat{\mathit{q}}(\mathbf{\Gamma },\mathit{n})=\overline{\mathit{q}}.$$

(42)

Comment on equation (54) and slip-system form

Equation (54) in the manuscript is the specialization of the micro-balance (equation (39)) written in the intermediate configuration ${\mathcal{B}}_i$ , with curvature measures projected on slip-system directions $({\mathit{m}}^i,{\mathit{s}}^i)$ and using the identity ${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathit{R}}_e\right){\mathit{R}}_e^{\top }$ . In that setting, the operator ${\mathrm{Curl}}^{*}$ acting on $\mathbf{\mu }$ produces the terms involving $\mathrm{Curl}({\mathit{m}}^i\times {\mathit{s}}^i)$ , thus justifying equation (32) as the micro-equilibrium of curvature-related (moment) stresses in ${\mathcal{B}}_i$ .

Having established the micro- and macro-balance equations through the principle of virtual power, the thermodynamic consistency of the framework is now addressed. In particular, the reduced dissipation inequality in the intermediate configuration is derived, ensuring that the constitutive choices introduced above are compatible with the second law of thermodynamics.

Reduced dissipation inequality in the intermediate configuration

Following standard nonequilibrium thermodynamics,^6,27 the local form of the second law in the intermediate configuration is obtained by combining the balance of energy, the entropy inequality, and the definition of the Helmholtz free energy. Specifically, the first law reads:

$$\tilde{\rho }\stackrel{·}{e}=\tilde{\mathit{\sigma }}:{\tilde{\mathbf{L}}}^{\top }-\mathrm{Div}\tilde{\mathbf{q}}+r,$$

(43)

while the entropy inequality is:

$$\tilde{\rho }\stackrel{·}{\eta }+\mathrm{Div}\left(\frac{\tilde{\mathbf{q}}}{\mathcal{T}}\right)\ge \frac{r}{\mathcal{T}}.$$

(44)

Introducing the Helmholtz free energy $\tilde{\psi }=e-\mathcal{T}\eta$ and substituting the heat supply $r$ from equation (43) into equation (44) leads to the reduced Clausius–Duhem inequality:

$$-\tilde{\rho }\stackrel{·}{\stackrel{~}{\psi }}-\tilde{\rho }\eta \stackrel{·}{\mathcal{T}}+\tilde{\mathit{\sigma }}:{\tilde{\mathbf{L}}}^{\top }-\frac{1}{\mathcal{T}}\tilde{\mathbf{q}}·\tilde{\nabla }\mathcal{T}\ge 0,$$

(45)

where $\tilde{\rho }$ is the mass density, $\tilde{\psi }$ the Helmholtz free energy, $\eta$ the entropy per unit mass, $\tilde{\mathit{\sigma }}$ the pulled-back Cauchy stress, $\tilde{\mathbf{L}}$ the velocity gradient, $\tilde{\mathbf{q}}$ the pulled-back heat flux, and $\mathcal{T}$ the absolute temperature.

The mapping of the Cauchy stress $\mathit{\sigma }$ to $\tilde{\mathit{\sigma }}$ is achieved via:

$$\tilde{\mathit{\sigma }}=\det \left({\mathbf{V}}^e\right){\left({\mathbf{V}}^e\right)}^{-1}\mathit{\sigma }{\left({\mathbf{V}}^e\right)}^{-1},$$

(46)

where ${\mathbf{V}}^e$ is the elastic stretch from ${\mathbf{F}}^e={\mathbf{V}}^e{\mathbf{R}}^e$ . ¹³ This ensures that $\tilde{\mathit{\sigma }}:{\tilde{\mathbf{L}}}^{\top }$ is work-conjugate to $\tilde{\mathbf{L}}$ , consistent with finite strain kinematics. ⁸ In the absence of couple stresses, $\tilde{\mathit{\sigma }}={\tilde{\mathit{\sigma }}}^{\top }$ , so only the symmetric part of $\tilde{\mathbf{L}}$ contributes to the stress power:

$$\tilde{\mathit{\sigma }}:{\tilde{\mathbf{L}}}^{\top }=\tilde{\mathit{\sigma }}:\mathrm{sym}\tilde{\mathbf{L}}=:\tilde{\mathit{\sigma }}:\tilde{\mathbf{D}}.$$

(47)

The stretching is split into elastic and dissipative parts, $\tilde{\mathbf{D}}={\tilde{\mathbf{D}}}^e+{\tilde{\mathbf{D}}}^p$ , and postulate the free energy:

$$\tilde{\psi }=\tilde{\psi }({\tilde{\mathbf{E}}}^e,{\epsilon }_{\mathrm{SS}},{\tilde{\mathbf{\kappa }}}_e,\mathrm{Curl}{\tilde{\mathbf{\kappa }}}_e,\mathcal{T}).$$

The chain rule yields:

$$\begin{array}{c}\tilde{\rho }\stackrel{·}{\stackrel{~}{\psi }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{E}}}^e}:{\stackrel{·}{\stackrel{~}{\mathbf{E}}}}^e+\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\epsilon }_{\mathrm{SS}}}{\stackrel{·}{\epsilon }}_{\mathrm{SS}}+\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{\kappa }}}_e}:{\stackrel{·}{\stackrel{~}{\mathbf{\kappa }}}}_e\hfill \\ \multicolumn{1}{c}{+\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial \mathrm{Curl}{\tilde{\mathbf{\kappa }}}_e}:\mathrm{Curl}{\stackrel{·}{\stackrel{~}{\mathbf{\kappa }}}}_e+\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial \mathcal{T}}\stackrel{·}{\mathcal{T}}.}\end{array}$$

(48)

Thermodynamic consistency requires the identification of reversible forces,

$$\begin{array}{c}\tilde{\mathit{\sigma }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{E}}}^e},\eta =-\frac{\partial \tilde{\psi }}{\partial \mathcal{T}},\tau =\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\epsilon }_{\mathrm{SS}}},\hfill \\ \multicolumn{1}{c}{\tilde{\mathbf{\chi }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial {\tilde{\mathbf{\kappa }}}_e},\tilde{\mathbf{\Gamma }}=\tilde{\rho }\frac{\partial \tilde{\psi }}{\partial \left(\mathrm{Curl}{\tilde{\mathbf{\kappa }}}_e\right)},}\end{array}$$

(49)

together with the conjugacy relation:

$$\tilde{\mathit{\sigma }}:{\tilde{\mathbf{D}}}^e=\tilde{\mathit{\sigma }}:{\stackrel{·}{\stackrel{~}{\mathbf{E}}}}^e.$$

(50)

Substituting equations (47) and (48) into equation (45), and applying equations (49) and (50), all reversible contributions cancel. The inequality reduces to the reduced dissipation inequality:

$$\mathcal{D}=\tilde{\mathit{\sigma }}:{\tilde{\mathbf{D}}}^p-\frac{1}{\mathcal{T}}\tilde{\mathbf{q}}·\tilde{\nabla }\mathcal{T}\ge 0$$

(51)

in the intermediate configuration.

Finally, the constitutive equation (49) couple naturally with the balance laws. The backstress $\tilde{\mathbf{\chi }}$ and moment stress $\tilde{\mathbf{\Gamma }}$ introduce internal length scales ${\ell }_1$ and ${\ell }_2$ , which regularize strain localization and enable size-dependent hardening. ¹⁸ These stresses also enter microforce balance equations derived from variational principles, ²⁶ ensuring equilibrium of internal defect fields. The framework thus guarantees thermodynamic admissibility while consistently bridging microscopic defect mechanics with macroscopic continuum response. ⁴

Crystal plasticity model

To formulate a thermodynamically consistent crystal plasticity model that accounts for size-dependent phenomena, the internal stresses derived in the preceding sections—namely, the backstress $\tilde{\mathbf{\chi }}$ and the higher-order moment stress $\tilde{\mathbf{\Gamma }}$ —are coupled to the evolution of crystallographic slip. The backstress $\tilde{\mathbf{\chi }}$ , proportional to the curvature tensor ${\tilde{\mathbf{\kappa }}}_e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T$ , ⁹ introduces a natural length scale associated with geometrically necessary dislocations (GNDs). Similarly, $\tilde{\mathbf{\Gamma }}$ , conjugate to the curvature-gradient measure $\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right)$ , embeds a second length scale linked to disclination effects. ¹⁷ This dual-length-scale approach distinguishes the proposed model from classical crystal plasticity frameworks, which lack intrinsic length scales and therefore cannot capture size-dependent hardening and the regularization of strain localization.^19,25

Recent advances in gradient-enhanced crystal plasticity have sought to incorporate length scales through various microstructural mechanisms. Acharya and Bassani ³ introduced the Curl of the elastic distortion into hardening laws to account for lattice incompatibility. Shu and Fleck ³⁷ proposed slip-gradient contributions to the critical shear stress, enhancing size effects in micro-scale deformation. Steinmann ³⁸ incorporated second-order gradients of the plastic deformation gradient to define internal stresses on slip systems, while Forest ²⁴ developed a micropolar crystal plasticity framework that accounts for micro-rotations. Acharya and Beaudoin ⁴ utilized lattice curvature to predict stage IV hardening, demonstrating the role of GNDs in large-strain deformation. Building on these works, the present model integrates curvature and curvature-gradient effects into a unified constitutive framework, leveraging the thermodynamic structure established earlier.

For a crystalline material with slip systems defined by unit normal ${\tilde{\mathit{m}}}^{\left(i\right)}$ and slip direction ${\tilde{\mathit{s}}}^{\left(i\right)}$ in the intermediate configuration, the plastic rate-of-deformation tensor ${\tilde{\mathit{D}}}_p$ and plastic spin tensor ${\tilde{\mathit{W}}}_p$ are expressed as the sum of contributions from all active slip systems:

$${\tilde{\mathit{D}}}_p=\sum _i{\stackrel{·}{\gamma }}^{\left(i\right)}\mathrm{sym}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})=\sum _i{\stackrel{·}{\gamma }}^{\left(i\right)}\frac{1}{2}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)}+{\tilde{\mathit{s}}}^{\left(i\right)}\otimes {\tilde{\mathit{m}}}^{\left(i\right)}),$$

(52)

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathit{W}}}_p=\sum _i{\stackrel{·}{\gamma }}^{\left(i\right)}\mathrm{skew}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})}\\ \multicolumn{1}{c}{=\sum _i{\stackrel{·}{\gamma }}^{\left(i\right)}\frac{1}{2}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)}-{\tilde{\mathit{s}}}^{\left(i\right)}\otimes {\tilde{\mathit{m}}}^{\left(i\right)}),}\end{array}$$

(53)

where ${\stackrel{·}{\gamma }}^{\left(i\right)}$ is the plastic slip rate on the $i$ th slip system, and $\mathrm{sym}(·)$ and $\mathrm{skew}(·)$ denote the symmetric and skew-symmetric parts of a tensor, respectively. These expressions, standard in crystal plasticity, ¹³ reflect the additive decomposition of the plastic velocity gradient ${\tilde{\mathit{L}}}_p={\tilde{\mathit{D}}}_p+{\tilde{\mathit{W}}}_p$ , where ${\tilde{\mathit{L}}}_p={\mathbf{R}}^e{\stackrel{·}{\mathbf{F}}}^p{\left({\mathbf{F}}^p\right)}^{-1}{\left({\mathbf{R}}^e\right)}^T$ is mapped to the intermediate configuration.

The slip rate ${\stackrel{·}{\gamma }}^{\left(i\right)}$ is governed by a power-law flow rule that incorporates both classical and gradient-dependent stresses to capture nonlocal effects:

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\gamma }}^{\left(i\right)}={\stackrel{·}{\gamma }}_0^{\left(i\right)}}\\ \multicolumn{1}{c}{{\left(\frac{|(\tilde{\mathit{\sigma }}-\tilde{\mathbf{\chi }}):({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})-\tilde{\tau }-\tilde{\mathbf{\Gamma }}:\mathrm{Curl}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})|}{{\sigma }_0}\right)}^{1/m},}\end{array}$$

(54)

where ${\stackrel{·}{\gamma }}_0^{\left(i\right)}$ is a reference slip rate, ${\sigma }_0$ is a reference stress, $m$ is the rate-sensitivity exponent, and $\mathrm{sgn}(·)$ ensures the correct direction of slip. The resolved shear stress on the $i$ th slip system is modified by three terms:

• $(\tilde{\mathit{\sigma }}-\tilde{\mathbf{\chi }}):({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})$ , representing the effective stress driving slip, where $\tilde{\mathit{\sigma }}$ is the Cauchy stress pulled back to the intermediate configuration, and $\tilde{\mathbf{\chi }}={\ell }_1^2\mu \left(\Theta \right){\tilde{\mathbf{\kappa }}}_e$ is the backstress due to GNDs. ⁹

• $\tilde{\tau }=\mu \left(\mathcal{T}\right)b\sqrt{{\rho }_{\mathrm{SS}}}$ , the isotropic hardening stress proportional to the density of statistically stored dislocations (SSDs) ${\rho }_{\mathrm{SS}}$ , with $b$ as the Burgers vector magnitude. ²⁶

• $\tilde{\mathbf{\Gamma }}:\mathrm{Curl}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})$ , a higher-order term where $\tilde{\mathbf{\Gamma }}={\ell }_2^4\mu \left(\mathcal{T}\right)\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}_e\right)$ is the moment stress conjugate to the second-order tensor $\mathrm{Curl}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})$ , capturing spatial variations in slip system orientation. ²⁴

The operator “:” denotes the double contraction. For two second-order tensors $\mathbf{A}$ and $\mathbf{B}$ , it is defined as:

$$\mathbf{A}:\mathbf{B}=A_{\mathrm{ij}}{B}_{\mathrm{ij}}.$$

The gradient term $\mathrm{Curl}({\tilde{\mathit{m}}}^{\left(i\right)}\otimes {\tilde{\mathit{s}}}^{\left(i\right)})$ quantifies the inhomogeneity of the slip system projector due to lattice curvature or GNDs, introducing a nonlocal energetic penalty that regularizes strain localization and enhances size-dependent hardening. ²⁷

To close the constitutive model, the evolution of the SSD density ${\rho }_{\mathrm{SS}}$ must be specified, as it determines $\tilde{\tau }$ . The accumulation of SSDs is governed by their mean free path $\lambda$ , which decreases with increasing dislocation density due to interactions between SSDs and GNDs. Following a parallel-process assumption, the effective mean free path is:

$$\lambda ={(\frac{1}{{\lambda }_{\mathrm{SS}}}+\frac{1}{{\lambda }_{\mathrm{gn}}})}^{-1},$$

(55)

where ${\lambda }_{\mathrm{SS}}$ and ${\lambda }_{\mathrm{gn}}$ are the mean free paths associated with SSDs and GNDs, respectively. These are assumed to scale inversely with the square root of their respective density:

$${\lambda }_{\mathrm{SS}}=\frac{1}{\sqrt{{\rho }_{\mathrm{SS}}}},{\lambda }_{\mathrm{gn}}=\frac{1}{\sqrt{{\rho }_{\mathrm{gn}}}},$$

(56)

where the GND density ${\rho }_{\mathrm{gn}}$ is related to the curvature tensor via:

$${\rho }_{\mathrm{gn}}=|({\tilde{\mathbf{\kappa }}}_e·{\tilde{\mathit{m}}}^{\left(i\right)})·{\tilde{\mathit{s}}}^{\left(i\right)}|,$$

(57)

approximating the scalar projection of ${\tilde{\mathbf{\kappa }}}_e$ onto the slip system, consistent with its role as a measure of lattice incompatibility.^9,35 The norm ensures a positive density, reflecting the magnitude of GND content.

The evolution of ${\rho }_{\mathrm{SS}}$ balances accumulation due to dislocation storage and annihilation due to dynamic recovery:

$${\stackrel{·}{\rho }}_{\mathrm{SS}}=\frac{c_1}{\lambda }\sum _i\left|{\stackrel{·}{\gamma }}^{\left(i\right)}\right|-c_2{\rho }_{\mathrm{SS}}\sum _i\left|{\stackrel{·}{\gamma }}^{\left(i\right)}\right|,$$

(58)

where $c_1$ and $c_2$ are material parameters controlling storage and recovery rates, respectively. Substituting $\lambda$ from equations (55) to (56) yields:

$$\frac{1}{\lambda }=\sqrt{{\rho }_{\mathrm{SS}}}+\sqrt{{\rho }_{\mathrm{gn}}},$$

(59)

so that:

$${\stackrel{·}{\rho }}_{\mathrm{SS}}=c_1(\sqrt{{\rho }_{\mathrm{SS}}}+\sqrt{{\rho }_{\mathrm{gn}}})\sum _i\left|{\stackrel{·}{\gamma }}^{\left(i\right)}\right|-c_2{\rho }_{\mathrm{SS}}\sum _i\left|{\stackrel{·}{\gamma }}^{\left(i\right)}\right|.$$

(60)

This evolution law introduces length-scale dependence through ${\rho }_{\mathrm{gn}}$ , which is tied to ${\tilde{\mathbf{\kappa }}}_e$ , enhancing hardening at large strains and capturing stage IV behavior.

The physical interpretation of this coupling mechanism between SSDs and GNDs is summarized below, clarifying how it differs from classical hardening models.

The interaction between SSDs and GNDs is central to the present framework. As GND density increases through lattice curvature, the mean free path of mobile dislocations decreases, which accelerates SSD storage. Consequently, regions with strong gradients accumulate both GNDs and SSDs, leading to enhanced isotropic hardening. This differs from classical models where SSDs evolve independently of GNDs. The proposed coupling mechanism captures stronger size effects and more pronounced hardening during large plastic deformation, consistent with experimental trends in small-scale crystals.

The backstress $\tilde{\mathbf{\chi }}$ depends on the curvature tensor ${\tilde{\mathbf{\kappa }}}_e$ , which requires the evolution of the elastic rotation ${\mathbf{R}}^e$ . From the kinematics, the elastic spin in the intermediate configuration is:

$${\tilde{\mathbf{\omega }}}_e={\stackrel{·}{\mathbf{R}}}^e{\left({\mathbf{R}}^e\right)}^T=\tilde{\mathit{W}}-{\tilde{\mathit{W}}}_e-{\tilde{\mathit{W}}}_p,$$

(61)

where $\tilde{\mathit{W}}$ is the total vorticity, ${\tilde{\mathit{W}}}_e={\left({\mathbf{V}}^e\right)}^{-1}{\stackrel{·}{\mathbf{V}}}^e$ is the elastic vorticity, and ${\tilde{\mathit{W}}}_p$ is given by equation (53). Thus:

$${\stackrel{·}{\mathbf{R}}}^e=(\tilde{\mathit{W}}-{\tilde{\mathit{W}}}_e-{\tilde{\mathit{W}}}_p){\mathbf{R}}^e,$$

(62)

allowing ${\mathbf{R}}^e$ to be updated. The curvature tensor is then computed as:

$${\tilde{\mathbf{\kappa }}}^e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T,$$

(63)

and the backstress follows from the constitutive relation:

$$\tilde{\mathbf{\chi }}={\ell }_1^2\mu \left(\mathcal{T}\right){\tilde{\mathbf{\kappa }}}_e,$$

(64)

where ${\ell }_1$ is an energetic length scale. ⁹ This completes the specification of the flow rule equation (54).

In summary, the proposed crystal plasticity model integrates gradient effects through the curvature tensor ${\tilde{\mathbf{\kappa }}}_e$ and its gradient, introducing two internal length scales that govern backstress and moment stress, respectively. The model captures size-dependent hardening and strain localization by coupling nonlocal stresses to the slip system evolution, aligning with the thermodynamic framework of Gurtin ²⁶ and the geometric interpretations of Forest. ²⁴ The curvature tensor, directly related to GNDs, ⁹ influences the elasto-plastic transition and enhances hardening at large strains, critical for predicting post-bifurcation responses and misoriented cell formation during severe plastic deformation. ¹⁸

Balance laws at multiple scales

The presence of geometrically necessary dislocations (GNDs) and disclinations in crystalline materials introduces internal stress and couple stress fields that classical continuum mechanics cannot capture.^24,26 These defects, arising from lattice incompatibilities, necessitate a multi-scale formulation of balance laws that incorporate contributions from the dislocation density tensor $\mathit{\alpha }$ and the disclination density tensor $\mathit{\theta }$ . This section develops a thermodynamically consistent framework, distinguishing between mesoscopic and macroscopic scales, to derive balance equations that embed length-scale-dependent stress measures. These measures ensure self-equilibration of internal fields, regularize strain localization, and enable the prediction of size-dependent phenomena such as enhanced hardening and subgrain formation.^9,18 The equations are formulated to be interconnected, self-contained, and aligned with the kinematic and thermodynamic frameworks established earlier, with references to prior work in gradient plasticity and defect mechanics providing theoretical support.

Mesoscopic scale

At the mesoscopic scale, balance laws are enhanced to account for incompatibilities in the elastic strain and curvature fields driven by GNDs and disclinations. The dislocation density tensor $\mathit{\alpha }$ , which quantifies the torsional incompatibility of the lattice, is defined in the intermediate configuration as:

$${\alpha }_{\mathrm{ij}}={ϵ}_{\mathrm{imn}}{\partial }_m{\epsilon }_{\mathrm{nj}}^e-{\kappa }_{\mathrm{ij}}^e+{\kappa }_{\mathrm{mm}}^e{\delta }_{\mathrm{ij}},$$

(65)

where ${\epsilon }_{\mathrm{ij}}^e$ is the elastic strain tensor, ${\kappa }_{\mathrm{ij}}^e=[\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T{]}_{\mathrm{ij}}$ is the curvature tensor derived from the elastic rotation tensor ${\mathbf{R}}^e\in \mathrm{SO}\left(3\right)$ , ${ϵ}_{\mathrm{imn}}$ is the Levi–Civita permutation symbol, and ${\delta }_{\mathrm{ij}}$ is the Kronecker delta. The term ${ϵ}_{\mathrm{imn}}{\partial }_m{\epsilon }_{\mathrm{nj}}^e$ captures the incompatibility of the elastic strain field, vanishing in a defect-free medium. ³⁵ The curvature terms $-{\kappa }_{\mathrm{ij}}^e+{\kappa }_{\mathrm{mm}}^e{\delta }_{\mathrm{ij}}$ adjust for rotational incompatibilities induced by lattice curvature, ensuring consistency with finite strain kinematics. ³ This definition generalizes the classical Kröner–Bilby construction to finite deformations, aligning with geometric interpretations of dislocations as torsional distortions in a Riemann–Cartan manifold. ³⁹ More details about the derivation are provided in the Appendices A and B in the back of this paper.

The internal stress ${\chi }_{\mathrm{ij}}$ , associated with the energetic cost of GNDs, is postulated to be linearly proportional to the dislocation density tensor:

$${\chi }_{\mathrm{ij}}=l_1^2\mu \left(\mathcal{T}\right){\alpha }_{\mathrm{ij}},$$

(66)

where $\mu$ is the shear modulus, and $l_1$ is a material parameter governing the strength of dislocation-induced stresses. ²⁷ This linear relation introduces a length scale associated with dislocation density, consistent with strain gradient plasticity frameworks. ¹⁹

The disclination density tensor $\mathit{\theta }$ , representing rotational incompatibilities due to non-uniform curvature, is defined as:

$${\theta }_{\mathrm{ij}}={ϵ}_{\mathrm{imn}}{\partial }_m{\kappa }_{\mathrm{nj}}^e,$$

(67)

where ${\partial }_m{\kappa }_{\mathrm{nj}}^e$ denotes the spatial gradient of the curvature tensor. This tensor quantifies disclinations, topological defects associated with curvature gradients in the lattice. ¹⁷ The corresponding couple stress ${\mathrm{\Gamma }}_{\mathrm{ij}}$ , work-conjugate to the disclination density, is given by:

$${\mathrm{\Gamma }}_{\mathrm{ij}}=l_2^4\mu \left(\mathcal{T}\right){\theta }_{\mathrm{ij}},$$

(68)

where $l_2$ is a material parameter introducing a second length scale for disclination effects. This formulation mirrors equation (66), ensuring a consistent treatment of defect-induced stresses. ²⁴

The mesoscopic balance laws are derived from the principle of virtual work, incorporating contributions from internal stresses and couple stresses. The linear momentum balance for the internal stress field is:

$${\partial }_j{\chi }_{\mathrm{ij}}=0,$$

(69)

ensuring that the dislocation-induced stress is divergence-free, thus self-equilibrated within the mesoscopic domain. ⁶ The angular momentum balance, accounting for the interaction between couple stresses and the skew-symmetric part of the internal stress, is expressed as:

$${\partial }_j{\mathrm{\Gamma }}_{\mathrm{ij}}+{ϵ}_{\mathrm{ijk}}{\chi }_{\mathrm{jk}}=0,$$

(70)

where ${ϵ}_{\mathrm{ijk}}{\chi }_{\mathrm{jk}}$ represents the moment contribution from the internal stress. This equation ensures that rotational incompatibilities due to disclinations are balanced by the couple stress field, consistent with micromorphic theories.^17,23

The interconnection between these equations is evident through the coupling of defect fields. Substituting ${\chi }_{\mathrm{ij}}=l_1^2\mu \left(\theta \right){\alpha }_{\mathrm{ij}}$ from equation (66) into equation (69) yields ${\partial }_j\left(l_1^2\mu \left(\theta \right){\alpha }_{\mathrm{ij}}\right)=0$ , which depends on the dislocation density ${\alpha }_{\mathrm{ij}}$ defined in equation (65). Similarly, substituting ${\mathrm{\Gamma }}_{\mathrm{ij}}=l_2^4\mu \left(\theta \right){\theta }_{\mathrm{ij}}$ from equation (68) and ${\theta }_{\mathrm{ij}}={ϵ}_{\mathrm{imn}}{\partial }_m{\kappa }_{\mathrm{nj}}^e$ from equation (67) into equation (70) yields ${\partial }_j\left(l_2^4\mu \left(\theta \right)b{ϵ}_{\mathrm{imn}}{\partial }_m{\kappa }_{\mathrm{nj}}^e\right)+{ϵ}_{\mathrm{ijk}}{\chi }_{\mathrm{jk}}=0$ , introducing higher-order gradients of the curvature tensor. The term ${ϵ}_{\mathrm{ijk}}{\chi }_{\mathrm{jk}}$ couples the dislocation and disclination fields via the curvature tensor ${\kappa }_{\mathrm{ij}}^e$ , which appears in both ${\alpha }_{\mathrm{ij}}$ and ${\theta }_{\mathrm{ij}}$ . This coupling ensures that the mesoscopic balance laws capture the interplay between translational and rotational defects, enhancing the model’s ability to describe complex defect interactions. ⁹

Macroscopic scale

At the macroscopic scale, the classical balance laws of continuum mechanics are recovered, as the effects of internal stresses and couple stresses are homogenized. The conservation of linear momentum is given by:

$${\partial }_j{\sigma }_{\mathrm{ij}}=0,$$

(71)

where ${\sigma }_{\mathrm{ij}}$ is the Cauchy stress tensor in the current configuration, governing the equilibrium of external forces. ¹³ The conservation of angular momentum enforces the symmetry of the Cauchy stress tensor: (In the present framework, objectivity of the internal power implies the symmetry of the macroscopic Cauchy stress, in agreement with Germain ³⁰ While non-symmetric macrostresses may be admitted in some generalized continuum theories such as Cosserat models, the proposed dual-length-scale formulation remains within the symmetric setting.)

$${\sigma }_{\mathrm{ij}}={\sigma }_{\mathrm{ji}},$$

(72)

ensuring that no net macroscopic couple stresses arise from external loading. ²⁷ These macroscopic balance laws involve not only the standard Cauchy stress but also implicit contributions from the higher-order fields $\mathbf{\chi }$ and $\mathbf{\Gamma }$ . While the latter are self-equilibrated in the sense that their divergence vanishes locally, they nevertheless enter the definition of the macrostress and thus contribute indirectly to macroscopic equilibrium.^27,30 The transition from mesoscopic to macroscopic scales is facilitated by the homogenization of defect-induced stresses, consistent with multi-scale plasticity models. ¹⁸

Remarks

The present formulation provides a consistent extension of higher-order crystal plasticity models, combining energetic and balance considerations in a way, that is, both mathematically structured and physically motivated. The internal stress $\mathbf{\chi }$ (equation (66)) and couple stress $\mathbf{\Gamma }$ (equation (68)) introduce two distinct length scales governed by ${\ell }_1$ and ${\ell }_2$ , respectively, which regularize the post-bifurcation response and enable predictions of size-dependent phenomena such as cell formation, lattice curvature, and enhanced hardening.^11,18 The mesoscopic balance laws (equations (69) and (70)) are derived consistently with the thermodynamic framework, ensuring compatibility with the free energy formulation. ⁶ The formulation is amenable to finite element implementation, as demonstrated in prior gradient plasticity models, and provides a robust foundation for studying localization, defect evolution, and gradient-enhanced plasticity in crystalline solids.

Analytical validation

To validate the proposed gradient crystal plasticity model with two internal length scales, an analytical solution for the simple shear of a single crystal strip is derived, this benchmark problem induces non-uniform deformation and engages the model’s gradient-enhanced features.^6,26 The strip, of height $H$ along the $x_2$ -axis ( $0\le x_2\le H$ ), is subjected to simple shear in the $x_1$ -direction, with a single slip system active and micro-hard boundary conditions that prevent plastic slip at the boundaries. This setup generates geometrically necessary dislocations (GNDs) and curvature gradients, allowing us to test the influence of the internal length scales ${\ell }_1$ and ${\ell }_2$ on size-dependent hardening and stress distributions. The solution is compared to classical crystal plasticity predictions and experimental observations of size effects in micro-scale shear.^11,25

Problem formulation

A single crystal strip, that is, infinite in the $x_1$ and $x_3$ directions, and of finite height $H$ in the $x_2$ -direction is considered. The deformation is assumed to be homogeneous along $x_1$ and $x_3$ , with displacement restricted to the $x_1$ -direction. Accordingly, the displacement field is taken as:

$$\mathbf{u}\left(x_2\right)=\left(\begin{array}{c}{u}_1\left(x_2\right)\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right),$$

(73)

from which the only non-zero component of the deformation gradient is the shear strain,

$$\gamma \left(x_2\right)={\partial }_2{u}_1\left(x_2\right).$$

(74)

A single active slip system is assumed in the intermediate configuration, characterized by the slip direction $\tilde{\mathbf{s}}=(1,0,0)$ and the slip-plane normal $\tilde{\mathbf{m}}=(0,1,0)$ . This choice corresponds to pure shear deformation aligned with crystallographic directions, adopted for analytical simplicity.

The boundary conditions are specified at the macroscopic and microscopic levels. At the macroscopic level, the displacement field satisfies:

$$u_1\left(0\right)=0,u_1\left(H\right)=\mathrm{\Gamma }H,$$

(75)

where $\mathrm{\Gamma }$ denotes the prescribed average shear strain. At the microscopic level, micro-hard conditions are imposed to prevent dislocation motion at the boundaries:

$$\stackrel{·}{\gamma }\left(0\right)=0,\stackrel{·}{\gamma }\left(H\right)=0.$$

(76)

The analysis is conducted under the assumptions of small elastic strain (i.e. ${\mathbf{V}}^e\approx \mathbf{I}$ ), isothermal conditions, and rate-independent plastic response. The kinematics follow the multiplicative split $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p$ , and the polar decomposition ${\mathbf{F}}^e={\mathbf{R}}^e{\mathbf{V}}^e$ is adopted, with the elastic rotation approximated for small angles by:

$${\mathbf{R}}^e\left(x_2\right)\approx \left(\begin{array}{ccc}1\hfill & -\omega \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{\omega \left(x_2\right)}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right).$$

(77)

(The approximation ${\mathbf{R}}^e\left(x_2\right)\approx \left(\begin{array}{ccc}1\hfill & -\omega \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{\omega \left(x_2\right)}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right)$ is based on the first-order Taylor expansion of the exact rotation matrix around the $x_3$ -axis:

$$\mathbf{R}\left(\omega \right)=\left(\begin{array}{ccc}\cos \omega \hfill & -\sin \omega \hfill & 0\hfill \\ \multicolumn{1}{c}{\sin \omega }& \cos \omega \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right).$$

For small rotation angles $\omega <<1$ , it holds that $\cos \omega \approx 1$ and $\sin \omega \approx \omega$ , so that $\mathbf{R}\left(\omega \right)\approx \left(\begin{array}{ccc}1\hfill & -\omega \hfill & 0\hfill \\ \multicolumn{1}{c}{\omega }& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right)$ This linearization is standard in small-strain crystal plasticity models where elastic distortions remain small, and it enables closed-form expressions for quantities such as lattice curvature and GND density.)

The lattice curvature is measured by the third-order tensor ${\tilde{\kappa }}^e$ , defined as:

$${\tilde{\kappa }}^e:=\mathrm{Curl}\left({\mathbf{R}}^e\right){\mathbf{R}}^{\mathrm{eT}}.$$

(78)

From this, the disclination density tensor is defined as:

$${\theta }_{\mathrm{ij}}:={\epsilon }_{\mathrm{imn}}{\partial }_m{\kappa }_{\mathrm{nj}}^e,$$

(79)

where ${\epsilon }_{\mathrm{imn}}$ is the Levi–Civita symbol, and summation is implied over repeated indices.

To account for size effects and internal resistance due to geometrically necessary dislocations (GNDs), the internal backstress tensor is introduced as:

$${\tilde{\chi }}_{\mathrm{ij}}:={\ell }_1^2\mu {\tilde{\kappa }}_{\mathrm{ij}}^e,$$

(80)

and the higher-order moment stress tensor is given by:

$${\tilde{\mathrm{\Gamma }}}_{\mathrm{ij}}:={\ell }_2^4\mu {\theta }_{\mathrm{ij}},$$

(81)

where ${\ell }_1$ and ${\ell }_2$ are internal length scales associated with the energetic cost of curvature and disclination density, respectively, and $\mu$ is the shear modulus.

The plastic flow rule is modified to incorporate these nonlocal effects, following the thermodynamically consistent frameworks proposed in Bammann ⁹ and Forest. ²⁴ In the absence of rate dependence, the resolved shear stress driving slip is given by:

$$\tau ={\tilde{\sigma }}_{12}-{\tilde{\chi }}_{12},$$

(82)

and is subject to a critical stress condition modified by the presence of statistically and geometrically stored dislocations, as developed in subsequent sections.

Analytical solution

The deformation is characterized by a shear displacement in the $x_1$ -direction, varying across the height $x_2$ . The corresponding deformation gradient is given by:

$$\mathbf{F}=\left(\begin{array}{ccc}1\hfill & \gamma \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right),\mathrm{with}\gamma \left(x_2\right)={\partial }_2{u}_1.$$

(83)

Applying the standard multiplicative decomposition $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^p$ , a (quantity/definition) is introduced:

$${\mathbf{F}}^p=\left(\begin{array}{ccc}1\hfill & {\gamma }^p\left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right),{\mathbf{F}}^e=\left(\begin{array}{ccc}1\hfill & {\gamma }^e\left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right),$$

(84)

where $\gamma ={\gamma }^e+{\gamma }^p$ . For a perfectly aligned slip system, one can take ${\gamma }^p=\gamma$ , with elastic strain vanishing in the shear direction.

The elastic rotation tensor, linearized for small angles, is approximated by:

$${\mathbf{R}}^e\approx \left(\begin{array}{ccc}1\hfill & -\omega \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{\omega \left(x_2\right)}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right).$$

(85)

From this, the curvature tensor is defined as:

$${\tilde{\kappa }}^e=\mathrm{Curl}\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^{\mathrm{T}}=\left(\begin{array}{ccc}0\hfill & 0\hfill & -\omega \prime \left(x_2\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \omega \prime \left(x_2\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right),$$

(86)

and the disclination density follows as:

$$\theta =\mathrm{Curl}\left({\tilde{\kappa }}^e\right)=\left(\begin{array}{ccc}0\hfill & -\omega \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \omega \left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right).$$

(87)

The backstress and couple stress tensors are thus given by:

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\chi }={\ell }_1^2\mu {\tilde{\kappa }}^e=\left(\begin{array}{ccc}0\hfill & 0\hfill & -{\ell }_1^2\mu \omega \prime \left(x_2\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\ell }_1^2\mu \omega \prime \left(x_2\right)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right),}\\ \multicolumn{1}{c}{\tilde{\mathrm{\Gamma }}={\ell }_2^4\mu \theta =\left(\begin{array}{ccc}0\hfill & -{\ell }_2^4\mu \omega ″\left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\ell }_2^4\mu \omega ″\left(x_2\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right).}\end{array}$$

(88)

The local flow rule, accounting for backstress and SSD hardening, reads:

$${\tilde{\sigma }}_{12}-{\ell }_1^2\mu \omega \prime \left(x_2\right)-c_{\tau }\mu b\sqrt{{\rho }_{\mathrm{SS}}}={\sigma }_0.$$

(89)

The mesoscopic balance laws yield the following nontrivial conditions:

$${\partial }_2{\chi }_{13}=-{\ell }_1^2\mu \omega ″\left(x_2\right)=0,$$

(90)

$${\partial }_2{\mathrm{\Gamma }}_{12}=-{\ell }_2^4\mu \omega ″\prime \left(x_2\right)=0.$$

(91)

These imply:

$$\omega \left(x_2\right)={\mathrm{Ax}}_2+B,\mathrm{and}\mathrm{hence}{\partial }_2\omega =A=\mathrm{const}.$$

The macroscopic equilibrium ${\partial }_2{\tilde{\sigma }}_{12}=0$ leads to a constant shear stress ${\tilde{\sigma }}_{12}=S$ , so equation (89) simplifies to:

$$S={\sigma }_0+{\ell }_1^2\mu A+c_{\tau }\mu b\sqrt{{\rho }_{\mathrm{SS}}}.$$

(92)

At steady state, the statistically stored dislocation density satisfies:

$$c_1(\sqrt{{\rho }_{\mathrm{SS}}}+\sqrt{\left|A\right|})=c_2{\rho }_{\mathrm{SS}}.$$

(93)

Solving the quadratic yields:

$$\sqrt{{\rho }_{\mathrm{SS}}}=\frac{c_1}{2c_2}(1+\sqrt{1+\frac{4c_2}{c_1}\sqrt{\left|A\right|}}),$$

(94)

and substituting back into the flow rule gives the closed-form expression:

$$S={\sigma }_0+{\ell }_1^2\mu A+c_{\tau }\mu b\frac{c_1}{2c_2}(1+\sqrt{1+\frac{4c_2}{c_1}\sqrt{\left|A\right|}}).$$

(95)

Approximating ${\partial }_2\omega \approx {\kappa }_0/H$ , where ${\kappa }_0$ is a characteristic curvature yields: (The quantity ∂₂ω is approximated as ∂₂ω ≈ κ₀/H, where κ₀ denotes the accumulated curvature across the strip height H, assuming a linear variation of ω(x₂).)

$$S={\sigma }_0+{\ell }_2^1\mu \frac{{\kappa }_0}{H}+c_{\tau }\mu b\frac{c_1}{2c_2}(1+\sqrt{1+\frac{4c_2}{c_1}\sqrt{\frac{{\kappa }_0}{H}}})$$

(96)

Finally, the average shear strain is composed of the elastic and plastic contributions:

$$\mathrm{\Gamma }=\frac{S}{\mu }+\frac{1}{H}{\int }_0^H{\gamma }^p\left(x_2\right){\mathrm{dx}}_2,$$

(97)

where the plastic strain field is taken to be:

$${\gamma }^p\left(x_2\right)={\mathrm{Cx}}_2(H-x_2),$$

satisfying the boundary conditions ${\gamma }^p\left(0\right)={\gamma }^p\left(H\right)=0$ .

Validation and discussion

The closed-form expression for the shear stress,

$$S={\sigma }_0+{\ell }_1^2\mu \frac{{\kappa }_0}{H}+c_{\tau }\mu b\left(\frac{c_1}{2c_2}(1+\sqrt{1+\frac{4c_2}{c_1}\sqrt{\frac{{\kappa }_0}{H}}})\right),$$

(98)

exhibits a clear dependence on the sample height $H$ , reflecting length-scale effects intrinsic to strain gradient plasticity. As $H\to 0$ , both the backstress term ${\ell }_1^2\mu {\kappa }_0/H$ and the gradient-induced contribution from statistically stored dislocations (SSDs) increase, leading to a pronounced size effect. This trend is consistent with experimental observations in micro-shear tests^20,25 and mirrors the phenomenology associated with the Hall–Petch effect.

In contrast, classical crystal plasticity models, which neglect gradient contributions (i.e. ${\ell }_1={\ell }_2=0$ ), predict the size-independent relation:

$$S={\sigma }_0+c_{\tau }\mu b\sqrt{{\rho }_{\mathrm{SS}}},$$

(99)

failing to capture the observed strengthening in micron-scale specimens. The inclusion of the backstress ${\ell }_1\mu {\partial }_2\omega$ , which arises from geometrically necessary dislocations (GNDs), accounts for curvature-induced internal stresses near constrained boundaries. This effect aligns with theoretical predictions from Acharya and Bassani’s field dislocation mechanics framework. ³

Moreover, although the disclination-related moment stress ${\ell }_2^4\mu {\partial }_2^2{\kappa }^e$ appears in the governing balance laws, its contribution vanishes in the present simple shear setting since ${\partial }_2^3\omega =0$ . This is consistent with the fact that shear involves only the curvature tensor ${\mathbf{\kappa }}^e$ , while the disclination density $\mathit{\theta }$ remains inactive. By contrast, in problems with non-uniform curvature distributions—such as torsional deformation with micro-clamping, bending of a crystal plate, or bicrystal tilt boundaries—higher-order incompatibilities arise, making $\mathit{\theta }$ nonzero and localized. ²⁴ These cases highlight the distinct role of the second internal length scale ${\ell }_2$ in controlling boundary layers and diffuse interfaces, and explicit analytical reductions of the present model for these problems are provided in Section 8.

The model’s predictions are in qualitative agreement with both numerical solutions of gradient-enhanced plasticity theories ⁶ and experimental studies reporting grain- and sample-size effects. This consistency confirms the model’s capability to embed intrinsic material length scales into a physically motivated continuum framework, thereby extending the descriptive power of classical plasticity into regimes where gradient effects cannot be neglected.

Additional analytical illustrations of curvature and disclination effects

To demonstrate in concrete terms how the proposed dual-length-scale model activates both curvature ( ${\kappa }^e$ ) and disclination density ( $\theta$ ), simplified boundary-value problems admitting closed-form solutions are presented. These examples are directly based on the free energy density:

$$\psi =\frac{1}{2}{\mathbf{E}}^e:\mathbb{C}:{\mathbf{E}}^e+\frac{1}{2}{\ell }_1^2{\kappa }^e:\mathbb{H}:{\kappa }^e+\frac{1}{2}{\ell }_2^4\theta :\mathbb{K}:\theta ,$$

(100)

with generalized stresses:

$$\mathit{\sigma }=\frac{\partial \psi }{\partial {\mathbf{E}}^e},\mathbf{\chi }=\frac{\partial \psi }{\partial {\mathbf{\kappa }}^e},\mathbf{\Gamma }=\frac{\partial \psi }{\partial \mathit{\theta }}.$$

The associated balance laws (derived in Section 6 of the manuscript) reduce to:

$$\mathrm{Div}\mathit{\sigma }+\mathit{b}=0,{\mathrm{Curl}}^{*}\mathbf{\chi }-\mathrm{Div}\mathbf{\Gamma }=0,$$

(101)

for macro- and micro-equilibrium, respectively.

Remark on one-dimensional reduction

In the full three-dimensional theory, the curvature and disclination tensors are defined as:

$${\mathbf{\kappa }}^e=\mathrm{Curl}\left({\mathit{R}}_e\right){\mathit{R}}_e^{\top },\mathit{\theta }=\mathrm{Curl}\left({\mathbf{\kappa }}^e\right),$$

with ${\mathit{R}}_e\in \mathrm{SO}\left(3\right)$ the elastic rotation. For analytical tractability, attention is restricted to high-symmetry kinematics (torsion, bending, and bicrystal tilt boundaries) in which the lattice rotates about a fixed axis. In such cases, ${\mathit{R}}_e$ can be parametrized by a single rotation angle $\alpha \left(x\right)$ depending on one spatial coordinate, so that:

$${\mathbf{\kappa }}^e~\alpha \prime \left(x\right),\mathit{\theta }~\alpha ″\left(x\right).$$

Here $\alpha \prime \left(x\right)$ plays the role of a scalar curvature (linked to GND density) and $\alpha ″\left(x\right)$ represents the scalar disclination density.

The reduced free energy then simplifies to:

$$\psi =\frac{1}{2}\mu {\ell }_1^2(\alpha \prime {)}^2+\frac{1}{2}\mu {\ell }_2^4(\alpha ″{)}^2,$$

(102)

which preserves the structure of the governing equations while allowing closed-form analytical solutions.

The governing Euler–Lagrange equation associated with the reduced one-dimensional energy functional:

$$\psi =\frac{1}{2}\mu {\ell }_1^2{(\alpha \prime )}^2+\frac{1}{2}\mu {\ell }_2^4{(\alpha ″)}^2$$

follows from stationarity with respect to variations in the rotation field $\alpha \left(x\right)$ . Integration by parts yields:

$$-\mu {\ell }_1^2\alpha ″+\mu {\ell }_2^4{\alpha }^{\left(4\right)}=0,$$

(103)

which can be recast in normalized form as:

$${\alpha }^{\left(4\right)}-k^2\alpha ″=0,k^2=\frac{{\ell }_1^2}{{\ell }_2^4}.$$

(104)

The characteristic polynomial of equation (104) has roots $r=0,0,\pm k$ , so that the general solution is a quadratic polynomial superposed with exponential modes $e^{\pm \mathrm{kx}}$ . These exponential contributions imply the presence of boundary layers whenever nontrivial boundary conditions are imposed.

It is convenient to introduce the emergent intrinsic length scale:

$$L_c=\frac{1}{k}=\frac{{\ell }_2^2}{{\ell }_1}.$$

(105)

Whereas ${\ell }_1$ and ${\ell }_2$ enter the free energy as independent penalties on curvature and curvature gradients, respectively, it is their ratio that controls the exponential decay of lattice curvature from boundaries. In other words, $L_c$ governs the thickness of boundary layers over which the disclination density $\theta ~\alpha ″$ concentrates.

From a physical perspective, $L_c$ can be interpreted as the effective penetration depth of disclination-induced incompatibilities into the bulk. Regions of strong orientation gradients (e.g. near free surfaces, grain boundaries, or interfaces) therefore develop localized boundary layers of width $O\left(L_c\right)$ , within which both curvature and disclination fields are non-negligible. This emergent scale provides the natural bridge between the microscopic parameters $({\ell }_1,{\ell }_2)$ and experimentally observable features such as the width of orientation-gradient boundary layers revealed by EBSD or Laue diffraction.

Torsion of a cylindrical bar

Kinematics

A cylindrical bar of length $L$ and radius $R$ subject to an imposed twist about its axis is considered. A simple kinematic approximation for the displacement field is:

$$u_{\phi }(r,z)=\gamma \left(z\right)r,u_r=0,u_z=0,$$

where $r$ and $z$ are the radial and axial coordinates, and $\gamma \left(z\right)$ is the twist per unit length, which may vary along $z$ . This deformation corresponds to an elastic rotation $\alpha \left(z\right)$ about the bar axis, with $\alpha \prime \left(z\right)=\gamma \left(z\right)$ . In the one-dimensional reduction,

$${\mathbf{\kappa }}^e~\alpha \prime \left(z\right),\mathit{\theta }~\alpha ″\left(z\right).$$

The reduced free energy density is:

$$\psi \left(\alpha \right)=\frac{1}{2}\mu {\ell }_1^2(\alpha \prime {)}^2+\frac{1}{2}\mu {\ell }_2^4(\alpha ″{)}^2,$$

and the Euler–Lagrange equation becomes:

$${\alpha }^{\left(4\right)}-k^2\alpha ″=0,k^2=\frac{{\ell }_1^2}{{\ell }_2^4},L_c=\frac{{\ell }_2^2}{{\ell }_1}.$$

Boundary conditions

Two representative cases are considered:

• Uniform torsion (free micro-tractions): prescribed end rotations $\alpha \left(0\right)=0$ , $\alpha \left(L\right)={\alpha }_0$ , with $\alpha ″\left(0\right)=\alpha ″\left(L\right)=0$ .

• Nonuniform torsion (micro-clamping): prescribed end rotations $\alpha \left(0\right)=0$ , $\alpha \left(L\right)={\alpha }_0$ , together with slope constraints $\alpha \prime \left(0\right)=\alpha \prime \left(L\right)=0$ .

Solution

For uniform torsion, the solution of the governing ODE is linear:

$$\alpha \left(z\right)=\frac{{\alpha }_0}{L}z,{\mathbf{\kappa }}^e=\alpha \prime =\frac{{\alpha }_0}{L},\mathit{\theta }=\alpha ″=0.$$

In this case, the bar develops a constant GND density but no disclinations.

For nonuniform torsion with micro-clamping, the higher-order ODE admits solutions containing hyperbolic functions. Imposing the boundary conditions yields:

$$\alpha ″\left(z\right)=\mathit{\theta }\left(z\right)=-{\alpha }_0{k}^2\frac{\cosh \left(\mathrm{kz}\right)-\cosh \left(k(L-z)\right)}{D},$$

$$\alpha \prime \left(z\right)={\mathbf{\kappa }}^e\left(z\right)=-{\alpha }_{0}k\frac{-\sinh \left(\mathrm{kL}\right)+\sinh \left(\mathrm{kz}\right)+\sinh \left(k(L-z)\right)}{D},$$

$$\begin{array}{c}\multicolumn{1}{c}{\alpha \left(z\right)={\alpha }_0}\\ \multicolumn{1}{c}{\frac{\mathrm{kz}\sinh \left(\mathrm{kL}\right)-\cosh \left(\mathrm{kL}\right)-\cosh \left(\mathrm{kz}\right)+\cosh \left(k(L-z)\right)+1}{D},}\end{array}$$

with denominator:

$$D=Lk\sinh \left(\mathrm{kL}\right)-2\cosh \left(\mathrm{kL}\right)+2.$$

In this case, $\mathit{\theta }\left(z\right)$ is nonzero near the bar ends and decays exponentially into the interior over a thickness $L_c={\ell }_2^2/{\ell }_1$ .

Interpretation

These two solutions illustrate the contrasting roles of the internal length scales. With free micro-tractions, torsion is uniform: GNDs are present but disclinations vanish. When micro-clamping is imposed, the slope is blocked at the ends, enforcing curvature gradients that activate disclination density $\mathit{\theta }$ . The resulting boundary layers have thickness $L_c$ , and their peak magnitude scales as ${\alpha }_0/\left({\mathrm{LL}}_c\right)$ . Thus, the second internal length scale ${\ell }_2$ directly controls the size-dependent hardening under torsional loading with higher-order boundary constraints.

Bending of a crystal plate

Kinematics

Consider a crystalline plate of thickness $h$ with the transverse coordinate $y\in [-h/2,h/2]$ . The plate is subjected to pure bending about the $z$ -axis, so that the lattice undergoes a tilt in the $x$ – $y$ plane. This rotation can be described by a scalar angle $\alpha \left(y\right)$ that varies across the thickness. Differentiation of $\alpha \left(y\right)$ provides direct access to the curvature and disclination density in this one-dimensional reduction.

Reduced problem

In this simplified setting, the curvature and disclination tensors reduce to:

$${\mathbf{\kappa }}^e\left(y\right)~\alpha \prime \left(y\right),\mathit{\theta }\left(y\right)~\alpha ″\left(y\right).$$

The reduced free energy density therefore reads:

$$\psi \left(y\right)=\frac{1}{2}\mu {\ell }_1^2{(\alpha \prime \left(y\right))}^2+\frac{1}{2}\mu {\ell }_2^4{(\alpha ″\left(y\right))}^2,$$

which penalizes both curvature (through ${\ell }_1$ ) and curvature gradients (disclinations, through ${\ell }_2$ ). The associated Euler–Lagrange equation is:

$$-\mu {\ell }_1^2\alpha ″\left(y\right)+\mu {\ell }_2^4{\alpha }^{\left(4\right)}\left(y\right)=0\iff {\alpha }^{\left(4\right)}\left(y\right)-k^2\alpha ″\left(y\right)=0,$$

with $k^2={\ell }_1^2/{\ell }_2^4$ . This fourth-order ODE reflects the fact that disclination effects ( $\alpha ″$ ) enter as higher gradients of the lattice rotation.

Boundary conditions

To model macroscopic bending, the following conditions are imposed:

• Macroscopic curvature: $\alpha \prime (\pm h/2)=\pm {\kappa }_0$ , ensuring that the average slope of the rotation field matches the imposed bending curvature ${\kappa }_0$ .

• Free micro-tractions: $\alpha ″(\pm h/2)=0$ , which express the absence of higher-order (disclination-related) tractions at the free surfaces.

Together, these four conditions provide the closure required to solve the fourth-order ODE.

Solution

The general solution of ${\alpha }^{\left(4\right)}-k^2\alpha ″=0$ can be written as:

$$\alpha \left(y\right)=A+\mathrm{By}+C\cosh \left(\mathrm{ky}\right)+D\sinh \left(\mathrm{ky}\right).$$

Enforcing symmetry about the mid-plane ( $\alpha (-y)=-\alpha \left(y\right)$ ) yields $A=0$ and $D=0$ , so that:

$$\alpha \left(y\right)=\mathrm{By}+C\cosh \left(\mathrm{ky}\right).$$

Differentiation gives:

$$\alpha \prime \left(y\right)=B+\mathrm{Ck}\sinh \left(\mathrm{ky}\right),\alpha ″\left(y\right)={\mathrm{Ck}}^2\cosh \left(\mathrm{ky}\right).$$

Applying the boundary conditions $\alpha \prime (\pm h/2)=\pm {\kappa }_0$ and $\alpha ″(\pm h/2)=0$ determines $B$ and $C$ , leading to the closed-form solution:

$$\alpha \prime \left(y\right)={\kappa }_0\frac{\cosh \left(\mathrm{ky}\right)}{\cosh (\mathrm{kh}/2)},\theta \left(y\right)=\alpha ″\left(y\right)={\kappa }_{0}k\frac{\sinh \left(\mathrm{ky}\right)}{\cosh (\mathrm{kh}/2)}.$$

Interpretation

The field $\alpha \prime \left(y\right)$ represents the through-thickness curvature distribution, while $\theta \left(y\right)$ quantifies the associated disclination density. In the bulk ( $\left|y\right|<<h/2$ ), the curvature is nearly uniform, recovering the classical bending prediction. Near the surfaces, however, $\alpha \prime \left(y\right)$ adjusts within boundary layers of thickness:

$$L_c=\frac{{\ell }_2^2}{{\ell }_1},$$

where $\theta \left(y\right)$ peaks. Thus, the second length scale ${\ell }_2$ directly controls the surface concentration of disclinations, producing enhanced size-dependent hardening in thin plates. As $h/L_c\to \mathrm{\infty }$ , the model converges to the classical prediction, while for $h/L_c=O\left(1\right)$ , strong size effects are observed.

Bicrystal tilt boundary

Kinematics

Consider two crystals meeting at $y=0$ with a misorientation $\mathrm{\Delta }\alpha$ about the $z$ -axis. Instead of an abrupt jump in orientation, the lattice rotation is allowed to vary smoothly across the boundary region, described by $\alpha \left(y\right)$ with far-field values:

$$\alpha \left(y\right)\to \pm \frac{\mathrm{\Delta }\alpha }{2}\mathrm{as}y\to \pm \mathrm{\infty }.$$

Reduced problem

As in the previous cases,

$${\mathbf{\kappa }}^e\left(y\right)~\alpha \prime \left(y\right),\mathit{\theta }\left(y\right)~\alpha ″\left(y\right),$$

and the free energy density is:

$$\psi \left(y\right)=\frac{1}{2}\mu {\ell }_1^2{(\alpha \prime \left(y\right))}^2+\frac{1}{2}\mu {\ell }_2^4{(\alpha ″\left(y\right))}^2.$$

Governing equation

The Euler–Lagrange equation is again:

$${\alpha }^{\left(4\right)}\left(y\right)-k^2\alpha ″\left(y\right)=0,k^2={\ell }_1^2/{\ell }_2^4.$$

The natural boundary conditions at $y\to \pm \mathrm{\infty }$ enforce the prescribed far-field orientations.

Solution

The minimization problem with these asymptotic boundary conditions yields a diffuse interface of characteristic thickness $\lambda$ :

$$\alpha \left(y\right)=\frac{\mathrm{\Delta }\alpha }{2}\tanh \left(\frac{y}{\lambda }\right),$$

with:

$$\begin{array}{c}\multicolumn{1}{c}{{\kappa }^e\left(y\right)=\alpha \prime \left(y\right)=\frac{\mathrm{\Delta }\alpha }{2\lambda }{\mathrm{sech}}^2\left(\frac{y}{\lambda }\right),\theta \left(y\right)=\alpha ″\left(y\right)}\\ \multicolumn{1}{c}{=-\frac{\mathrm{\Delta }\alpha }{{\lambda }^2}{\mathrm{sech}}^2\left(\frac{y}{\lambda }\right)\tanh \left(\frac{y}{\lambda }\right).}\end{array}$$

The interface thickness is controlled by:

$$\lambda =\frac{{\ell }_2^2}{{\ell }_1}.$$

Interpretation

The bicrystal misorientation is accommodated by a localized GND distribution ${\kappa }^e\left(y\right)$ centered at the boundary, accompanied by an antisymmetric disclination density $\theta \left(y\right)$ . The model predicts that the diffuse boundary width is not arbitrary but determined by the intrinsic length ratio ${\ell }_2^2/{\ell }_1$ . This provides a direct physical interpretation: EBSD or Laue diffraction measurements of subgrain boundary widths can be used to calibrate ${\ell }_1$ and ${\ell }_2$ .

Summary

These analytical examples show that: (i) simple shear or uniform torsion activate only ${\kappa }^e$ (GNDs), (ii) torsion with axial gradients, bending, and grain boundaries activate $\theta$ (disclinations), and (iii) the internal lengths ${\ell }_1,{\ell }_2$ set the magnitude of GND-related hardening and the thickness of disclination-rich boundary layers. These connections provide direct pathways to calibrate ${\ell }_1,{\ell }_2$ from EBSD or Laue diffraction measurements of misorientation gradients and boundary widths.

Numerical implementation of the multi-length scale model

Constitutive model summary

This subsection summarizes the constitutive equations and the local update procedure implemented at each Gauss point in the finite-element discretization. The presentation is structured to reflect the sequence of operations required in a computational implementation.

State variables

At each integration point, the constitutive state is characterized by the set:

$$\mathcal{Z}=\{{\tilde{\mathbf{E}}}^e,{\mathbf{R}}^e,{\tilde{\mathbf{\kappa }}}^e,\mathit{\theta },\gamma ,{\rho }_{\mathrm{SS}}\},$$

(106)

where ${\tilde{\mathbf{E}}}^e$ is the elastic strain measure, ${\mathbf{R}}^e$ the elastic rotation, ${\tilde{\mathbf{\kappa }}}^e$ the lattice curvature tensor, $\mathit{\theta }=\mathrm{Curl}\left({\tilde{\mathbf{\kappa }}}^e\right)$ the curvature-gradient (disclination density), $\gamma$ the accumulated slip (single slip), and ${\rho }_{\mathrm{SS}}$ the statistically stored dislocation (SSD) density.

Free energy and thermodynamic forces

The Helmholtz free-energy density is taken in quadratic form as:

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\rho }\tilde{\psi }=\mu \parallel {\tilde{\mathbf{E}}}^e{\parallel }^2+\frac{K}{2}{\left({\tilde{\mathbf{E}}}^e\right)}^2+\frac{1}{2}{c}_{\tau }\mu {\epsilon }_{\mathrm{SS}}^2+\frac{1}{2}\mu {\ell }_1^2\parallel {\tilde{\mathbf{\kappa }}}^e{\parallel }^2}\\ \multicolumn{1}{c}{+\frac{1}{2}\mu {\ell }_2^4\parallel \mathit{\theta }{\parallel }^2,}\end{array}$$

(107)

with ${\epsilon }_{\mathrm{SS}}=b\sqrt{{\rho }_{\mathrm{SS}}}$ . Differentiation of equation (107) yields the constitutive forces:

$$\tilde{\mathit{\sigma }}=2\mu \left({\tilde{\mathbf{E}}}^e\right)+K\left({\tilde{\mathbf{E}}}^e\right)\mathbf{I},$$

(108)

$$\tilde{\mathbf{\chi }}=\mu {\ell }_1^2{\tilde{\mathbf{\kappa }}}^e,$$

(109)

$$\tilde{\mathbf{\Gamma }}=\mu {\ell }_2^4\mathit{\theta },$$

(110)

$$\tilde{\tau }=c_{\tau }\mu {\epsilon }_{\mathrm{SS}}.$$

(111)

Here $\tilde{\mathit{\sigma }}$ is the Cauchy stress, $\tilde{\mathbf{\chi }}$ the curvature-induced backstress, $\tilde{\mathbf{\Gamma }}$ the higher-order moment stress associated with curvature gradients, and $\tilde{\tau }$ the isotropic hardening stress.

Resolved shear stress and flow rule

For a single active slip system with slip direction $\mathbf{s}$ and slip-plane normal $\mathbf{m}$ , the effective resolved shear stress is defined as:

$${\tau }_{\mathrm{eff}}=(\tilde{\mathit{\sigma }}-\tilde{\mathbf{\chi }}):(\mathbf{m}\otimes \mathbf{s})-\tilde{\tau }-\tilde{\mathbf{\Gamma }}:(\mathbf{m}\otimes \mathbf{s}).$$

(112)

Plastic flow is governed by a rate-dependent power-law relation,

$$\stackrel{·}{\gamma }={\stackrel{·}{\gamma }}_0{\left(\frac{\left|{\tau }_{\mathrm{eff}}\right|}{{\sigma }_0}\right)}^{1/m}\left({\tau }_{\mathrm{eff}}\right),$$

(113)

where ${\stackrel{·}{\gamma }}_0$ is a reference slip rate, ${\sigma }_0$ a reference stress, and $m$ the rate-sensitivity exponent. The use of a viscoplastic regularization ensures robust numerical convergence of the local update.

SSD evolution

The evolution of the SSD density is driven by slip activity and accounts for interactions with geometrically necessary dislocations (GNDs). The GND density is estimated from the curvature tensor as:

$${\rho }_{\mathrm{gn}}=|({\tilde{\mathbf{\kappa }}}^e·\mathbf{m})·\mathbf{s}|.$$

(114)

The SSD density then evolves according to:

$${\stackrel{·}{\rho }}_{\mathrm{SS}}=c_1(\sqrt{{\rho }_{\mathrm{SS}}}+\sqrt{{\rho }_{\mathrm{gn}}})\left|\stackrel{·}{\gamma }\right|-c_2{\rho }_{\mathrm{SS}}\left|\stackrel{·}{\gamma }\right|,$$

(115)

where $c_1$ and $c_2$ are material parameters controlling dislocation storage and dynamic recovery, respectively.

Update of elastic rotation and curvature

The elastic rotation evolves from the elastic spin relation:

$${\stackrel{·}{\mathbf{R}}}^e={\mathbf{\omega }}^e{\mathbf{R}}^e,{\mathbf{\omega }}^e=\mathbf{W}-{\mathbf{W}}^e-{\mathbf{W}}^p,$$

(116)

where $\mathbf{W}$ is the total spin and ${\mathbf{W}}^p$ the plastic spin associated with slip. The curvature tensor is then obtained as:

$${\tilde{\mathbf{\kappa }}}^e=\left({\mathbf{R}}^e\right){\left({\mathbf{R}}^e\right)}^T,\mathit{\theta }=\left({\tilde{\mathbf{\kappa }}}^e\right),$$

(117)

providing the quantities required to evaluate $\tilde{\mathbf{\chi }}$ and $\tilde{\mathbf{\Gamma }}$ through equations (109) and (110).

Algorithmic structure (Gauss point update)

For a given time increment $\mathrm{\Delta }t$ and Gauss point, the constitutive update reads:

Algorithm 1. Local constitutive update at a Gauss point

Require: FE solution at $t_{n+1}$ , state variables at $t_n$ , time step $\mathrm{\Delta }t$ Ensure: Updated stresses $(\tilde{\mathit{\sigma }},\tilde{\mathbf{\chi }},\tilde{\mathbf{\Gamma }})$ , internal variables, and algorithmic tangent ${\mathbb{C}}_{\mathrm{alg}}$  1: Compute trial state: compute the trial elastic strain (or deformation gradient) from the FE solution.  2: Slip integration: evaluate ${\tau }_{\mathrm{eff}}$ and integrate the slip rate (113) by backward Euler to obtain ${\gamma }_{n+1}$ .  3: SSD update: update ${\rho }_{\mathrm{SS}}$ using (115).  4: Kinematics update: update ${\mathbf{R}}^e$ , ${\tilde{\mathbf{\kappa }}}^e$ , and $\mathit{\theta }$ via (116)–(117).  5: Return to FE: compute and return $\tilde{\mathit{\sigma }}$ , $\tilde{\mathbf{\chi }}$ , $\tilde{\mathbf{\Gamma }}$ , and the algorithmic tangent ${\mathbb{C}}_{\mathrm{alg}}$ to the global FE solver.

Remark

This local constitutive structure is modular and directly compatible with the mixed finite-element formulation described in Section 9.2. It provides a clear separation between the integration-point physics (slip, hardening, curvature) and the global solution of the coupled macro–micro balance equations, facilitating both validation studies and future extensions to multiple slip systems and finite-strain kinematics.

Finite element implementation

This section summarizes the finite-element (FE) implementation adopted to solve the coupled macro–micro boundary value problem associated with the proposed dual-length-scale gradient crystal plasticity model. The objective is to provide a reproducible computational procedure suitable for validation studies and for subsequent coupling with microstructure-evolution models.

Mixed formulation and primary unknowns

The free-energy density contains a curvature contribution and a curvature-gradient contribution through the disclination density:

$$\theta =\kappa ,$$

(118)

which would introduce second spatial derivatives if $\kappa$ were eliminated in favor of an underlying rotation field. To avoid $C^1$ interpolation, a mixed $C^0$ formulation is adopted in which the curvature tensor $\kappa$ is treated as an independent FE field in addition to the displacement field.

In the two-dimensional plane-strain setting used for the illustrative examples, the global unknowns are:

$$\mathbf{u}=(u_1,u_2),\kappa =\left(\begin{array}{cc}{\kappa }_{11}\hfill & {\kappa }_{12}\hfill \\ \multicolumn{1}{c}{{\kappa }_{21}}& {\kappa }_{22}\hfill \end{array}\right),$$

(119)

while all internal variables associated with crystal plasticity (slip, hardening, etc.) are stored and updated at Gauss points.

Spatial discretization

The computational domain $\mathrm{\Omega }\subset {\mathbb{R}}^2$ is discretized by isoparametric bilinear quadrilateral elements (Q4). Standard $C^0$ Lagrange interpolation is used for both fields:

$${\mathbf{u}}_h\left(\mathbf{x}\right)=\sum _{a=1}^{n_n}{N}_a\left(\mathbf{x}\right){\mathbf{u}}_a,{\kappa }_h\left(\mathbf{x}\right)=\sum _{a=1}^{n_n}{N}_a\left(\mathbf{x}\right){\kappa }_a,$$

(120)

where $N_a$ are the usual Q4 shape functions, $n_n$ is the number of nodes per element, and $\{{\mathbf{u}}_a,{\kappa }_a\}$ denote nodal degrees of freedom. Element integrals are evaluated with a $2\times 2$ Gauss rule.

For plane strain, the standard small-strain kinematics are employed in the prototype implementation:

$$\epsilon \left(\mathbf{u}\right)=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\top }),$$

(121)

and the Cauchy stress $\sigma$ is obtained from the algorithmic crystal plasticity update described in Section 9.2.5.

Discrete Curl operator in two dimensions

The curvature-gradient term is driven by the row-wise Curl of $\kappa$ . In 2D plane strain, define $\theta \in {\mathbb{R}}^2$ by:

$${\theta }_1={\partial }_x{\kappa }_{12}-{\partial }_y{\kappa }_{11},{\theta }_2={\partial }_x{\kappa }_{22}-{\partial }_y{\kappa }_{21}.$$

(122)

With $C^0$ interpolation for ${\kappa }_h$ , ${\theta }_h$ involves only first derivatives of the FE shape functions and is therefore straightforward to assemble element-wise.

Weak form and residual statements

The discrete problem is written in residual form. Let $\delta \mathbf{u}$ and $\delta \kappa$ be admissible test functions. The macroforce balance is enforced by the standard virtual-work statement:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{R}}_u({\mathbf{u}}_h,{\kappa }_h;\delta \mathbf{u}):={\int }_{\mathrm{\Omega }}\sigma :{\nabla }^s\delta \mathbf{u}\mathrm{d}V-{\int }_{\mathrm{\Omega }}\mathbf{b}·\delta \mathbf{u}\mathrm{d}V}\\ \multicolumn{1}{c}{-{\int }_{\partial {\mathrm{\Omega }}_t}\overline{\mathbf{t}}·\delta \mathbf{u}\mathrm{d}S=0.}\end{array}$$

(123)

The curvature field is governed by the energetic microforce balance associated with the curvature part of the free energy,

$${\psi }_{\kappa }(\kappa ,\theta )=\frac{1}{2}\mu {\ell }_1^2\parallel \kappa {\parallel }^2+\frac{1}{2}\mu {\ell }_2^4\parallel \theta {\parallel }^2,\theta =\kappa ,$$

(124)

leading to the weak form:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{R}}_{\kappa }({\mathbf{u}}_h,{\kappa }_h;\delta \kappa ):={\int }_{\mathrm{\Omega }}\mu {\ell }_1^2\kappa :\delta \kappa \mathrm{d}V+{\int }_{\mathrm{\Omega }}\mu {\ell }_2^4\left(\kappa \right)·\left(\delta \kappa \right)\mathrm{d}V}\\ \multicolumn{1}{c}{-{\int }_{\mathrm{\Omega }}S:\delta \kappa \mathrm{d}V-{\int }_{\partial {\mathrm{\Omega }}_q}\overline{q}:\delta \kappa \mathrm{d}S=0.}\end{array}$$

(125)

Here $S$ is a source term that couples the curvature field to nonuniform plastic flow. In the single-slip plane-strain benchmarks, a parsimonious coupling is used, where only the component conjugate to the dominant shear mechanism is driven by the slip gradient. Specifically, for a slip system with slip direction $s$ and slip-plane normal $m$ , the scalar slip $\gamma$ is recovered nodally from Gauss-point values and the source is taken as:

$$S_{12}=c_{\kappa }{\partial }_y\gamma ,S_{\mathrm{ij}}=0\mathrm{otherwise},$$

(126)

with a coupling parameter $c_{\kappa }>0$ . This choice is sufficient to activate the length-scale effects through ${\ell }_1$ and ${\ell }_2$ in the validation tests.

Constitutive update and internal variables

Crystal plasticity is integrated at Gauss points through an operator-split algorithm. For each time increment $t_n\to t_{n+1}$ , the FE solution provides the strain (or deformation gradient in a finite-strain extension), which is passed to a local return/update procedure that evolves: (i) the slip $\gamma$ (single slip in the present benchmarks), (ii) the isotropic hardening variables, and (iii) any additional state variables needed by the model.

A rate-dependent flow rule is used in the implementation to regularize the local problem and ensure robust convergence:

$$\stackrel{·}{\gamma }={\stackrel{·}{\gamma }}_0{\left(\frac{\left|{\tau }_{\mathrm{eff}}\right|}{{\tau }_c}\right)}^{1/m}\mathrm{sign}\left({\tau }_{\mathrm{eff}}\right),$$

(127)

where ${\stackrel{·}{\gamma }}_0$ is a reference slip rate, $m$ is the rate-sensitivity exponent, and ${\tau }_c$ is the current slip resistance (including hardening). The effective resolved shear stress is corrected by the curvature-driven backstress:

$${\tau }_{\mathrm{eff}}=\tau -{\chi }_s,{\chi }_s={\alpha }_{\chi }\mu {\ell }_1^2{\kappa }_{12},$$

(128)

where ${\alpha }_{\chi }$ is a dimensionless coupling coefficient. The local update is performed by backward Euler integration of equation (127) with a small Newton solve for ${\gamma }_{n+1}$ , after which the stress and algorithmic tangent are returned to the global assembly.

Solution strategy and Newton iterations

The coupled discrete problem is solved by Newton iterations on the global unknown vector. For the Q4 plane-strain implementation, each node carries:

$$\underset{2\mathrm{DOFs}}{\underbrace{(u_1,u_2)}}\oplus \underset{4\mathrm{DOFs}}{\underbrace{({\kappa }_{11},{\kappa }_{12},{\kappa }_{21},{\kappa }_{22})}},$$

(129)

leading to a block-structured linearized system. In the baseline implementation, a staggered quasi-Newton approach is employed for robustness:

• 1. Update the crystal plasticity variables locally at Gauss points given the current ${\mathbf{u}}_{\mathit{h}}$ and ${\kappa }_{\mathit{h}}$ .

• 2. Assemble and solve the macro problem for $\mathbf{\Delta }\mathbf{u}$ .

• 3. Assemble and solve the curvature problem for $\mathrm{\Delta }\kappa$ using the updated nodal slip-gradient source (equation (126)).

This strategy yields stable convergence for the benchmark problems and is sufficient for the validation section. A fully monolithic Newton method can be obtained by retaining the off-diagonal Jacobian blocks resulting from the dependence of $\sigma$ on $\kappa$ (through ${\chi }_s$ ) and of $S$ on the plastic variables.

Boundary conditions and benchmark loading

Standard Dirichlet and Neumann boundary conditions are prescribed on the displacement field. For the curvature field, natural boundary conditions emerge from equation (125) via the integration by parts of the Curl term and are interpreted as micro-tractions $\overline{q}$ on $\partial {\mathrm{\Omega }}_q$ . In the numerical examples, both micro-free (natural) and micro-hard (Dirichlet) conditions can be considered; the reported results employ the boundary conditions stated explicitly in the corresponding benchmark subsections.

Remark

The present implementation summary is intended as a computationally tractable realization of the proposed theoretical framework. The mixed formulation retains the dual internal length scales $({\ell }_1,{\ell }_2)$ and enables systematic comparisons with classical (no-length-scale) crystal plasticity and with single-length-scale gradient variants in the validation studies.

Numerical results: Micro-bending and size effects

This section presents numerical illustrations of the proposed gradient-enhanced crystal plasticity model with two internal length scales. The objective is twofold: (i) to demonstrate the well-posedness and numerical robustness of the formulation, and (ii) to highlight its added value compared to classical crystal plasticity through size-dependent responses in micro-bending.

Problem setup and numerical implementation

A two-dimensional cantilever beam of length $L_x=1$ and thickness $L_y$ , modeled under plane strain conditions is considered. The beam is discretized using structured bilinear quadrilateral (Q4) finite elements. Each node carries the standard displacement degrees of freedom $(u_1,u_2)$ as well as the components of the curvature tensor $\mathbf{\kappa }=\{{\kappa }_{11},{\kappa }_{12},{\kappa }_{21},{\kappa }_{22}\}$ .

The left edge ( $x=0$ ) is fully clamped, while the right edge ( $x=L_x$ ) is subjected to a downward point load applied at the top corner, inducing a bending-dominated deformation. This loading configuration is deliberately chosen to activate curvature gradients and to prevent trivial relaxation of the plastic slip through homogeneous shear.

The constitutive response follows a single-slip viscoplastic crystal plasticity model, augmented by a curvature-dependent backstress of the form:

$${\chi }_s={\alpha }_{\chi }\mu {\ell }_1^2{\kappa }_{12},$$

(130)

and governed by a curvature evolution equation involving both a mass-like term ${\ell }_1^2\mathbf{\kappa }$ and a Curl-regularization term ${\ell }_2^4\mathrm{Curl}\mathbf{\kappa }$ . Natural boundary conditions are imposed on $\mathbf{\kappa }$ , unless otherwise stated.

The coupled problem is solved using a staggered Newton scheme, where the macroscopic equilibrium and the curvature field are updated iteratively, while the local slip variables are integrated at the Gauss-point level using an implicit viscoplastic update. Unless otherwise specified, the simulations reported below use identical material parameters and differ only in the values of the internal length scales.

Reaction–deflection response

Figure 1 shows the load–deflection curves obtained from the cantilever micro-bending test for three models: (i) classical crystal plasticity, obtained by suppressing the curvature field ( $\mathbf{\kappa }\equiv 0$ ), (ii) a one-length formulation with ${\ell }_1>0$ and ${\ell }_2=0$ , and (iii) the proposed two-length formulation with ${\ell }_1>0$ and ${\ell }_2>0$ .

Figure 1.

Load–deflection response for the cantilever micro-bending test. Comparison between classical crystal plasticity and gradient-enhanced formulations.

All models exhibit a nonlinear response associated with the activation of viscoplastic slip. However, the gradient-enhanced models display a systematically higher effective stiffness compared to the classical formulation. This behavior reflects the additional energetic cost associated with curvature development and the resulting backstress opposing plastic flow.

Size effect induced by internal length scales

To assess the size sensitivity of the model, the micro-bending test is repeated for three beam thicknesses, $L_y=\{0.20,0.10,0.05\}$ , while keeping all material parameters unchanged. Figure 2 reports the normalized response in terms of the dimensionless quantities:

$$\frac{-u_2^{\mathrm{tip}}}{L_y}\mathrm{and}\frac{P}{L_y^2},$$

(131)

which would collapse onto a single curve for a classical scale-free continuum.

Figure 2.

Normalized load–deflection response under tip point loading. The two-length model exhibits a clear thickness-dependent response, in contrast to the classical crystal plasticity formulation.

The classical crystal plasticity model shows an almost perfect collapse of the normalized curves, indicating the absence of intrinsic size effects. In contrast, the gradient-enhanced model exhibits a clear separation between the responses corresponding to different thicknesses. Thinner beams display an increased resistance to bending, consistent with the physical interpretation of curvature-induced hardening.

This size effect originates from the competition between the internal length scales ${\ell }_1$ and ${\ell }_2$ and the characteristic structural dimension $L_y$ . When $L_y$ becomes comparable to the internal lengths, curvature gradients become energetically significant, leading to a non-classical strengthening mechanism.

Role of the second internal length scale

A key feature of the proposed formulation is the introduction of a second internal length scale ${\ell }_2$ through the Curl-regularization of the curvature tensor. While the length ${\ell }_1$ primarily controls the magnitude of the curvature-induced backstress, the parameter ${\ell }_2$ governs the spatial regularity of the curvature field.

Numerically, this manifests as smoother curvature distributions and a reduced localization of $\mathbf{\kappa }$ in regions of high bending gradients. The resulting response differs quantitatively from the one-length model, particularly for small thicknesses, thereby demonstrating that ${\ell }_2$ is not a redundant parameter but contributes independently to the macroscopic behavior.

Discussion and implications

These numerical results highlight the added value of the proposed two-length crystal plasticity model. Unlike classical formulations, the model naturally captures size effects in bending-dominated problems without introducing ad hoc boundary conditions or nonlocal averaging procedures.

The curvature-based framework provides a physically transparent mechanism for size dependence, rooted in the energetics of lattice curvature and dislocation structures. Moreover, the formulation is compatible with multiscale extensions, such as coupling with phase-field descriptions of microstructure evolution or discrete dislocation dynamics, while retaining a well-defined continuum structure.

Overall, the micro-bending simulations demonstrate that the proposed model is both computationally tractable and capable of reproducing key non-classical features observed in small-scale plasticity.

Comparative discussion with strain-gradient and dislocation-based crystal plasticity

The micro-bending simulations reported in Section 10 exhibit a thickness-dependent strengthening, that is, absent in classical (scale-free) crystal plasticity. This behavior is consistent with a large body of experimental and theoretical evidence on “smaller is stronger” effects in bending, torsion, and indentation at the micron scale. In this subsection, the present two-length curvature-based formulation is positioned relative to canonical strain-gradient plasticity (SGP) and gradient crystal plasticity (GCP) frameworks.

Comparison with phenomenological SGP (Fleck–Hutchinson)

The Fleck–Hutchinson family of phenomenological strain-gradient plasticity theories introduces additional hardening through dependence on plastic strain gradients, typically through an effective measure combining classical plastic strain and one or more invariants of $\nabla {\mathbf{\epsilon }}^p$ .^21,22 In bending-dominated problems, these theories predict increased resistance as the structural dimension decreases, due to the energetic or dissipative cost associated with nonuniform plastic deformation. The present model reproduces the same qualitative signature in micro-bending; however, the mechanism is formulated in terms of a curvature tensor $\mathbf{\kappa }$ and its Curl, rather than directly through $\nabla {\mathbf{\epsilon }}^p$ . This distinction is important for two reasons.

First, the curvature variable provides a more direct kinematic proxy for geometrically necessary dislocations (GNDs) and lattice curvature effects, which are known to dominate size dependence in single crystals. Second, the model naturally yields a two-parameter regularization: ${\ell }_1$ controls the amplitude of curvature-induced backstress, while ${\ell }_2$ controls the spatial regularity of the curvature field. In contrast, many phenomenological SGP formulations employ a single dominant length scale, with the regularity primarily governed by the choice of gradient invariants and the boundary conditions imposed on higher-order microtractions.

Comparison with gradient crystal plasticity of Gurtin–Anand type

The Gurtin–Anand framework introduces gradients of slip (or slip-system hardening variables) and associated microforce balances, producing higher-order stresses conjugate to slip gradients.^28,29 This class of theories is physically attractive because it is formulated at the slip-system level and provides a thermodynamically consistent microforce structure, including boundary conditions that can model micro-free or micro-hard surfaces. The present formulation shares the same thermodynamic philosophy: an energetic penalty is associated with nonuniform internal fields and gives rise to additional stresses (backstresses) that oppose slip.

The key difference lies in the choice of the gradient variable and the resulting boundary data. In Gurtin–Anand-type GCP, the gradient typically acts on slip ${\gamma }^{\alpha }$ (or related measures), so the length scale enters through terms like ${\ell }^2{|\nabla {\gamma }^{\alpha }|}^2$ and the associated higher-order stresses produce boundary conditions on ${\gamma }^{\alpha }$ (micro-hard) or on microtractions (micro-free). In the present two-length model, the regularization is expressed through curvature: a mass-like penalty ${\ell }_1^2{\left|\mathbf{\kappa }\right|}^2$ and a Curl penalty ${\ell }_2^4{\left|\mathrm{Curl}\mathbf{\kappa }\right|}^2$ . Consequently, natural boundary conditions correspond to vanishing higher-order fluxes associated with $\mathbf{\kappa }$ , which can be interpreted as a curvature microtraction condition. This provides a direct route to encode surface constraints on lattice curvature, which is particularly relevant in micro-bending where curvature gradients concentrate near clamped boundaries and load application zones.

Comparison with dislocation-density and Nye-tensor based approaches

A substantial body of GCP work defines the defect content through a dislocation density tensor (often related to the Nye tensor), built from gradients of plastic distortion, and uses energetic penalties in terms of this density.^5,7,35 In many such formulations, the internal length scale emerges through the energetic cost of GND density, yielding size effects in torsion and bending. The present model is compatible with this spirit: the curvature tensor $\mathbf{\kappa }$ and its Curl are intended to capture the influence of dislocation/curvature mechanisms on plastic response. The two length scales can be interpreted as separating contributions associated with (i) the local curvature magnitude (through ${\ell }_1$ ) and (ii) the spatial organization or interaction range of curvature gradients (through ${\ell }_2$ ). Numerically, this separation is observed in the difference between one-length and two-length responses in micro-bending, especially at small thickness where curvature localization would otherwise be more pronounced.

Implications for micro-bending and “added value” of the two-length structure

From a modeling perspective, the two-length model refines the description of size effects by decoupling two roles often merged in single-length formulations: (a) strengthening through an internal backstress proportional to curvature (amplitude control) and (b) regularization through Curl-type smoothing (spatial structure control). This is analogous in spirit to introducing both energetic and higher-gradient regularization channels, but here achieved within a curvature-based architecture tailored to crystal plasticity.

In the micro-bending simulations, the classical CP response collapses under normalization, as expected for a scale-free theory, whereas the gradient-enhanced response exhibits a thickness dependence. This qualitative signature aligns with the predictive trends of Fleck–Hutchinson SGP and Gurtin–Anand GCP, while offering an alternative route based on curvature kinematics, that is, well suited for coupling with defect-transport or phase-field microstructure models.

Limitations and outlook

The present numerical examples are illustrative and intended to validate implementation and demonstrate qualitative size dependence. A quantitative validation against experimental micro-bending data and/or benchmark simulations from established GCP models (e.g. slip-gradient formulations) is the natural next step. Such a study will also enable calibration of ${\ell }_1$ and ${\ell }_2$ and will clarify the regimes in which the two-length separation yields a measurable advantage over single-length theories. Such studies are conducted in the section that follows.

Numerical results: Numerical simulations versus Demir et al. ¹⁵ ’ experiments

In this section, the numerical micro-bending response predicted by the present finite-element implementation is compared with the cantilever micro-bending experiments reported by Demir et al. ¹⁵ (see Figure 3) on single-crystal copper micro-beams. The experiments were performed at three representative thickness levels, ∼ $t\approx 1.0$ , $1.7$ , and $4.2\mu \mathrm{m}$ and reported as force–displacement curves exhibiting a pronounced size effect and serrated plasticity.

Figure 3.

Experimental micro-bending force–displacement curves for single-crystal copper micro-beams reported by Demir et al. ¹⁵ (see Figure 3 therein). The results exhibit a pronounced thickness dependence and an intermittent serrated response. The figure is adapted from Demir et al. ¹⁵ These experiments provide a qualitative reference for the size-effect trend captured by the numerical simulations in Figure 4.

Numerical set-up and force definition

The simulations consider a 2D plane-strain cantilever bending configuration discretized with bilinear Q4 elements and a single-slip viscoplastic crystal plasticity update. The vertical tip deflection $\delta$ is prescribed on the free end (right edge), and the reaction measure $P\left(\delta \right)$ is computed as the sum of the nodal reactions on the loaded edge vertical degrees of freedom. Since the finite-element problem is formulated in 2D, the computed reaction is a force per unit out-of-plane thickness. For comparison with the experimental total force $F$ , one may convert the numerical reaction via:

$$F_{\mathrm{num}}\left(\delta \right)=P_{2\mathrm{D}}\left(\delta \right)w,$$

(132)

where $w$ is the micro-beam width reported in Demir et al. ¹⁵ (Table 1). This mapping is standard for plane-strain idealizations and allows the numerical curves to be placed on the same force scale as the experiments, provided that the width used in the conversion matches the tested specimen.

Two constitutive descriptions are contrasted: (i) a classical crystal plasticity model (no intrinsic length scale) and (ii) the proposed two-length formulation. In the latter, ${\ell }_1$ primarily controls the magnitude of the curvature-induced backstress entering the slip driving force, whereas ${\ell }_2$ governs the spatial regularity of the curvature field through Curl-regularization.

Force–deflection response and thickness dependence

Figure 4 reports the computed reaction–deflection curves for the three thicknesses, comparing the classical CP and the two-length model. The simulations exhibit a clear thickness dependence: thinner beams display a higher apparent bending resistance, that is, a size effect. Moreover, the two-length model predicts a systematically stronger size effect than classical CP, which is consistent with the additional energetic penalty associated with curvature gradients and with the strengthening mechanisms underlying strain-gradient-type theories.

Figure 4.

Numerical micro-bending response: reaction force $P$ versus imposed tip deflection $\delta$ for three thicknesses ( $t\approx 1.0,1.7,4.2\mu$ m), comparing classical crystal plasticity (classic CP) and the proposed two-length model (two-length). In a plane-strain idealization, the total force can be obtained as $F_{\mathrm{num}}=P_{2\mathrm{D}}w$ using the specimen width $w$ from Demir et al. ¹⁵ The two-length formulation predicts a stronger size effect and higher bending resistance with decreasing thickness.

Figure 3 reproduces the experimental force–displacement curves of Demir et al. The experiments show (i) a strong thickness dependence and (ii) intermittent force drops/serrations, typically attributed to discrete dislocation activity and localized plastic events at the microscale. While the present single-slip, monotonic-loading prototype does not aim to reproduce the full hysteretic and serrated signature quantitatively, it captures the primary experimental trend: decreasing thickness leads to an increased bending resistance.

Scope and limitations of the validation

The present comparison is intended as a first validation in the sense of: (i) verifying that the proposed two-length model reproduces the experimentally observed trend of micro-bending strengthening with decreasing thickness and (ii) demonstrating that the additional length scale ${\ell }_2$ contributes non-redundantly to the global response by modifying both the magnitude and regularity of the curvature field, hence affecting the macroscopic bending curve. A fully quantitative match of the complete experimental loops in Demir et al. ¹⁵ Figure 3 (including unloading/reloading and serrations) is outside the scope of this initial numerical study and would require additional modeling ingredients, for example, multi-slip activation, cyclic loading capability, and/or stochastic/discrete dislocation effects to reproduce the intermittent plasticity observed experimentally.

Conclusion

This paper has developed a gradient-enhanced finite-strain crystal plasticity framework that introduces two intrinsic length scales to represent, within a single continuum theory, both lattice-curvature effects and curvature-gradient effects. The primary objective has been to advance the predictive modeling of size-dependent plastic deformation in crystalline materials by embedding defect-sensitive mechanisms directly into the kinematic, energetic, and balance-law structure of the model. The resulting formulation goes beyond classical, scale-free crystal plasticity by providing a transparent separation between (i) curvature-amplitude strengthening associated with dislocation-type incompatibilities and (ii) curvature-gradient regularization associated with higher-order incompatibilities of disclination type. This separation is not only conceptually important; it is also essential for capturing the distinct spatial patterns and boundary-layer phenomena observed in small-scale tests and in boundary-constrained crystalline deformation.

A first major contribution of the work is the systematic organization of the theoretical framework, from kinematics to constitutive structure to balance laws. The paper establishes the defect-sensitive measures in a finite-strain setting and clarifies their physical meaning through a differential-geometric viewpoint. This perspective provides a coherent interpretation of the additional fields and operators that appear in the theory, and it shows how the model naturally connects to established notions from generalized continua and curvature-based plasticity. Within this kinematic setting, a thermodynamically consistent internal-variable structure is then posed to define the material response. The free-energy statement is used as the unifying mechanism that generates the full set of work-conjugate stress measures, including both classical stresses and higher-order stresses associated with curvature and curvature gradients. As a consequence, the balance structure of the theory arises in a systematic way and yields higher-order equilibrium relations that are tightly linked to the defect measures. This approach ensures that the model is not an ad hoc regularization but a constitutive framework grounded in a clear energetic interpretation.

A second major contribution is the explicit coupling between the curvature-driven fields and crystallographic slip. The paper specifies how defect-sensitive internal stresses enter the driving force for slip and how the evolution of internal variables is linked to defect content and plastic flow. This coupling is essential for translating the geometric content of the theory into observable macroscopic effects such as strengthening, size dependence, and boundary-layer formation. In particular, by distinguishing the roles of the two length scales, the framework provides a direct means to control and interpret (i) the amplitude of curvature-related internal resistance and (ii) the spatial organization and smoothing of curvature gradients. This feature is important because many size effects in crystals are governed not only by the magnitude of incompatibility but also by how it localizes near boundaries and interfaces.

The analytical studies developed in the paper provide direct validation and physical insight. The simple-shear single-crystal strip problem serves as a benchmark demonstrating that the model produces enhanced hardening under constrained deformation, together with the emergence of boundary layers that localize orientation gradients and defect-related fields. This canonical setting clarifies the role of curvature-induced internal stresses in generating a backstress-like resistance that increases with decreasing specimen dimension. Beyond this primary validation, the analytical program is extended to additional canonical boundary-value problems—torsion of a cylindrical bar, bending of a crystal plate, and a bicrystal tilt boundary. These problems are not merely illustrative; collectively, they probe the theory across distinct classes of deformation modes and boundary conditions. Torsion highlights the interaction between rotation gradients and geometry-induced incompatibilities; bending reveals the localization of orientation gradients through thickness and the emergence of characteristic boundary-layer structures; and the bicrystal tilt boundary provides a structurally meaningful setting in which incompatibility localizes at an interface, thereby emphasizing the relevance of higher-order fields in modeling orientation mismatch. Across all these cases, the analysis underscores how the two length scales govern boundary-layer thickness, orientation-gradient localization, and the distribution of defect-related internal stresses.

A third major contribution is the establishment of a practical computational route and its use to produce structural-scale simulations that complement the analytical benchmarks. The paper describes a finite-element implementation strategy compatible with the higher-order structure of the governing equations, including the local integration-point constitutive update and the mixed weak form required to represent curvature-related fields. This implementation discussion is central to the utility of the theory: the presence of higher-order stresses and additional balance relations places nontrivial demands on discretization, linearization, and solution strategy, and the paper’s algorithmic presentation provides a clear pathway for reproducing and extending the results. The numerical studies reported for cantilever micro-bending quantify thickness-dependent strengthening and demonstrate how the curvature-gradient length scale influences not only the overall reaction–deflection behavior but also the spatial organization of curvature-related fields. In this sense, the computational results serve a dual purpose: they verify that the formulation can be implemented robustly in a finite-element framework, and they show how the two-length structure manifests in a representative structural configuration. The comparison with reported micro-beam experiments on single-crystal copper further supports the relevance of the proposed framework for capturing experimentally observed thickness dependence in micro-bending.

Taken together, the theoretical development, analytical validations, and numerical simulations establish a comprehensive two-length-scale gradient crystal plasticity framework that provides a defect-informed continuum route to size-dependent strengthening and localization control. The theory unifies dislocation-type and disclination-type mechanisms within a single energetic setting and clarifies their distinct contributions through both canonical closed-form reductions and structural-scale computations. In doing so, the paper provides a consistent bridge between defect geometry and macroscopic response, while also delivering an implementation-ready formulation for continued exploration in computational mechanics.

Several immediate extensions follow naturally from the present results. On the analytical side, the established approach can be further applied to additional boundary-value problems that probe different combinations of constraint, loading mode, and orientation structure, including tension/compression, multiaxial loading, and mixed boundary conditions. On the computational side, the finite-element route opens the door to simulations of more complex specimen geometries and to studies of localization and pattern formation in settings where boundary layers interact with structural features. In addition, the two-length-scale structure offers a natural basis for systematic parameter identification using size-effect data across multiple tests, since the two length scales can be linked to distinct signatures in boundary-layer thickness and curvature-gradient localization. Finally, because the framework is posed in a general internal-variable and energetic format, it provides a clear foundation for coupling with additional physics (e.g. microstructural evolution mechanisms) whenever such couplings are of interest.

Overall, this work provides a robust theoretical and computational foundation for gradient-enhanced crystal plasticity with two intrinsic length scales. By integrating curvature-based strengthening and curvature-gradient regularization within a finite-strain, thermodynamically consistent framework and by supporting the formulation with analytical benchmarks, a detailed implementation strategy, and micro-bending simulations, the paper advances a mechanism-based continuum route for predicting crystalline deformation and hardening at reduced length scales.