# Strain Localization Modes within Single Crystals Using Finite Deformation Strain Gradient Crystal Plasticity

Lei Cai, Mohamed Jebahi, Farid Abed-Meraim

► To cite this version:

Lei Cai, Mohamed Jebahi, Farid Abed-Meraim. Strain Localization Modes within Single Crystals Using Finite Deformation Strain Gradient Crystal Plasticity. Crystals, 2021, 11 (10), pp.1235. 10.3390/cryst11101235. hal-03439390

HAL Id: hal-03439390

https://hal.science/hal-03439390

Submitted on 22 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Strain Localization Modes within Single Crystals Using Finite Deformation Strain Gradient Crystal Plasticity

Lei Cai ID, Mohamed Jebahi *ID and Farid Abed-Meraim ID

Arts et Metiers Institute of Technology, CNRS, Université de Lorraine, LEM3, F-57000 Metz, France;
lei.CAI@ensam.eu (L.C.); Farid.ABEDMERAIM@ensam.eu (F.A.-M.)

* Correspondence: mohamed.jebahi@ensam.eu

**Abstract:** The present paper aims at providing a comprehensive investigation of the abilities and limitations of strain gradient crystal plasticity (SGCP) theories in capturing different kinds of localization modes in single crystals. To this end, the small deformation Gurtin-type SGCP model recently proposed by the authors, based on non-quadratic defect energy and the uncoupled dissipation assumption, is extended to finite deformation. The extended model is then applied to simulate several single crystal localization problems with different slip system configurations. These configurations are chosen in such a way as to obtain idealized slip and kink bands as well as general localization bands, i.e., with no particular orientation with respect to the initial crystallographic directions. The obtained results show the good abilities of the applied model in regularizing various kinds of localization bands, except for idealized slip bands. Finally, the model is applied to reproduce the complex localization behavior of single crystals undergoing single slip, where competition between kink and slip bands can take place. Both higher-order energetic and dissipative effects are considered in this investigation. For both effects, mesh-independent results are obtained, proving the good capabilities of SGCP theories in regularizing complex localization behaviors. The results associated with higher-order energetic effects are in close agreement with those obtained using a micromorphic crystal plasticity approach. Higher-order dissipative effects led to different results with dominant slip banding.

**Keywords:** strain gradient crystal plasticity; finite deformation; strain localization modes; slip bands; kink bands; general localization bands

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Cai, L.; Jebahi, M.; Abed-Meraim, F. Strain Localization Modes within Single Crystals Using Finite Deformation Strain Gradient Crystal Plasticity. *Crystals* **2021**, *11*, 1235. https://doi.org/10.3390/cryst11101235

Academic Editor: Pavel Lukáč

Received: 27 August 2021

Accepted: 6 October 2021

Published: 13 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Strain localization is an important plastic instability process occurring prior to fracture. It is usually observed in the form of narrow bands of intense plastic shear strain in deformed bodies undergoing severe inhomogeneous deformation. This process has received strong scientific interest since the pioneering work of Considère [1], and numerous studies have been published on the subject [2–9]. Several types of localization modes have been revealed by these studies, depending on the microstructure of materials. Considering a single crystal with multiple activated slip systems, it has been shown that strain localization may occur along general shear bands, which can have an arbitrary orientation with respect to the crystallographic directions [3,10,11].

When only one slip system is activated, two types of shear bands, known as slip and kink bands, may occur at incipient plasticity according to the seminal work of Asaro and Rice [4] bifurcation analysis. Slip bands are sharp structures formed by intense dislocation slip on a few neighboring crystallographic planes. They are parallel to the active slip plane and correspond to the main signature of plastic slip in single crystals. Kink bands, which were first observed by Orowan [12], are localization bands orthogonal to the slip direction. They are accompanied by a strong lattice rotation in the bands compared to the surrounding regions. These bands are bounded by high lattice curvature boundaries formed by geometrically necessary dislocations (GNDs) walls.
These bands present a much larger width than slip bands. They are known to occur when strain incompatibility arises and there are not enough slip systems available to accommodate this incompatibility. According to the experimental evidence of Di Gioacchino and Quinta da Fonseca [13], such bands are characterized by the presence of a dense distribution of short slip bands along the slip plane normal direction. An interesting discussion about the occurrence of these bands can be found in [9].

Conventional crystal plasticity (CCP) theories have widely been applied in the literature to study strain localization within single crystals [5,14–17]. Although these theories are able to capture several kinds of localization modes, including slip and kink bands, they present mesh-dependence difficulties. Localization bands obtained using such theories present mesh-dependent widths and plastic slip peak values [6,18]. In addition, CCP theories identically predict slip and kink bands, which appear in this framework as equivalent bifurcation modes.

This result was demonstrated by Asaro and Rice [4] using bifurcation analysis for single crystals undergoing single slip. A confirmation of such results for poly-crystals was recently provided by Marano et al. [8] using Fast Fourier Transform (FFT) poly-crystal simulations. In these simulations, systematically comparable distributions of slip and kink bands were obtained. This seems not to be in accordance with strain localization in real materials where overwhelming dominant slip bands are commonly observed. In the context of CCP theories, selection between slip and kink bands is only due to structural effects.

From a physical point of view, the equivalence between these bands is questioned as they rely on quite different formation processes. Compared to slip bands, more dislocation sources are generally required to form kink bands. Consequently, CCP theories are not suitable to study localization phenomena in single crystals. These theories include no internal length scale(s) allowing for stabilizing localization, which theoretically will occur in a set of zero measures.

A solution to overcome the limitations of the aforementioned theories consists in applying nonlocal plasticity approaches. Including internal length scale(s), these approaches provide a natural framework to capture nonlocal effects. Their promising features have motivated the development of various nonlocal models for both single and poly-crystal structures over the last two decades [6,8,19–28]. Although most of these models have been developed to investigate the size effect at small scales, a few works applying them to localization problems exist in the literature.

In the context of single crystals, Forest [6] proposed a small deformation Cosserat crystal plasticity model to investigate the formation of slip and kink bands. Taking into account lattice rotation, the proposed model allows for breaking the equivalence between the two bands. Furthermore, it enables the regularization of kink bands, which present a finite width. However, this is not the case for slip bands, which are not necessarily accompanied by lattice curvature. Recently, a reduced micromorphic crystal plasticity model relying on the gradient of a cumulative scalar slip variable was proposed by Ling et al. [29].

This model presents the unique feature of regularizing both slip and kink bands. However, the same internal length scale is involved in the regularization of these bands, which could be problematic considering their different physical features and formation processes. The application of the proposed model to study strain localization within single crystals undergoing single slip revealed a complex localization behavior with competition between kink and slip bands. Experimental evidence of this behavior can be found in [17].

Another class of nonlocal approaches, which presents auspicious features to capture localization phenomena in single crystals, is the class of strain gradient crystal plasticity (SGCP) theories. This class has been the subject of a large number of recent works mostly focusing on size effects [11, 24, 30–34]. However, only a few works applying SGCP theories to study localization phenomena can be found in the literature. Using a finite deformation SGCP model, Kuroda [11] investigated the effects of the finite element type on the formation of localization bands within single crystals. Scherer et al. [35] proposed a reduced SGCP
model based on the work of Ling et al. [29] to study the formation of localization bands in voided irradiated materials.

This model relies on the gradient of a cumulative scalar slip variable, and this wa implemented using an approach combining micromorophic and Lagrange multiplier methods. Recently, Marano et al. [9] applied a small deformation version of the model of Gurtin [21], based on the full GND density tensor, to study the influence of the strain gradient plasticity on the modeling of localization bands within poly-crystal solids. In almost all existing studies of localization phenomena in single crystals, only higher-order energetic effects have been considered. Higher-order dissipative effects on these phenomena, particularly the competition between slip and kink bands, have not yet been explored.

The present paper aims at tackling these tasks by providing a comprehensive investigation of the abilities and limitations of SGCP theories in capturing different localization modes while considering both higher-order energetic and dissipative effects. To this end, the small deformation Gurtin-type SGCP model presented in [34], based on non-quadratic defect energy and uncoupled dissipation assumption, is extended to a finite deformation framework following the procedure detained in [21,23].

The extended model is then applied to simulate various localization problems implying single crystals with different slip system configurations. A focus is given to the ability of this model to capture the complex localization behavior characterizing single crystals undergoing single slip. The effects of higher-order energetic and dissipative stresses on this behavior are discussed.

Following this introduction, the paper is organized as follows. Section 2 details the theoretical formulation of the proposed Gurtin-type finite deformation SGCP model. A two-dimensional (2D) version of this model is derived and implemented in Section 3. Section 4 validates this implementation using simple shear tests of a constrained layer, for which small and finite deformation results exist in the literature. In Section 5, the 2D version of the proposed SGCP model is applied to key localization problems in single crystals. Finally, the paper ends with our concluding remarks in Section 6.

## 2. Finite Deformation Strain Gradient Crystal Plasticity (SGCP) Model

In this section, the small deformation Gurtin-type SGCP model proposed by Jebahi et al. [34], based on generalized defect energy and uncoupled dissipation, is extended to the finite deformation framework. Following [21,23], the finite deformation model is formulated with respect to the reference configuration to simplify its numerical implementation.

### 2.1. Kinematics and Single Crystal Hypothesis

As in [21,23], within the framework of finite deformation, the elastic-plastic kinematics are described by the Kröner–Lee decomposition of the deformation gradient $F$ [36,37]:

$$F = F^e F^p \qquad (1)$$

where $F^e$ represents the elastic distortion including stretch and rotation of the lattice and $F^p$ represents the plastic distortion due to plastic deformation caused by the flow of defects. The plastic incompressibility is also assumed with $det(F^p) = 1$. Considering (1), the velocity gradient tensor can be written:

$$L = \dot{F}F^{-1} = L^e + F^e L^p F^{e-1} \qquad (2)$$

with the elastic and plastic distortion rate tensors $L^e$ and $L^p$ defined by:

$$L^e = \dot{F}^e F^{-e} \quad \text{and} \quad L^p = \dot{F}^p F^{p-1} \qquad (3)$$

The single crystal hypothesis requires $L^p$ to be governed by slip rates on the individual slip systems. These systems are defined in the intermediate space, which is also called
lattice or structural space. In this space, the slip direction and the slip plane normal direction associated with a slip system $\alpha$, respectively noted $s^\alpha$ and $m^\alpha$, are constant lattice vectors:

$$s^\alpha \cdot m^\alpha = 0 \quad \text{and} \quad |s^\alpha| = |m^\alpha| = 1 \qquad (4)$$

These vectors as well as the glide direction vector ($l^\alpha = m^\alpha \wedge s^\alpha$) are expressed in the reference configuration as [38]:

$$s_r^\alpha = F^{p-1}s^\alpha, \quad m_r^\alpha = F^{pT}m^\alpha, \quad l_r^\alpha = F^{p-1}l^\alpha \qquad (5)$$

Using single crystal hypothesis, the plastic distortion rate tensor $L^p$ can be written within the intermediate space as:

$$L^p = \sum_{\alpha=1}^{q} \dot{\gamma}^{\alpha} S^{\alpha} \qquad (6)$$

where $S^\alpha = s^\alpha \otimes m^\alpha$ denotes the Schmid tensor associated with the slip system $\alpha$, $\dot{\gamma}^\alpha$ represents the plastic slip rate on this system, and $q$ is the number of slip systems. For the sake of consistency with the small deformation SGCP model [34], the same symbol ($\dot{\gamma}^\alpha$) is used to designate the plastic slip rate on the slip system $\alpha$. However, within the finite deformation framework, $\dot{\gamma}^\alpha$ cannot be interpreted as the material time derivative of a physically-meaningful field variable $\gamma^\alpha$. Indeed, slips (as fields on individual slip systems) are not generally well-defined quantities from a physical point of view [23].

## 2.2. Principle of Virtual Power and Governing Equations

Considering the generalized power density of internal forces [39,40], the internal and external virtual power expenditures within an arbitrary subregion $V_0$ of boundary $S_0$ can be expressed in the reference configuration as (neglecting body forces):

$$\mathcal{P}_{\text{int}} = \int_{V_0} \mathbf{P} : \dot{\mathbf{F}} \, d\mathbf{v}_0 + \sum_{\alpha=1}^{q} \int_{V_0} \kappa^\alpha \dot{\gamma}^\alpha \, d\mathbf{v}_0 + \sum_{\alpha=1}^{q} \int_{V_0} \xi^\alpha \cdot \mathbf{grad}_X (\dot{\gamma}^\alpha) \, d\mathbf{v}_0 \qquad (7)$$

$$\mathcal{P}_{\text{ext}} = \int_{S_0} \mathbf{t} \cdot \dot{\mathbf{u}} \, d\mathbf{s}_0 + \sum_{\alpha=1}^{q} \int_{S_0} \chi^\alpha \dot{\gamma}^\alpha \, d\mathbf{s}_0 \qquad (8)$$

where $\mathbf{P} = J\sigma\mathbf{F}^{T-1}$ is the first Piola–Kirchhoff stress tensor (power-conjugate to $\dot{\mathbf{F}}$), $J = \det(\mathbf{F})$ is the volumetric Jacobian, $\sigma$ is the Cauchy stress tensor, $\kappa^\alpha$ is a scalar microscopic stress (power-conjugate to $\dot{\gamma}^\alpha$), $\xi^\alpha$ is a vector microscopic stress (power-conjugate to $\mathbf{grad}_X(\dot{\gamma}^\alpha)$), $\mathbf{t}$ is a macroscopic surface traction, $\chi^\alpha$ is a microscopic surface traction, and $\mathbf{grad}_X(\bullet)$ is the Lagrangian gradient operator (i.e., w.r.t. the Lagrangian coordinates). This operator can be related to the lattice gradient operator $\mathbf{grad}_{\#}(\bullet)$ and the Eulerian gradient operator $\mathbf{grad}_X(\bullet)$ by:

$$\mathbf{grad}_{\#}(\bullet) = \mathbf{grad}_{X}(\bullet) \cdot \mathbf{F}^{p-1} = \mathbf{grad}_{X}(\bullet) \cdot \mathbf{F}^{e} \qquad (9)$$

Introducing the first and second Piola elastic stress tensors, which are, respectively, defined by:

$$\mathbf{P}^{e} = J\sigma\mathbf{F}^{eT-1} \quad \text{and} \quad \mathbf{S}^{e} = J\mathbf{F}^{e-1} \cdot \sigma \cdot \mathbf{F}^{eT-1} \qquad (10)$$

and defining the resolved shear (Schmid) stress $\tau^\alpha$ as:

$$\tau^\alpha = C^e S^e : S^\alpha$$

with $C^e = \mathbf{F}^{eT}\mathbf{F}^e$ is the right Cauchy–Green elastic deformation tensor. The term $\mathbf{P} : \dot{\mathbf{F}}$ in the expression of $\mathcal{P}_{\text{int}}$ (7) can be rewritten as:

$$\mathbf{P} : \dot{\mathbf{F}} = \mathbf{P}^e : \dot{\mathbf{F}}^e + \sum_{\alpha=1}^{q} \tau^\alpha \dot{\gamma}^\alpha \qquad (11)$$
which leads to:

$$P_{int} = \int_{V_0} \mathbf{P}^e : \dot{\mathbf{F}}^e dv_0 + \sum_{\alpha=1}^{q} \int_{V_0} [\pi^\alpha \dot{\gamma}^\alpha + \xi^\alpha \cdot \mathbf{grad}_X(\dot{\gamma}^\alpha)] dv_0$$

with the scalar microscopic stress $\pi^\alpha = \kappa^\alpha + \tau^\alpha$. Application of the generalized principle of virtual power [41,42] results in the following macroscopic and microscopic balance equations and boundary conditions:

$$\begin{gathered} \mathbf{div}_X(\mathbf{P}) = \mathbf{0} \quad \text{in } V_0 \\ \mathbf{P} \cdot \mathbf{n}_0 = \mathbf{t} \quad \text{on } S_0 \end{gathered} \tag{12}$$

$$\begin{gathered} \mathbf{div}_X(\xi^\alpha) + \tau^\alpha - \pi^\alpha = 0 \quad \text{in } V_0 \\ \xi^\alpha \cdot \mathbf{n}_0 = \chi^\alpha \quad \text{on } S_0 \end{gathered} \tag{13}$$

with $\mathbf{n}_0$ the outward unit normal to $S_0$ and $\mathbf{div}_X(\bullet)$ the Lagrangian divergence operator.

## 2.3. Free Energy Imbalance

The free energy imbalance (second thermodynamic law) requires that the temporal increase of free energy in an arbitrary subregion $V_0$ be less than or equal to the external power expended on $V_0$. Mathematically, this can be formulated as:

$$\int_{V_0} \Psi dv_0 - \int_{V_0} \left\{ \mathbf{P}^e : \dot{\mathbf{F}}^e + \sum_{\alpha=1}^{q} (\pi^\alpha \dot{\gamma}^\alpha + \xi^\alpha \cdot \mathbf{grad}_X(\dot{\gamma}^\alpha)) \right\} dv_0 \le 0 \tag{14}$$

where $\Psi$ denotes the free energy measured in the reference configuration. Defining the elastic Green–Lagrange strain tensor $E^e$ as:

$$E^e = \frac{1}{2}(C^e - I) \tag{15}$$

it can easily be shown that $\mathbf{P}^e : \dot{\mathbf{F}}^e = \mathbf{S}^e : \dot{\mathbf{E}}^e$. Since $V_0$ is arbitrary, the local free energy imbalance can be obtained:

$$\Psi - \mathbf{S}^e : \dot{\mathbf{E}}^e - \sum_{\alpha=1}^{q} (\pi^\alpha \dot{\gamma}^\alpha + \xi^\alpha \cdot \mathbf{grad}_X(\dot{\gamma}^\alpha)) \le 0 \tag{16}$$

## 2.4. Constitutive Laws

### 2.4.1. Energetic Constitutive Laws

The free energy in (16) consists of an elastic strain energy $\Psi^e$ and a defect energy $\Psi^\rho$. Building on the work of Gurtin [23], $\Psi^\rho$ is assumed to be a function of dislocation densities:

$$\rho = \{ \rho_{\perp}^1 \cdots \rho_{\perp}^q, \rho_{\odot}^1 \cdots \rho_{\odot}^q \}$$

where $\rho_{\perp}^\alpha$ and $\rho_{\odot}^\alpha$ are, respectively, the edge and screw dislocation densities associated with the slip system $\alpha$, and are, respectively, defined by:

$$\dot{\rho}_{\perp}^{\alpha} = -\mathbf{s}^{\alpha} \cdot \mathbf{grad}_{\#}(\dot{\gamma}^{\alpha}) \quad \text{and} \quad \dot{\rho}_{\odot}^{\alpha} = \mathbf{l}^{\alpha} \cdot \mathbf{grad}_{\#}(\dot{\gamma}^{\alpha}) \tag{17}$$

These densities are in units of length$^{-1}$ and may be positive or negative. Bearing in mind Equations (5) and (9), it can easily be demonstrated that:

$$\dot{\rho}_{\perp}^{\alpha} = -\mathbf{s}_{r}^{\alpha} \cdot \mathbf{grad}_{X}(\dot{\gamma}^{\alpha}) \quad \text{and} \quad \dot{\rho}_{\odot}^{\alpha} = \mathbf{l}_{r}^{\alpha} \cdot \mathbf{grad}_{X}(\dot{\gamma}^{\alpha}) \tag{18}$$
The elastic strain energy $\Psi^e$ is assumed to be a quadratic form of $E^e$:

$$ \Psi^e = \frac{1}{2} E^e : \mathbb{C} : E^e $$

where $\mathbb{C}$ is a symmetric and positive-definite elasticity tensor. In rate form, this quantity can be written as:

$$ \Psi^e = \frac{\partial \Psi^e}{\partial E^e} : \dot{E}^e \qquad (19) $$

As in [34], a generalized power-law defect energy with adjustable order-controlling index $n$ is adopted in this work:

$$ \Psi^\rho = \frac{1}{n} X_0 l_{en}^n \sum_{\alpha=1}^{q} (|\rho_{\perp}^\alpha|^n + |\rho_{\odot}^\alpha|^n) \qquad (20) $$

where $X_0$ is a constant representing the energetic slip resistance, and $l_{en}$ is an energetic length scale. To ensure the convexity of $\Psi^\rho$, the defect energy index $n$ must be greater than or equal to 1 ($n \ge 1$). Its rate form is given by:

$$ \begin{aligned} \dot{\Psi}^\rho &= \sum_{\alpha=1}^{q} \left( \frac{\partial \Psi^\rho}{\partial \rho_{\perp}^\alpha} \dot{\rho}_{\perp}^\alpha + \frac{\partial \Psi^\rho}{\partial \rho_{\odot}^\alpha} \dot{\rho}_{\odot}^\alpha \right) \\ &= X_0 l_{en}^n \sum_{\alpha=1}^{q} \left( -|\rho_{\perp}^\alpha|^{n-2} \rho_{\perp}^\alpha s_r^\alpha + |\rho_{\odot}^\alpha|^{n-2} \rho_{\odot}^\alpha l_r^\alpha \right) \cdot \mathbf{grad}_X(\dot{\gamma}^\alpha) \end{aligned} \qquad (21) $$

Substituting (19) and (21) into (16), the local free energy imbalance can be rewritten as:

$$ \left( \frac{\partial \Psi^e}{\partial E^e} - S^e \right) : \dot{E}^e + \sum_{\alpha=1}^{q} \left[ X_0 l_{en}^n \left( -|\rho_{\perp}^{\alpha}|^{n-2} \rho_{\perp}^{\alpha} s_r^{\alpha} + |\rho_{\odot}^{\alpha}|^{n-2} \rho_{\odot}^{\alpha} l_r^{\alpha} \right) - \zeta^{\alpha} \right] \cdot \mathbf{grad}_X(\dot{\gamma}^{\alpha}) - \sum_{\alpha=1}^{q} \pi^{\alpha} \dot{\gamma}^{\alpha} \le 0 $$

In the following, it is assumed that $S^e$ is fully energetic, $\pi^\alpha$ is fully dissipative, and $\zeta^\alpha$ is divided into energetic and dissipative parts:

$$ S^e = \frac{\partial \Psi^e}{\partial E^e}, \quad \pi^\alpha = \pi_{dis}^\alpha, \quad \text{and} \quad \zeta^\alpha = \zeta_{en}^\alpha + \zeta_{dis}^\alpha $$

It should be noted that higher-order dissipation in the context of strain gradient (crystal) plasticity theories still remains an open question. Although ignored in some works, several authors consider this kind of dissipation in their models [34,43–46]. As reported by Fleck and Willis [45], direct experimental measurement shows that the core energy of dislocations stored during plastic deformation is much smaller than the plastic work dissipated in dislocation motion.

Consequently, both statistically stored and geometrically necessary dislocations would contribute more to plastic dissipation than to a change in free energy. The source of the dissipation process can be attributed to the motion of the mobile dislocations. The motion of dislocations yields energy dissipation, which, in turn, results in resistance to the dislocation motion and increases in the yield strength due to size effects [44].

The energetic constitutive laws can then be obtained:

$$ \begin{aligned} S^e &= \mathbb{C} : E^e \\ \zeta_{en}^\alpha &= X_0 l_{en}^n \left( -|\rho_{\perp}^\alpha|^{n-2} \rho_{\perp}^\alpha s_r^\alpha + |\rho_{\odot}^\alpha|^{n-2} \rho_{\odot}^\alpha l_r^\alpha \right) \end{aligned} \qquad (22) $$

It can be verified from (22)$_2$ that the energetic microscopic stress $\zeta_{en}^\alpha$ is tangent to the $\alpha$th slip plane. Using (22), the dissipation inequality becomes:

$$ \sum_{\alpha=1}^{q} (\pi_{dis}^{\alpha} \dot{\gamma}^{\alpha} + \zeta_{dis}^{\alpha} \cdot \mathbf{grad}_{X}(\dot{\gamma}^{\alpha})) \ge 0 \qquad (23) $$
$\xi_{dis}^{\alpha}$ characterizes the dissipative microscopic forces associated with the evolution of dislocations in the a th slip plane. As the motion of such dislocations is tangent to this plane, $\xi_{dis}^{\alpha}$ is also required to be tangential. Therefore, without losing generality, the gradient term in (23) can be replaced by the corresponding tangential gradient defined by:

$$ \mathbf{grad}_{X}^{\alpha}(\dot{\gamma}^{\alpha}) = (I - \bar{m}_{r}^{\alpha} \otimes \bar{m}_{r}^{\alpha}) \cdot \mathbf{grad}_{X}(\dot{\gamma}^{\alpha}), \quad \text{with} \quad \bar{m}_{r}^{\alpha} = \frac{m_{r}^{\alpha}}{|m_{r}^{\alpha}|} $$

This yields:

$$ \sum_{\alpha=1}^{q} (\pi_{dis}^{\alpha} \dot{\gamma}^{\alpha} + \xi_{dis}^{\alpha} \cdot \mathbf{grad}_{X}^{\alpha}(\dot{\gamma}^{\alpha})) \geq 0 \quad (24) $$

### 2.4.2. Dissipative Constitutive Laws

Referring to [34], uncoupled effective plastic strain measures are used to describe the first- and higher-order dissipative effects. These measures, which are referred to hereafter as first- and higher-order measures, are, respectively, defined by:

$$ \begin{gathered} \dot{e}_{\pi}^{\alpha} = |\dot{\gamma}^{\alpha}|, \quad e_{\pi}^{\alpha} = \int \dot{e}_{\pi}^{\alpha} dt \\ \dot{e}_{\zeta}^{\alpha} = l_{dis} \|\mathbf{grad}_{X}^{\alpha}(\dot{\gamma}^{\alpha})\|, \quad e_{\zeta}^{\alpha} = \int \dot{e}_{\zeta}^{\alpha} dt \end{gathered} $$

where $l_{dis}$ is dissipative internal length scale. To ensure positive dissipation, the scalar and vector dissipative stresses associated with slip system $\alpha$ are taken as:

$$ \pi_{dis}^{\alpha} = S_{\pi}^{\alpha} \left( \frac{\dot{e}_{\pi}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\dot{\gamma}^{\alpha}}{e_{\pi}^{\alpha}} \quad (25) $$

$$ \xi_{dis}^{\alpha} = S_{\zeta}^{\alpha} l_{dis}^{2} \left( \frac{\dot{e}_{\zeta}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\mathbf{grad}_{X}^{\alpha}(\dot{\gamma}^{\alpha})}{e_{\zeta}^{\alpha}} \quad (26) $$

where $\dot{\gamma}_{0}^{\alpha} > 0$ is a constant strain rate representative of the flow rates of interest, $m > 0$ is a constant characterizing the rate-sensitivity of the considered material, and $S_{i}^{\alpha}$ ($i \in \{\pi, \zeta\}$) are stress-dimensional internal state variables ($S_{\pi}^{\alpha} > 0$ and $S_{\zeta}^{\alpha} \geq 0$). These variables are called, hereafter, dissipative slip resistances. Their evolutions are assumed to be governed by:

$$ \dot{S}_{i}^{\alpha} = \sum_{\beta=1}^{q} h_{i}^{\alpha\beta} (S_{i}^{\beta}) \dot{e}_{i}^{\beta} \quad \text{with} \quad S_{i}^{\alpha}(0) = S_{i0} \quad (27) $$

$S_{i0}$ ($i \in \{\pi, \zeta\}$) are initial dissipative slip resistances ($S_{\pi0} > 0$ and $S_{\zeta0} \geq 0$), and $h_{i}^{\alpha\beta}$ ($i \in \{\pi, \zeta\}$) are positive hardening moduli.

### 2.5. Flow Rule

Based on the expressions of the microscopic stresses, the microscopic balance equation associated with slip system $\alpha$ can be reformulated as follows:

$$ \tau^{\alpha} + \mathbf{div}_{X}(\xi_{en}^{\alpha}) + \mathbf{div}_{X}(\xi_{dis}^{\alpha}) - \pi_{dis}^{\alpha} = 0 \quad (28) $$

Substituting (22), (25) and (26) into (28), the flow rule associated with slip system $\alpha$ can be obtained:

$$ \tau^{\alpha} - X_{0} l_{en}^{n} \mathbf{div}_{X} \left( |\rho_{\perp}^{\alpha}|^{n-2} \rho_{\perp}^{\alpha} s_{r}^{\alpha} - |\rho_{\odot}^{\alpha}|^{n-2} \rho_{\odot}^{\alpha} l_{r}^{\alpha} \right) = \underbrace{S_{\pi}^{\alpha} \left( \frac{\dot{e}_{\pi}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\dot{\gamma}^{\alpha}}{e_{\pi}^{\alpha}}}_{\text{II}} - \underbrace{\mathbf{div}_{X} \left( S_{\zeta}^{\alpha} l_{dis}^{2} \left( \frac{\dot{e}_{\zeta}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\mathbf{grad}_{X}^{\alpha}(\dot{\gamma}^{\alpha})}{e_{\zeta}^{\alpha}} \right)}_{\text{III}} \quad (29) $$

Term I in the above equation, being energetic represents a backstress, which leads to kinematic hardening (Bauschinger) effects. Terms II and III are dissipative and lead
to dissipative hardening effects (a reversal of the flow direction, i.e., $\dot{\gamma}^{\alpha} \to -\dot{\gamma}^{\alpha}$, simply changes the sign of these terms).

## 2.6. Microscopic Boundary Conditions

Within the framework of higher-order gradient theories, the presence of higher-order stresses can result in an expenditure of microscopic power on the boundary of the studied domain ($S_0$). For each slip system $\alpha$, the microscopic power expended, per unit area, on this boundary is given by:

$$ (\xi^{\alpha} \cdot n_0) \dot{\gamma}^{\alpha} \qquad (30) $$

The present work focuses on simple microscopic boundary conditions (BCs) that result in a null power expenditure on $S_0$:

$$ (\xi^{\alpha} \cdot n_0) \dot{\gamma}^{\alpha} = 0 \quad \text{on } S_0 \qquad (31) $$

Specifically, two limit cases of microscopic BCs are considered in this work:

$$ \begin{aligned} \dot{\gamma}^{\alpha} &= 0 \quad \text{on } S_h \\ \xi^{\alpha} \cdot n_0 &= 0 \quad \text{on } S_f \end{aligned} \qquad (32) $$

with $S_h$ and $S_f$ complementary subsurfaces of $S_0$. The micro-hard BCs (32)1 correspond to impenetrable surfaces with no slip, whereas the micro-free BCs (32)2 represent free surfaces without slip resistance.

# 3. Numerical Implementation

In this section, a two-dimensional version of the finite deformation SGCP model is implemented under plane strain conditions with planar slip systems, using a total Lagrangian finite element formulation. An extension to the generalized three-dimensional case can easily be derived.

## 3.1. Two-Dimensional Version of the Proposed Model

Assuming plane strain conditions with planar slip systems, it can easily be demonstrated that screw dislocation densities vanish [34,43,47]:

$$ \dot{\rho}_{\odot}^{\alpha} = l^{\alpha} \cdot \mathbf{grad}_{\#}(\dot{\gamma}^{\alpha}) = 0 \qquad (33) $$

Therefore, the energetic microscopic stresses become:

$$ \xi_{en}^{\alpha} = -X_0 l_{en}^{n} (|\rho_{\perp}^{\alpha}|^{n-2} \rho_{\perp}^{\alpha} s_{r}^{\alpha}). $$

The dissipative microscopic stresses are given by (25) and (26), with the slip resistances $S_{\pi}^{\alpha}$ and $S_{\xi}^{\alpha}$ assumed, hereafter, to evolve linearly with the effective plastic strain measures:

$$ \dot{S}_{i}^{\alpha} = \sum_{\beta=1}^{q} H_{i} \dot{e}_{i}^{\beta} \quad \text{with} \quad S_{i}^{\alpha}(0) = S_{i0}, \quad H_{i} = \text{constant} \geq 0 \qquad (34) $$

Note that $i \in \{\pi, \xi\}$, $S_{\pi0} > 0$ and $S_{\xi0} \geq 0$. Using the above expressions of the microscopic stresses, a simplified form of the flow rule (29) can be obtained:

$$ \tau^{\alpha} - X_0 l_{en}^{n} \mathbf{div}_X \left( |\rho_{\perp}^{\alpha}|^{n-2} \rho_{\perp}^{\alpha} s_r^{\alpha} \right) = S_{\pi}^{\alpha} \left( \frac{\dot{e}_{\pi}^{\alpha}}{\dot{\gamma}_0^{\alpha}} \right)^{m} \frac{\dot{\gamma}_0^{\alpha}}{\dot{e}_{\pi}^{\alpha}} - \mathbf{div}_X \left( S_{\xi}^{\alpha} l_{dis}^{2} \left( \frac{\dot{e}_{\xi}^{\alpha}}{\dot{\gamma}_0^{\alpha}} \right)^{m} \mathbf{grad}_{X}^{\alpha}(\dot{\gamma}_0^{\alpha}) \right) \qquad (35) $$

In summary, key equations characterizing the simplified two-dimensional SGCP model are recalled below:
$$
\left\{
\begin{aligned}
& F = F^e F^p, \quad L^p = \sum_{\alpha=1}^{q} \dot{\gamma}^{\alpha} S^{\alpha}, \quad S^{\alpha} = s^{\alpha} \otimes m^{\alpha} \\
& \dot{F}^p = L^p F^p, \quad F^e = FF^{p-1}, \quad C^e = F^{eT} F^e, \quad E^e = \frac{1}{2}(C^e - I) \\
& S^e = \mathbb{C} : E^e, \quad P = F^e S^e F^{pT-1}, \quad \sigma = J^{-1} F^e S^e F^{eT} \\
& \tau^{\alpha} + \mathbf{div}_X(\xi^{\alpha}) - \pi_{dis}^{\alpha} = 0 \\
& \pi_{dis}^{\alpha} = S_{\pi}^{\alpha} \left( \frac{\dot{e}_{\pi}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\dot{\gamma}^{\alpha}}{\dot{e}_{\pi}^{\alpha}} \\
& \xi^{\alpha} = X_0 l_{en}^n \left( -|\rho_{\perp}^{\alpha}|^{n-2} \rho_{\parallel}^{\alpha} s_r^{\alpha} \right) + S_{\xi}^{\alpha} l_{dis}^2 \left( \frac{\dot{e}_{\xi}^{\alpha}}{\dot{\gamma}_{0}^{\alpha}} \right)^{m} \frac{\mathbf{grad}_X^{\alpha}(\dot{\gamma}^{\alpha})}{\dot{e}_{\xi}^{\alpha}} \\
& \dot{S}_{\pi}^{\alpha} = \sum_{\beta=1}^{q} H_{\pi} (S_{\pi}^{\beta}) \dot{e}_{\pi}^{\beta}, \quad S_{\pi 0}^{\alpha} > 0 \\
& \dot{S}_{\xi}^{\alpha} = \sum_{\beta=1}^{q} H_{\xi} (S_{\xi}^{\beta}) \dot{e}_{\xi}^{\beta}, \quad S_{\xi 0}^{\alpha} \geq 0 \\
& \rho_{\parallel}^{\alpha} = -s_r^{\alpha} \cdot \mathbf{grad}_X(\dot{\gamma}^{\alpha}) \\
& \dot{e}_{\pi}^{\alpha} = |\dot{\gamma}^{\alpha}|, \quad \dot{e}_{\xi}^{\alpha} = l_{dis} \| \mathbf{grad}_X(\dot{\gamma}^{\alpha}) \|
\end{aligned}
\right.
\quad (36)
$$

These equations will be solved numerically within a total Lagrangian finite element framework. Several methods are used in the literature to solve these equations. Among them, one can cite the approach proposed by Scherer et al. [48] combining the advantages of micromorphic and Lagrange multiplier methods. In this approach, a scalar micromorphic variable is linked to a cumulative plastic slip measure by a Lagrange multiplier to ensure uncoupling between constitutive nonlinearly and spatial nonlocality.

Other methods based on an independent thermodynamic variable for each slip system can be found, allowing for more flexible control of gliding on the slip systems. Two common examples of these methods can be cited: the ρ-method [11,49–52] and the γ-method [34,43,53]. The former considers dislocation densities or their rates as primary variables (additional degrees of freedom), whereas plastic slips or plastic slip rates are considered as primary variables in the latter.

According to Klusemann et al. [50], the two methods generally lead to similar results. However, γ-methods provide more flexibility in constraining the dislocation flow on the boundaries. Furthermore, it does not require a local loop (at the integration points) to calculate the plastic slip rates γ^a, which are given by solving the global system of equations. This method is chosen in the present work. In addition to the displacement field u, the scalar field variables defined on each slip system α by:

$$
\gamma^{\alpha} = \int_{0}^{t} \dot{\gamma}^{\alpha} dt, \quad \text{with} \quad \gamma^{\alpha}(t=0) = 0 \qquad (37)
$$

are taken as primary variables, in consistency with the implemented small deformation version of the proposed SGCP model [34]. It is worth recalling that these variables, which can be interpreted in small deformation as plastic slips on the slip systems, have no real physical interpretation in finite deformation (only their increments are physically relevant). As will be seen in the validation section, this choice of primary variables has no impact on the numerical results, as the integration of the constitutive laws is only based on the increments of these variables.

3.2. *Integration of Constitutive Laws*

Numerical integration of constitutive laws consists in updating the values of stresses and internal variables, given the increments of kinematic variables. Knowing the increments of the primary field variables (Δu and Δγ^a), the displacement, the deformation gradient, and the plastic distortion rate quantities can be calculated at a given point as:
$$ \mathbf{u}_{n+1} = \mathbf{u}_n + \Delta\mathbf{u}, \quad \mathbf{F}_{n+1} = \mathbf{grad}_X(\mathbf{u}_{n+1}) + \mathbf{I}, \quad \text{and} \quad \mathbf{L}_{n+1}^p = \sum_{\alpha=1}^q \left( \frac{\Delta\gamma^\alpha}{\Delta t} \mathbf{S}^\alpha \right) \quad (38) $$

with $\mathbf{I}$ as the second order unity tensor. The plastic distortion tensor can be determined based on Equation (3)$_2$. Using the implicit approach, this equation can be written in discretized form as:

$$ \mathbf{F}_{n+1}^p - \mathbf{F}_n^p = \left[ \sum_{\alpha=1}^q \Delta\gamma^\alpha \mathbf{S}^\alpha \right] \mathbf{F}_{n+1}^p \quad (39) $$

which leads to:

$$ \mathbf{F}_{n+1}^{p-1} = \mathbf{F}_n^{p-1} \left[ \mathbf{I} - \sum_{\alpha=1}^{q} \Delta\gamma^\alpha \mathbf{S}^\alpha \right] \quad (40) $$

To ensure plastic incompressibility (det $\mathbf{F}^p = 1$), a unimodular operator is applied to the term between brackets:

$$ \mathbf{F}_{n+1}^{p-1} = \mathbf{F}_n^{p-1} \text{uni} \left( \mathbf{I} - \sum_{\alpha=1}^{q} \Delta\gamma^\alpha \mathbf{S}^\alpha \right) \quad (41) $$

where, for a second order tensor $\mathbf{A}$, the term uni($\mathbf{A}$) is defined by:

$$ \text{uni}(\mathbf{A}) = [\det(\mathbf{A})]^{-\frac{1}{3}} \mathbf{A} \quad (42) $$

Considering (1), the elastic part of deformation gradient is calculated by:

$$ \mathbf{F}_{n+1}^e = \mathbf{F}_{n+1} \mathbf{F}_{n+1}^{p-1} \quad (43) $$

Using the above equation, the elastic Green–Lagrange strain tensor $\mathbf{E}_{n+1}^e$ can be obtained:

$$ \mathbf{E}_{n+1}^e = \frac{1}{2} (\mathbf{C}_{n+1}^e - \mathbf{I}), \quad \text{with} \quad \mathbf{C}_{n+1}^e = \mathbf{F}_{n+1}^{eT} \mathbf{F}_{n+1}^e \quad (44) $$

Knowing the updated kinematic quantities, the stresses and internal variables can be updated. Algorithm 1 summarizes the steps of integration of constitutive laws, where the subscripts “n” and “n + 1” are omitted for clarity.

### 3.3. Weak Forms and Finite Element Discretization

Assume the virtual fields $\delta\dot{\mathbf{u}}$ and $\delta\dot{\gamma}^\alpha$, which are kinematically admissible to zero on the portions of the boundary of the studied domain on which macroscopic and microscopic Dirichlet boundary conditions are respectively imposed. The weak forms of the macroscopic and microscopic balance equations (first lines of (12) and (13)) can be obtained as follows:

$$ \begin{aligned} G_u &= \int_{V_0} \mathbf{grad}_X(\delta\dot{\mathbf{u}}) : \mathbf{P} \, dv - \int_{S_0^t} \delta\dot{\mathbf{u}} \cdot \mathbf{t} \, ds \\ G_\gamma &= \sum_{\alpha=1}^q \left( \int_{V_0} \mathbf{grad}_X(\delta\dot{\gamma}^\alpha) \cdot \mathbf{\zeta}^\alpha \, dv + \int_{V_0} \delta\dot{\gamma}^\alpha \pi^\alpha \, dv - \int_{V_0} \delta\dot{\gamma}^\alpha \tau^\alpha \, dv \right) - \sum_{\alpha=1}^q \int_{S_0^{\chi^\alpha}} \delta\dot{\gamma}^\alpha \chi^\alpha \, ds \end{aligned} \quad (45) $$

where $S_0^t$ and $S_0^{\chi^\alpha}$ are the portions of the domain boundary on which macroscopic and microscopic traction forces, noted, respectively, as $\mathbf{t}$ and $\chi^\alpha$, are imposed. These weak forms are solved numerically using the total Lagrangian finite element method. To this end, a User-Element (UEL) subroutine is implemented within the commercial finite element package ABAQUS/Standard. In this UEL, displacement and plastic scalar fields ($\mathbf{u}$ and $\gamma^\alpha$) are considered as primary variables (degrees of freedom).

Several finite element types can be used to approximate these fields. In this work, an element type employing eight-noded (quadratic) interpolation of $\mathbf{u}$, four-noded (linear) interpolation of $\gamma^\alpha$ and $2 \times 2$ Gauss integration is adopted (Figure 1). Based on Kuroda [11], this element type is well suited to study localization problems. Within each element, dis-
placements $u_i$ ($i = 1, 2$) and plastic scalar variables $\gamma^\alpha$ ($\alpha = 1, 2, \dots, q$) can be approximated based on the associated nodal values as follows:

$$u^i = \sum_{k=1}^{8} N_u^k U_k^i, \quad \gamma^\alpha = \sum_{k=1}^{4} N_\gamma^k \Gamma_k^\alpha \qquad (46)$$

where $U_k^i$ and $\Gamma_k^\alpha$ are nodal values of displacements and plastic scalar variables, and $N_u^k$ and $N_\gamma^k$ are the associated interpolation functions.

**Algorithm 1 Integration of constitutive equations.**

**Inputs:** $u, F^p, \rho_{\perp}^\alpha, \Delta u, \Delta \gamma^\alpha$,

**Update:** $u \leftarrow u + \Delta u$

**Calculate:** $F, F^{-p}, F^e, C^e$ and $E^e$

**Calculate:** $J = \det(F)$

**Calculate:** $S^e = \mathbb{C} : E^e, P = F^e \cdot S^e \cdot F^{pT-1}, \sigma = J^{-1} F^e \cdot S^e \cdot F^{eT}$

**Do $\alpha = 1$ to $q$:**

**Calculate:** $\Delta e_{\pi}^{\alpha} = |\Delta \gamma^{\alpha}|, \quad \Delta e_{\zeta}^{\alpha} = l_{dis} \|\text{grad}_{X}^{\alpha}(\Delta \gamma^{\alpha})\|$

**end do**

**Do $\alpha = 1$ to $q$:**

**Calculate:** $\Delta \rho_{\perp}^{\alpha} = -s_r^{\alpha} \cdot \text{grad}_{X}^{\alpha}(\Delta \gamma^{\alpha}), \Delta S_{\pi}^{\alpha} = \sum_{\beta=1}^{q} H_{\pi} \Delta e_{\pi}^{\beta}, \Delta S_{\zeta}^{\alpha} = \sum_{\beta=1}^{q} H_{\zeta} \Delta e_{\zeta}^{\beta}$

**Update:** $\rho_{\perp}^{\alpha} \leftarrow \rho_{\perp}^{\alpha} + \Delta \rho_{\perp}^{\alpha}, S_{\pi}^{\alpha} \leftarrow S_{\pi}^{\alpha} + \Delta S_{\pi}^{\alpha}, S_{\zeta}^{\alpha} \leftarrow S_{\zeta}^{\alpha} + \Delta S_{\zeta}^{\alpha}$

**Calculate:** $\pi_{dis}^{\alpha} = S_{\pi}^{\alpha} \left( \frac{\Delta e_{\pi}^{\alpha}}{\dot{\gamma}_{0}^{\alpha} \Delta t} \right)^{m} \frac{\Delta \gamma^{\alpha}}{\Delta e_{\pi}^{\alpha}}$

**Calculate:** $\zeta_{n+1}^{\alpha} = X_0 l_{en}^n \left( -|\rho_{\perp}^{\alpha}|^{n-2} \rho_{\perp}^{\alpha} s_r^{\alpha} \right) + S_{\zeta}^{\alpha} l_{dis}^2 \left( \frac{\Delta e_{\zeta}^{\alpha}}{\dot{\gamma}_{0}^{\alpha} \Delta t} \right)^{m} \frac{\text{grad}_{X}^{\alpha}(\Delta \gamma^{\alpha})}{\Delta e_{\zeta}^{\alpha}}$

**end do**

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The finite element used in the present work.</figcaption></figure>

Using the above field approximations, the macroscopic and microscopic weak forms can be written within a representative finite element, in a simplified matrix form, as follows:

$$G_u^e = \left( \delta \dot{U}^e \right)^T \cdot \left( \int_{V_{e0}} \underline{B}_u^T \cdot \underline{P} \, dv - \int_{V_{e0}} \underline{N}_u^T \cdot \underline{t} \, ds \right) \\ G_\gamma^e = \left( \delta \dot{\underline{\Gamma}}^e \right)^T \cdot \left( \int_{V_{e0}} \underline{B}_\gamma^T \cdot \underline{\zeta} \, dv + \int_{V_{e0}} \underline{N}_\gamma^T \cdot \underline{\pi} \, dv - \int_{V_{e0}} \underline{N}_\gamma^T \cdot \underline{\tau} \, dv - \int_{S_{e0}^\chi} \underline{N}_\gamma^T \cdot \underline{\chi} \, ds \right) \qquad (47)$$

where $\underline{N}_u$ and $\underline{B}_u$ are the interpolation and gradient matrices associated with the displacement field, $\underline{N}_\gamma$ and $\underline{B}_\gamma$ are the interpolation and gradient matrices associated with the plastic scalar fields, $\underline{P}, \underline{\zeta}, \underline{\pi}$ are vector representations of the macroscopic stress and microscopic stresses on all slip systems, $\underline{\tau}$ is vector of resolved shear stresses on all slip systems, $\underline{t}$ and $\underline{\chi}$ are vectors of macroscopic and microscopic traction forces.
Considering the principle of virtual power ($G_u^e = 0$ and $G_\gamma^e = 0$, for any arbitrary $\delta\dot{U}^e$ and $\delta\dot{\Gamma}^e$), the associated residuals can be written:

$$ \underline{R}_u^e = \int_{V_{e0}} \underline{B}_u^T \cdot \underline{P} \, dv - \int_{S_{e0}^t} \underline{N}_u^T \cdot \underline{t} \, ds = 0 \quad (48) $$

$$ \underline{R}_\gamma^e = \int_{V_{e0}} \underline{B}_\gamma^T \cdot \underline{\zeta} \, dv + \int_{V_{e0}} \underline{N}_\gamma^T \cdot \underline{\pi} \, dv - \int_{V_{e0}} \underline{N}_\gamma^T \cdot \underline{\tau} \, dv - \int_{S_{e0}^\chi} \underline{N}_\gamma^T \cdot \underline{\chi} \, ds = 0 $$

Linearization of these equations with respect to small variations of $\underline{U}^e$ and $\underline{\Gamma}^e$ leads to an elementary system of equations, which can be written in matrix form as:

$$ \left[ \begin{array}{cc} \underline{K}_{uu}^e & \underline{K}_{u\gamma}^e \\ \underline{K}_{\gamma u}^e & \underline{K}_{\gamma\gamma}^e \end{array} \right] \left[ \begin{array}{c} \Delta \underline{U}^e \\ \Delta \underline{\Gamma}^e \end{array} \right] = - \left[ \begin{array}{c} \underline{R}_u^e \\ \underline{R}_\gamma^e \end{array} \right] \quad (49) $$

with

$$ \underline{K}_{uu} = \frac{\partial \underline{R}_u^e}{\partial \Delta \underline{U}^e}, \quad \underline{K}_{u\gamma} = \frac{\partial \underline{R}_u^e}{\partial \Delta \underline{\Gamma}^e}, \quad \underline{K}_{\gamma u} = \frac{\partial \underline{R}_\gamma^e}{\partial \Delta \underline{U}^e}, \quad \text{and} \quad \underline{K}_{\gamma\gamma} = \frac{\partial \underline{R}_\gamma^e}{\partial \Delta \underline{\Gamma}^e} \quad (50) $$

The global system of linear equations can be obtained by assembling all the elementary systems of equations associated with the overall finite elements. This system is solved by means of a Newton-Raphson iterative solution scheme. At each iteration, updated values of the overall increments of the displacement and plastic scalar fields ($\Delta U$ and $\Delta \Gamma$, respectively) are obtained and used to numerically solve the constitutive equations at the Gauss points (Algorithm 1).

### 4. Validation of the Model Implementation

To ensure that the finite deformation SGCP model is correctly implemented, it is applied to simulate a simple shear problem of a constrained layer under plane strain conditions. The considered layer is first subjected to small shear loading, allowing for comparison with the small deformation version of the proposed model [34]. This makes it possible to draw conclusions about the effectiveness of the model implementation to correctly reproduce different energetic and dissipative size effects, for which large deformation results are not necessarily available in the literature.

Then, large shear loading is applied for a particular case of quadratic defect energy ($n = 2$) with no higher-order dissipation ($l_{dis} = 0$). For this case, the present model is equivalent to the model of Kuroda and Tvergaard [49], for which large deformation shear results are available. These results are used to assess the ability of the model implementation to correctly capture finite deformation effects.

Figure 2 presents the geometrical model used in this section. A strip of infinite width in the $x_1$-direction and height $H$ in the $x_2$-direction is subjected to simple shear loading under plane strain conditions. Two slip systems symmetrically oriented with respect to the $x_1$-direction ($\theta_1 = 60^\circ$ and $\theta_2 = -60^\circ$) are considered. The macroscopic boundary conditions are given by:

$$ \left. \begin{array}{l} u_1 = 0, \quad u_2 = 0 \quad \text{along} \quad x_2 = 0 \\ u_1 = H\Gamma, \quad u_2 = 0 \quad \text{along} \quad x_2 = H \end{array} \right\} \quad (51) $$

where $\Gamma$ is the prescribed shear strain. Micro-hard boundary conditions are imposed on the top and bottom edges to hinder dislocation flow at them (impenetrable to dislocations):

$$ \dot{\gamma}^1 = \dot{\gamma}^2 = 0 \quad \text{along} \quad x_2 = 0, H $$

To model the infinite width of the strip, periodic boundary conditions are imposed at its left and right edges:

$$ \left. \begin{array}{l} u_i(0, x_2, t) = u_i(W, x_2, t) \quad \text{for} \quad i = 1, 2 \\ \gamma^\alpha(0, x_2, t) = \gamma^\alpha(W, x_2, t) \quad \text{for} \quad \alpha = 1, 2 \end{array} \right\} $$
With these conditions, the strip width *W* does not affect the results and can arbitrarily be chosen. To save computation time, *W* is chosen to have only one finite element along the *x*₁-direction. A numerical sample consisting of a column of 100 finite elements in the *x*₂-direction is then adopted.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Simple shear of a constrained strip with two active slip systems.</figcaption></figure>

For the case of small shear loading, a complete loading cycle is imposed (i.e., $0 \rightarrow \Gamma \rightarrow -\Gamma \rightarrow \Gamma$), with the prescribed shear strain restricted to $\Gamma = 0.02$. The material parameters are chosen as in [34]. These parameters are presented in Table 1. Figure 3 compares the associated results obtained using the present finite deformation model with those obtained using the small deformation one [34], for different energetic and dissipative parameters. Solid line curves correspond to the finite deformation results and those with markers correspond to the small deformation ones.

For comparison purposes, results obtained using conventional small and finite deformation crystal plasticity with no hardening ($H_\pi = 0$) are also presented in Figure 3 (curves named “Perfect CP”). It can be seen that the stress–strain responses obtained using the finite deformation formulation show very good agreement with those obtained using the small deformation one. Note that the stress–strain curves obtained using $n = 1.7$ present inflection points, reflecting nonlinear kinematic hardening type III (KIII) of Asaro [54]. Further discussions about this hardening type can be found in [34,55,56]. Applying small shear loading, the nonlinear geometrical effects are insignificant. In this case, the small and finite deformation formulations must lead to similar results, which is the case according to Figure 3. This validates the capabilities of the implemented finite deformation model in correctly reproducing the energetic and dissipative effects, at least in the range of the prescribed shear loading.

<table><caption>**Table 1.** The material parameters implied for simple shear of constrained layer undergoing small deformation, taken from [34] for comparison purposes.</caption><thead><tr><th>Material Parameter Name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i></td><td>260</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i></td><td>0.3</td><td>—</td></tr><tr><td>Reference slip rate</td><td><i>γ</i><sub>0</sub><sup>α</sup></td><td>0.04</td><td>s<sup>-1</sup></td></tr><tr><td>Rate-sensitivity parameter</td><td><i>m</i></td><td>0.05</td><td>—</td></tr><tr><td>Higher-order energetic slip resistance</td><td><i>X</i><sub>0</sub></td><td>50</td><td>MPa</td></tr><tr><td>Initial first-order dissipative slip resistance</td><td><i>S</i><sub>π0</sub></td><td>50</td><td>MPa</td></tr><tr><td>Initial higher-order dissipative slip resistance</td><td><i>S</i><sub>ξ0</sub></td><td>50</td><td>MPa</td></tr><tr><td>First-order hardening modulus</td><td><i>H</i><sub>π</sub></td><td>0</td><td>MPa</td></tr><tr><td>Higher-order hardening modulus</td><td><i>H</i><sub>ξ</sub></td><td>0</td><td>MPa</td></tr><tr><td>Energetic length scale</td><td><i>l</i><sub>en</sub></td><td>case study</td><td>μm</td></tr><tr><td>Dissipative length scale</td><td><i>l</i><sub>dis</sub></td><td>case study</td><td>μm</td></tr><tr><td>Defect energy index</td><td><i>n</i></td><td>case study</td><td>—</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Small shear loading results for different energetic and dissipative parameters: comparison between the small and finite deformation models (small deformation results are taken from [34]).</figcaption></figure>

To examine the effectiveness of the present model implementation with regard to nonlinear geometrical effects, the simple shear example of [49], which implies large shear loading, is reproduced in this work. In this example, the same geometrical sample as in Figure 2 is used, with imposed shear loading $\Gamma = 0.24$. The defect energy index and the dissipative length scale are set as $n = 2$ and $l_{dis} = 0$, respectively. With these choices, the present finite deformation model is equivalent to that of [49], in which only energetic back stress is considered.

The other material parameters, which are taken from [49], are presented in Table 2. Figure 4 compares the results of the present model with those obtained by [49] for two values of energetic length scale $l_{en}$. Very good agreement between the results is obtained, regardless of the value of $l_{en}$. This shows the ability of the present model implementation to correctly reproduce finite deformation effects.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Large shear loading results in terms of the (a) stress–strain curves and (b) distributions of the total plastic variables $\gamma^\alpha$ along the $x_2$ direction: comparison between the present model and the model of [49] for different values of $l_{en}/H$.</figcaption></figure>

As can be seen from Figure 4a, in the range of small shear strains, until approximately $\Gamma = 0.025$, nonlinear geometrical effects are insignificant and the shear stress $\sigma_{12}$ increases almost linearly with $\Gamma$ after yielding. Beyond approximately $\Gamma = 0.025$, nonlinear stress–strain dependence takes place and becomes more striking at larger strains, due to nonlinear geometrical effects. Figure 4b presents the distributions of the total plastic variables $\gamma^\alpha$ associated with the considered slip systems along $x_2$-direction, as obtained using the present and the [49] models.
As shown, the nonlinear geometrical effects lead to asymmetric evolution of these variables, despite the initial symmetry of the slip systems with respect to the loading direction. A distinct deviation between the $\gamma^1$ and $\gamma^2$ profiles, which coincide with those obtained by [49], is observed. Such a deviation was also reported by Bittencourt [32] using a similar simple shear problem under monotonic large shear loading. This can be explained by an apparent hardening increase in one system and an apparent hardening decrease in the other system due to lattice rotation effects [32,49].

<table><caption>**Table 2.** Material parameters implied for the simple shear of a constrained layer undergoing finite deformation, inspired by [49].</caption><thead><tr><th>Material Parameter Name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i></td><td>130</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i></td><td>0.3</td><td>—</td></tr><tr><td>Reference slip rate</td><td>γ̇<sub>0</sub><sup>α</sup></td><td>0.001</td><td>s<sup>-1</sup></td></tr><tr><td>Rate-sensitivity parameter</td><td><i>m</i></td><td>0.02</td><td>—</td></tr><tr><td>Higher-order energetic slip resistance</td><td><i>X</i><sub>0</sub></td><td>50</td><td>MPa</td></tr><tr><td>Initial first-order dissipative slip resistance</td><td><i>S</i><sub>π0</sub></td><td>50</td><td>MPa</td></tr><tr><td>Initial higher-order dissipative slip resistance</td><td><i>S</i><sub>ξ0</sub></td><td>50</td><td>MPa</td></tr><tr><td>First-order hardening modulus</td><td><i>H</i><sub>π</sub></td><td>250</td><td>MPa</td></tr><tr><td>Higher-order hardening modulus</td><td><i>H</i><sub>ξ</sub></td><td>0</td><td>MPa</td></tr><tr><td>Energetic length scale</td><td><i>l</i><sub>en</sub></td><td>case study</td><td>μm</td></tr><tr><td>Dissipative length scale</td><td><i>l</i><sub>dis</sub></td><td>0</td><td>μm</td></tr><tr><td>Defect energy index</td><td><i>n</i></td><td>2</td><td>—</td></tr></tbody></table>

## 5. Application to Strain Localization in Single Crystals

Several works studying strain localization in single crystals undergoing single and multiple slips can be found in the literature [3,4,6,8–10]. Using multiple slip systems, it has been shown that strain localization may occur along narrow shear bands, which can be arbitrarily oriented with respect to the initial crystallographic directions [3,10,11]. In the case of a single slip system, two types of shear bands, known as slip and kink bands, can take place [4]. Slip bands correspond to slip localization along a crystallographic slip plane, whereas kink bands are orthogonal to the slip direction (Figure 5).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of slip and kink bands within a single crystal having only one slip system.</figcaption></figure>

Due to lack of internal length scale(s), conventional crystal plasticity theories are not able to correctly capture localization bands. As an alternative, Forest [6] has applied a small deformation Cosserat single crystal plasticity theory to investigate the formation of
slip and kink bands within single crystals. As shown by the author, this theory allows for regularizing kink bands that will have a finite thickness. However, it does not regularize slip bands, which are not necessarily associated with the lattice curvature [6,57].

Although not well investigated in the literature, SGCP theories, which are classically based on GND densities, are also expected to present difficulties in regularizing idealized slip bands. Indeed, a slip gradient along the normal to the slip planes does not result in the storage of GNDs. Ling et al. [29] proposed a finite deformation micromorphic crystal plasticity model that is capable of regularizing both kink and slip bands.

Application of this model to strain localization in single crystals presenting one slip system has brought to light a complex localization behavior with competition between kink and slip bands. Recently, Marano et al. [9] applied a small deformation version of the SGCP model of Gurtin [21], which ignores higher-order dissipation to study localization phenomena within poly-crystal solids. Interesting discussion on the capabilities of the applied model in capturing the competition between slip and kink bands was provided in this study, using fully higher-order energetic effects.

In almost all existing studies of localization phenomena in single crystals, only higher-order energetic effects were considered. Higher-order dissipative effects on these phenomena, particularly the competition between slip and kink bands, have not yet been explored. The objective of the present section is to provide an in-depth investigation of the capabilities and limitations of the implemented finite deformation SGCP model in predicting different kinds of localization modes within single crystals, considering both higher-order energetic and dissipative effects. Particular focus will be given to the ability of such a model to reproduce the complex localization behavior encountered considering a single slip system, with competition between slip and kink bands.

## 5.1. Regularization of Localization Bands by SGCP Theories

Before applying the implemented SGCP model to study complex localization phenomena, the present subsection aims at reviewing its ability to capture different kinds of localization bands, including idealized slip and kink bands as well as general shear bands (i.e., with no particular orientation with respect to the initial crystallographic directions).

### 5.1.1. Simple Shear with Single Slip System: Regularization of Idealized Slip/Kink Bands

Simple shear simulations using single crystals with single slip are performed in the present subsection. The geometrical model used in these simulations is similar to that presented in Figure 2, except for the slip systems, which are replaced by only one system. To obtain idealized slip and kink bands, this system is oriented such that the slip plane is either parallel (horizontal slip case) or normal (vertical slip case) to the shear direction.

The same microscopic and macroscopic boundary conditions as in Figure 2 are assumed. The material parameters are those presented in Table 1, except for $H_π$, which is taken as $H_π = -500$ MPa to model softening. To trigger strain localization, a material imperfection is considered in a small central region of height $H/11$. In this region, the initial first-order slip resistance $S_{π0}$ is linearly decreased to 0.9 $S_{π0}$ at the center.

Considering the case of horizontal slip, Figure 6 presents the associated results in terms of the stress-strain curves and distributions of plastic scalar variable $\gamma$ along the $x_2$-direction for different mesh refinements. Although these results are obtained considering only higher-order energetic effects ($l_{en}/H = 0.1$ and $l_{dis} = 0$), it has been verified that activation of higher-order dissipative effects ($l_{dis} \neq 0$) do not affect them. The SGCP model failed to regularize the obtained idealized slip band, which shows mesh-dependence. As can be seen in Figure 6b, the width of the obtained band is set by the mesh spacing and decreases with increasing the number of finite elements.

For the considered slip configuration, the present SGCP model behaves as a conventional crystal plasticity model. Indeed, slip gradients in the direction normal to the slip plane do not correspond to GNDs [6,23,43]. Consequently, no size effects are predicted by GND-based SGCP theories for the development of such normal slip gradients. This is also
the case for Cosserat crystal plasticity theories in which size effects are only triggered by
lattice curvature [6].

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Results of the horizontal slip case: (a) stress–strain curves and (b) distributions of plastic scalar variable $\gamma$ along the $x_2$-direction at $\Gamma = 0.0032$ for different mesh refinements ($l_{en}/H = 0.1$, $l_{dis} = 0$, $H_{\pi} = -500$ MPa, $n = 2$).</figcaption></figure>

Figure 7 presents the shear results associated with the case of vertical slip in terms of the stress–strain responses and distributions of plastic scalar variable $\gamma$ along the */++$x_2$–direction for different mesh refinements. Using a sufficiently fine mesh, the stress–strain response becomes nearly mesh-independent (Figure 7a). Note that the slope inversion and the plateau-shaped form observed on the stress–strain curves beyond $\Gamma = 0.008$ is due to the competition between softening caused by negative $H_{\pi}$ and hardening caused by higher-order energetic effects ($l_{en} \neq 0$). As shown in Figure 7b, thick kink bands that spread over several finite elements are obtained. For a relatively fine mesh, the thickness of these bands is nearly mesh-insensitive.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Results of the vertical slip case: (a) stress–strain curves and (b) distributions of plastic scalar variable $\gamma$ along $x_2$-direction at $\Gamma = 0.01$ for different mesh refinements ($l_{en}/H = 0.1$, $l_{dis} = 0$, $H_{\pi} = -500$ MPa, $n = 2$).</figcaption></figure>

For Cosserat crystal plasticity models [6], the proposed SGCP model enables mesh-independent simulation of kink bands, which imply slip gradients in the slip direction. On the contrary, it fails to regularize idealized slip bands, which are not accompanied by the development of GNDs. However, the idealized slip conditions used to obtain the latter bands are not commonly encountered in real material deformation. The following subsection discusses the ability of the present SGCP model to correctly capture strain localization within a more practical problem, for which shear bands do not necessarily have particular orientations with respect to the initial crystallographic directions.
### 5.1.2. Uniaxial Tension with Two Slip Systems: Regularization of General Localization Bands

A uniaxial tension problem with two slip systems symmetrically oriented with respect to the tension direction is considered in this subsection. Figure 8 presents the associated geometrical model, which is a plate of dimensions $2L \times 2W = 20 \times 10$ mm having two slip systems inclined at 60° from both sides of $x_1$-direction ($\theta_1 = -\theta_2 = 60°$). For this problem, according to bifurcation analysis [3,10], it is expected that localization will occur along two crossed shear bands whose orientations are close to those of the slip planes. Due to symmetry, only one quarter (upper right part) of the plate is analyzed. The macroscopic boundary conditions applied on this portion of the plate are set as:

$$ \left\{ \begin{array}{ll} u_1 = 0 & \text{for } x_1 = 0 \\ u_2 = 0 & \text{for } x_2 = 0 \\ u_2 = U & \text{for } x_2 = L \end{array} \right. \qquad (52) $$

Furthermore, all the edges of this plate portion are assumed to be micro-free (transparent to dislocations). As in the previous simple shear problem, $H_\pi$ is assumed to be negative to model softening ($H_\pi = -500$ MPa). The other material constants are given in Table 3, following several previous works designated for single crystal tension problems [11,15]. To trigger localization, a material imperfection is considered within a square region of 0.25 mm side located at the bottom left-hand corner of the considered plate portion (Figure 8). In this region, the initial first-order slip resistance $S_{\pi 0}$ is replaced by 0.9 $S_{\pi 0}$. To quantify strain localization, the accumulated plastic slip variable defined by $\gamma_{acc} = \sum_{\alpha} |\gamma^{\alpha}|$ is adopted.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Uniaxial tension of a single crystal plate with two slip systems symmetrically oriented with respect to the tension direction.</figcaption></figure>

<table><caption>**Table 3.** The material parameters implied for uniaxial tension of a single crystal plate, inspired by [11,15].</caption><thead><tr><th>Material Parameter Name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i></td><td>60.84</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i></td><td>0.3</td><td>—</td></tr><tr><td>Reference slip rate</td><td>γ̇<sub>0</sub><sup>α</sup></td><td>0.001</td><td>s<sup>-1</sup></td></tr><tr><td>Rate-sensitivity parameter</td><td><i>m</i></td><td>0.005</td><td>—</td></tr><tr><td>Higher-order energetic slip resistance</td><td><i>X</i><sub>0</sub></td><td>60.84</td><td>MPa</td></tr><tr><td>Initial first-order dissipative slip resistance</td><td><i>S</i><sub>π0</sub></td><td>60.84</td><td>MPa</td></tr><tr><td>Initial higher-order dissipative slip resistance</td><td><i>S</i><sub>ξ0</sub></td><td>case study</td><td>MPa</td></tr><tr><td>First-order hardening modulus</td><td><i>H</i><sub>π</sub></td><td>case study</td><td>MPa</td></tr><tr><td>Higher-order hardening modulus</td><td><i>H</i><sub>ξ</sub></td><td>case study</td><td>MPa</td></tr><tr><td>Energetic length scale</td><td><i>l</i><sub>en</sub></td><td>case study</td><td>μm</td></tr><tr><td>Dissipative length scale</td><td><i>l</i><sub>dis</sub></td><td>case study</td><td>μm</td></tr><tr><td>Defect energy index</td><td><i>n</i></td><td>2</td><td>—</td></tr></tbody></table>
Figure 9 presents the associated results in terms of nominal stress–strain curves and accumulated plastic slip variable distributions along the $x_2$-direction for different mesh refinements. The nominal stress is calculated as the reaction force $F$ on the top edge over the initial width $W_0$ of the studied plate portion (the thickness is assumed to be 1) and the nominal strain as the ratio between the elongation $\Delta L$ and the initial length $L_0$. The values of $\gamma_{acc}$ in Figure 9b are determined at nodes belonging to the center line of the studied portion of the plate (dashed blue line in Figure 8).

Figure 10 presents the deformed meshes and contours of the accumulated plastic slip variable $\gamma_{acc}$ for different mesh refinements. Concerning the nominal stress–strain response, mesh convergence is obtained as soon as a sufficiently fine mesh is applied (Figure 9a). In terms of strain localization, it seems that finer meshes lead to a larger $\gamma_{acc}$ peak value in the obtained shear band. However, this band presents a finite width, which converges toward an intrinsic value covering several finite elements.

Based on Figure 10, such a band is oriented at an angle of approximately 54° with respect to the $x_1$-direction, which is slightly smaller than the angle of the first slip direction ($\theta_1 = 60^\circ$). This is in agreement with the bifurcation and numerical results of Kuroda [11]. Using a similar configuration of activated slip systems with $\theta_1 = 54.7^\circ$, these authors demonstrated that strain localization occurs along a line inclined at approximately 49° from the $x_1$-direction.

To better analyze the obtained results, the lattice rotation is approximated based on the polar decomposition of the elastic deformation gradient $F^e = R^e U^e$, where $R^e$ and $U^e$ are, respectively, the elastic rotation and elastic stretch tensors. Neglecting the elastic distortion described by $U^e$, $R^e$ can be interpreted as the lattice rotation, and the corresponding angle $\theta$ is calculated using [8]:

$$ \theta = \arccos\left[\frac{1}{2}(\text{tr} R^e - 1)\right] \qquad (53) $$

The angle sign is determined considering the sign of the off-diagonal terms of $R^e$. By doing so, we find that strain localization is accompanied by a lattice rotation of approximately $-5.5^\circ$. Therefore, the localization band obtained in the considered portion of the plate is not an idealized slip band, as it is associated with lattice curvature, which explains its regularization by the present SGCP model. However, the obtained band could be considered as a slip band since it is nearly parallel to the rotated first slip direction within the band ($60^\circ - 5.5^\circ = 54.5^\circ$).

Although GND-based SGCP theories fail to regularize idealized slip bands, the present investigation shows that they are capable of regularizing general localization bands. The next subsection aims at investigating the ability of such theories to correctly capture complex localization behaviors with competition between the kink and slip bands.

<figure><img src="image_12.png" /><figcaption>Figure 9. The results of uniaxial tension with two slip systems: (a) nominal stress–strain curves and (b) accumulated plastic slip variable distributions along the $x_2$-direction at $\Delta L/L_0 = 0.1$ for different mesh refinements ($l_{en}/W = 0.2$, $l_{dis} = 0$, $H_\pi = -500$ MPa).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10. Deformed meshes and contours of accumulated plastic slip variable $\gamma_{acc}$ at nominal strain $\Delta L/L_0 = 0.1$ for different mesh refinements (a) mesh 30 × 60 (b) mesh 40 × 80 and (c) mesh 50 × 100 ($l_{en}/W = 0.2$, $l_{dis} = 0$, $H_{\pi} = -500$ MPa).</figcaption></figure>

## 5.2. Investigation of Localization Modes in Single Crystals Undergoing Single Slip

The present finite deformation SGCP model is applied, hereafter, to reproduce the uniaxial tension problem of a single crystal undergoing single slip, which has been studied by Ling et al. [29] using a micromorphic crystal plasticity approach. This problem is characterized by complex localization behavior with competition between the kink and slip bands.

The geometrical model is similar to that used by Ling et al. [29]. It consists of a single crystal plate of dimensions W and L in, respectively, the $x_1$- and $x_2$- directions with an aspect ratio of $L/W = 6$ and a single slip system oriented to an angle of $\theta = -33.7^\circ$ with respect to $x_1$-direction. This plate is subjected to the following macroscopic boundary conditions:

$$ \left. \begin{array}{l} u_1 = 0, u_2 = 0 \quad \text{for} \quad x_2 = 0 \\ u_1 = 0, u_2 = U \quad \text{for} \quad x_2 = L \end{array} \right\} \qquad (54) $$

Furthermore, all the plate edges are assumed to be micro-free (transparent to dislocations). The material parameters are the same as those used in the previous uniaxial tension problem (Table 3), except for $H_{\pi}$, which is set to a small positive value ($H_{\pi} = 1$ MPa). For this problem, it is not needed to introduce material softening or imperfection to trigger strain localization. Indeed, the latter can occur due to geometrical softening induced by cross-section reduction and lattice rotation.

Although the material parameters used in the present and [29] investigations are not the same, comparison between the associated results is possible from a qualitative point of view, which is sufficient to meet the objectives of this subsection. As in [29], the plastic scalar variable defined by (37) on the considered slip system is adopted to investigate the strain localization and formation of slip and kink bands. Since the main difference between these bands is that kink bands are accompanied by a lot of lattice rotation, the lattice rotation angle (53) is used to detect the presence of the latter bands.

### 5.2.1. Regularization of the Mesh-Dependence of the Stress-Strain Response

To check the mesh-regularization capability of the present SGCP model regarding the stress-strain response, different mesh refinements are considered: 8 × 48, 16 × 96, and 24 × 144 meshes. Figure 11 presents the associated results obtained using different constitutive configurations. In this figure, the CP, SGCP-EN, and SGCP-DIS curves are, respectively, associated with the conventional CP ($l_{en} = l_{dis} = 0$), SGCP including only
higher-order energetic effects ($l_{en}/W = 0.25$ and $l_{dis}/W = 0$), and SGCP including only higher-order dissipative effects ($l_{en}/W = 0$ and $l_{dis}/W = 0.25$). For conventional crystal plasticity, the results are highly mesh-dependent and display spurious oscillations. Similar results were obtained by [29].

Activation of higher-order effects, regardless of their energetic or dissipative nature, allows for overcoming these mesh difficulties. No spurious oscillations were observed for the results associated with SGCP-EN and SGCP-DIS cases. Furthermore, mesh convergence is obtained for both cases as soon as a sufficiently refined mesh is used. The results associated with 16 × 96 and 24 × 144 meshes are in very close agreement for both these cases. The apparent increase of yield stress obtained for the SGCP-DIS case is due to the strengthening effects of the dissipative length scale $l_{dis}$, as widely recognized in the literature [34,43]. In the following subsections, the considered plate is discretized using a 24 × 144 mesh.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Uniaxial tension results obtained using conventional CP (CP: $l_{en} = l_{dis} = 0$), SGCP including only higher-order energetic effects (SGCP-EN: $l_{en}/W = 0.25$, $l_{dis} = 0$) and SGCP including only higher-order dissipative effects (SGCP-DIS: $l_{en} = 0$, $l_{dis}/W = 0.25$) for different meshes, with $S_{\zeta 0} = 60.84$ MPa and $H_{\zeta} = 0$.</figcaption></figure>

### 5.2.2. Localization Modes Using Conventional Crystal Plasticity

For comparison purposes, localization modes in the considered plate are first revisited using conventional CP ($l_{en} = l_{dis} = 0$). Figure 12a shows the evolution of the plastic slip variable $\gamma$ (37) within this plate at different nominal strains. To investigate the formation of kink bands, evolution of the lattice rotation angle within the plate is presented in Figure 12b. At the beginning of straining ($\Delta L/L_0 \approx 0.005$), both slip and kink bands start to emerge (crossed bands near the top and bottom edges of the plate). This is in agreement with the bifurcation analysis of Asaro and Rice [4].

According to this analysis, slip and kink bands can simultaneously appear at incipient plasticity when using conventional CP. By increasing $\Delta L/L_0$, the kinking phenomenon becomes dominant, leading to the formation of two pronounced kink bands. These bands are accompanied by considerable lattice rotation, which results in a local increase of the Schmid factor and then a local decrease of the resolved shear stress (RSS) for the considered slip system. This geometrical softening mechanism acting in the kink regions promotes the development of kink bands.

As straining increases, strain localizes in the kink band located near the bottom edge of the plate. In addition, the direction of localization starts to shift gradually, which is due to formation of a slip band crossing the initial kink region. Due to hardening effects, the RSS in the kink region increases and can reach a level for which the glide along the “unrotated” slip direction becomes easier, hence, resulting in the appearance of the slip band.

It can be observed from Figure 12b that formation of this band does not give rise to lattice rotation. Finally, the formed slip band intensifies with straining, leading to final necking of the plate.
<figure><img src="image_15.png" /><figcaption>(a) Evolution of the plastic slip variable $\gamma$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Evolution of the lattice rotation angle $\theta$ in radians.</figcaption></figure>

**Figure 12.** Evolution of the plastic slip variable $\gamma$ and lattice rotation angle $\theta$ at different nominal tensile strains predicted by conventional CP ($l_{en} = l_{dis} = 0$).

### 5.2.3. Localization Modes Using Purely Higher-Order Energetic Effects

The above complex localization behavior is studied, hereafter, using the present SGCP model with only higher-order energetic effects. Figure 13 presents the contours of the plastic slip variable $\gamma$ and the lattice rotation angle $\theta$ within the studied plate at different nominal strains, with the energetic length scale $l_{en} = 0.25$ W. Compared to conventional crystal plasticity results, two broader kink bands appear at the beginning of straining ($\Delta L/L_0 = 0.006$ in Figure 13a). These bands rapidly aggregate to form a larger kink region with concentrated lattice rotation at the center of the studied plate. In the meanwhile, a central slip band appears and starts to interact with this region as shown at straining $\Delta L/L_0 = 0.015$.
The kink-slip transition region shows a much smoother form compared to that obtained using conventional CP (Figure 12a). As straining increases, the formed slip band becomes more and more intensified, leading to final necking of the plate.

<figure><img src="image_17.png" /><figcaption>(a) Evolution of the plastic slip variable $\gamma$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Distribution of the lattice rotation angle $\theta$ in radians.</figcaption></figure>

**Figure 13.** Evolution of the plastic slip variable $\gamma$ and lattice rotation angle $\theta$ at different nominal tensile strains predicted by SGCP including only higher-order energetic effects ($l_{en}/W = 0.25$, $l_{dis} = 0$).

To investigate the influence of the higher-order energetic effects on the present localization behavior, different values of $l_{en}$ were tested. Figure 14 presents the associated results in terms of plastic slip contours in the studied plate at straining $\Delta L/L = 0.006$, which corresponds to the emergence of kink bands. The width of the obtained bands as well as the distance between them are functions of $l_{en}$. By increasing this parameter, larger and closer kink bands are obtained. For sufficiently large $l_{en}$, a single kink band located at the center of the considered plate is formed.
Figure 15 presents the nominal tensile stress–strain curves and the plastic slip variable distributions along the plate length direction for different energetic length scales. Although no change in the initial yield is observed, the stress level after yielding increases with increasing $l_{en}$ due to strain gradient hardening effects (Figure 15a). In addition, larger $l_{en}$ leads to a larger final localization band with lower peak value of plastic slip $\gamma$ in this band (Figure 15b). This is in accordance with several investigations based on nonlocal approaches [26,29,53,58,59]. The energetic length scale behaves as a “localization limiter” and governs the width of localization bands.

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Influence of energetic length scale $l_{en}$ on the formation of the kink bands at the beginning of straining $\Delta L/L_0 = 0.006$ ($l_{dis} = 0$).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Influence of energetic length scale $l_{en}$: (a) nominal tensile stress–strain curves and (b) distribution of plastic slip variable $\gamma$ along the $x_2$ direction at $\Delta L/L_0 = 0.2$ ($l_{dis} = 0$).</figcaption></figure>

### 5.2.4. Localization Modes Using Purely Higher-Order Dissipative Effects

The present subsection aims at investigating the higher-order dissipative effects on the localization behavior in the considered plate, while neglecting the higher-order energetic effects ($l_{en} = 0$). One distinctive feature of the present SGCP model is that the higher-order dissipative effects are uncoupled from the first-order ones. The advantages of this uncoupled constitutive assumption have already been discussed in [34]. It is worth recalling here that such an assumption makes it possible to flexibly choose the initial value of the higher-order dissipative slip resistance $S_{\zeta 0}$.
This parameter can also be set to zero in order to obtain progressive activation of the higher-order dissipative effects, as with the energetic ones involved in the previous subsection. In this case, a nonzero higher-order hardening modulus $H_\zeta$ is necessary to activate such effects according to (34). Using $S_{\zeta0} = 0$ and $H_\zeta = 60.84$ MPa, the obtained results (not presented in the paper) are qualitatively quite similar to those obtained in the above subsection.

Applying nonzero $S_{\zeta0}$, finite higher-order dissipative stresses take place and instantaneously activate with the emergence of strain gradients. The influence of these stresses is investigated, hereafter, assuming $S_{\zeta0} = S_{\pi0} = 60.84$ MPa, as in classical SGCP theories with coupled dissipation, while keeping $H_\zeta = 60.84$ MPa. Figure 16 presents the associated results in terms of the plastic slip variable and lattice rotation distributions in the plate at different nominal tensile strains.

These results are quite different from those obtained using conventional CP and SGCP, including progressive activation of higher-order effects. Instead of two kink bands, a large central region parallel to the slip direction appears at the beginning of straining ($\Delta L/L_0 = 0.0035$) with a nearly plateau-shaped distribution of $\gamma$. The width of this region progressively decreases with straining.

A less-pronounced kink region appears at the plate center and interacts with the initial slip region, as can be shown by analyzing the lattice rotation distributions in the plate (Figure 16b). The width of both these regions continues to reduce with increasing $\Delta L/L_0$ until the final necking of the plate. The final localization band is oriented to an angle close to that of the slip direction. Compared to those obtained in the previous subsections, this band is much more diffuse. Using higher-order dissipative effects with nonzero $S_{\zeta0}$, slip localization seems to outweigh kinking, although the latter phenomenon is not completely removed.

<figure><img src="image_21.png" /><figcaption>**Figure 16.** *Cont.*</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 16.** Evolution of the plastic slip variable $\gamma$ and lattice rotation angle $\theta$ at different nominal tensile strains predicted by SGCP, including only higher-order dissipative effects with nonzero initial higher-order slip resistance $S_{\tilde{\xi}0}$ ($l_{en} = 0$, $l_{dis}/W = 0.25$, $S_{\tilde{\xi}0} = 60.84$ MPa, $H_{\tilde{\xi}} = 60.84$ MPa).</figcaption></figure>

Effects of $l_{dis}$ on the localization response of the considered plate are given in Figure 17, which is obtained assuming $S_{\tilde{\xi}0} = S_{\pi 0} = 60.84$ MPa and $H_{\tilde{\xi}} = 60.84$ MPa. This parameter leads to similar effects as for $l_{en}$ regarding the increase of hardening after yielding and the growth of the width of the final localization band. In addition to these effects, $l_{dis}$ presents strengthening effects, which are illustrated by an increase of the initial yield.

<figure><img src="image_23.png" /><figcaption>**Figure 17.** Influence of dissipative length scale $l_{dis}$: (a) nominal tensile stress–strain response and (b) distribution of the plastic slip variable $\gamma$ along the $x_2$ direction at $\Delta L/L_0 = 0.4$ ($l_{en} = 0$, $S_{\tilde{\xi}0} = 60.84$ MPa, $H_{\tilde{\xi}} = 60.84$ MPa).</figcaption></figure>

**6. Conclusions**

This contribution focused on the capabilities of strain gradient crystal plasticity (SGCP) in modeling localization phenomena in single crystals, considering both higher-order energetic and dissipative effects. To this end, the Gurtin-type SGCP model pre-
sented in Jebahi et al. [34] was extended to finite deformation following the procedure of Gurtin [21,23].

The capabilities of this model in regularizing idealized slip and kink bands were first revisited using a simple shear problem with a single slip system. The obtained results show the effectiveness of the model in capturing kink bands. However, this is not the case for idealized slip bands. Indeed, slip gradients in the direction normal to the slip plane do not correspond to GNDs.

Consequently, the present model behaves as a conventional crystal plasticity model with regard to idealized slip bands. To investigate the effectiveness of the present model in dealing with more general localization bands, it was applied to simulate a uniaxial tension problem with double slip systems symmetrically oriented with respect to the tension direction. For this problem, regularized localization bands slightly rotated with respect to the slip direction were obtained, showing the ability of the applied model to regularize general slip bands.

The good capabilities of the applied model in dealing with the latter bands motivated the application of this model to study the complex localization behavior with competition between slip and kink bands, which is encountered within single crystals undergoing single slip. A uniaxial tension problem of a single crystal plate having a single slip system was simulated using the present SGCP model. Different constitutive configurations, including purely higher-order energetic effects and purely higher-order dissipative effects, were considered. The associated results show that mesh-dependence difficulties, encountered using conventional crystal plasticity, can be removed by including higher-order effects, either energetic or dissipative.

Considering purely higher-order energetic effects, the results of the competition between slip and kink bands are qualitatively quite similar to those obtained by [29] using the micromorphic approach. At the beginning of straining, two thick kink bands appear and rapidly aggregate to form a larger kink region with concentrated lattice rotation at the center of the studied plate. The width and the distance between the formed kink bands depend on the energetic length scale $l_{en}$. For sufficiently large $l_{en}$, a single kink band at the center of the plate is formed at the beginning of straining. As straining increases, a central slip band appears and starts to interact with the kink region. With further straining, the formed slip band becomes increasingly intensified, leading to final necking of the plate.

Concerning higher-order dissipative effects, distinct results were obtained depending on the value of the initial higher-order dissipative slip resistance $S_{\zeta0}$. The latter parameter can easily and independently be controlled in the present SGCP model due to the implied uncoupled dissipation assumption. Using zero $S_{\zeta0}$ (representing the progressive activation of the higher-order dissipative effects, as it is the case for the higher-order energetic ones), the obtained results are qualitatively quite similar to those obtained using purely higher-order energetic effects.

However, the use of a nonzero $S_{\zeta0}$ (as in classical Gurtin-type SGCP theories including higher-order dissipation) leads to quite different results. This choice leads to a localization behavior that is overwhelmingly dominated by slip banding, although analysis of lattice rotation in the studied plate shows the coexistence of kinking.

**Author Contributions:** Conceptualization, L.C., M.J. and F.A.-M.; methodology, L.C., M.J. and F.A.-M.; investigation, L.C. and M.J.; analysis of the results and discussion: L.C., M.J. and F.A.-M.; writing—original draft preparation, L.C. and M.J.; writing—review and editing, M.J. and F.A.-M.; supervision, M.J. and F.A.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received a financial support from the China Scholarship Council (CSC201706280201).

**Conflicts of Interest:** The authors declare no conflict of interest.
## References

1. Considère, A. Memoire sur l'emploi du fer et de l'acier dans les constructions. *Ann. Ponts Chaussee* **1885**, *9*, 574-775.

2. Hill, R. On discontinuous plastic states, with special reference to localized necking in thin sheets. *J. Mech. Phys. Solids* **1952**, *1*, 19-30. [CrossRef]

3. Rice, J.R. The localization of plastic deformation. In *Proceedings of the 14th International Congress on Theoretical and Applied Mechanics*, Delft, The Netherlands, 30 August-4 September 1976; pp. 207-220.

4. Asaro, R.J.; Rice, J.R. Strain localization in ductile single crystals. *J. Mech. Phys. Solids* **1977**, *25*, 309-338. [CrossRef]

5. Forest, S.; Cailletaud, G. Strain localization in single crystals: Bifurcation analysis, effects of boundaries and interfaces. *Eur. J. Mech. A Solids* **1995**, *14*, 747-771.

6. Forest, S. Modeling slip, kink and shear banding in classical and generalized single crystal plasticity. *Acta Mater.* **1998**, *46*, 3265-3281. [CrossRef]

7. Forest, S.; Lorentz, E. Localization phenomena and regularization methods. *Local Approach Fract.* **2004**, *1*, 311-371.

8. Marano, A.; Gélébart, L.; Forest, S. Intragranular localization induced by softening crystal plasticity: Analysis of slip and kink bands localization modes from high resolution FFT-simulations results. *Acta Mater.* **2019**, *175*, 262-275. [CrossRef]

9. Marano, A.; Lionel, G.; Gélébart, L.; Forest, S. FFT-based simulations of slip and kink bands formation in 3D polycrystals: Influence of strain gradient crystal plasticity. *J. Mech. Phys. Solids* **2021**, 104295. [CrossRef]

10. Needleman, A.; Tvergaard, V. Analyses of plastic flow localization in metals. *Appl. Mech. Rev.* **1992**, *45*, S3-S18. [CrossRef]

11. Kuroda, M. On large-strain finite element solutions of higher-order gradient crystal plasticity. *Int. J. Solids Struct.* **2011**, *48*, 3382-3394. [CrossRef]

12. Orowan, E. A type of plastic deformation new in metals. *Nature* **1942**, *149*, 643-644. [CrossRef]

13. Di Gioacchino, F.; Quinta da Fonseca, J. An experimental study of the polycrystalline plasticity of austenitic stainless steel. *Int. J. Plast.* **2015**, *74*, 92-109. [CrossRef]

14. Peirce, D.; Asaro, R.J.; Needleman, A. An analysis of nonuniform and localized deformation in ductile single crystals. *Acta Metall.* **1982**, *30*, 1087-1119. [CrossRef]

15. Peirce, D.; Asaro, R.J.; Needleman, A. Material rate dependence and localized deformation in crystalline solids. *Acta Metall.* **1983**, *31*, 1951-1976. [CrossRef]

16. Zhang, M.; Bridier, F.; Villechaise, P.; Mendez, J.; McDowell, D.L. Simulation of slip band evolution in duplex Ti-6Al-4V. *Acta Mater.* **2010**, *58*, 1087-1096. [CrossRef]

17. Takagi, K.; Mayama, T.; Mine, Y.; Chiu, Y.L.; Takashima, K. Extended ductility due to kink band formation and growth under tensile loading in single crystals of Mg-Zn-Y alloy with 18R-LPSO structure. *J. Alloys Compd.* **2019**, *806*, 1384-1393. [CrossRef]

18. Needleman, A. Material rate dependence and mesh sensitivity in localization problems. *Comput. Methods Appl. Mech. Eng.* **1988**, *67*, 69-85. [CrossRef]

19. Aifantis, E.C. on the Microstructural Origin of Certain Inelastic Models. *Am. Soc. Mech. Eng.* **1984**, *106*, 326-330. [CrossRef]

20. Acharya, A.; Bassani, J.L. Lattice incompatibility and a gradient theory of crystal plasticity. *J. Mech. Phys. Solids* **2000**, *48*, 1565-1595. [CrossRef]

21. Gurtin, M.E. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations. *J. Mech. Phys. Solids* **2002**, *50*, 5-32. [CrossRef]

22. Gurtin, M.E. A theory of grain boundaries that accounts automatically for grain misorientation and grain-boundary orientation. *J. Mech. Phys. Solids* **2008**, *56*, 640-662. [CrossRef]

23. Gurtin, M.E. A finite-deformation, gradient theory of single-crystal plasticity with free energy dependent on densities of geometrically necessary dislocations. *Int. J. Plast.* **2008**, *24*, 702-725. [CrossRef]

24. Pourayevali, H.; Xu, B. A Hardening Description based on a Finite-Deformation Gradient Crystal Plasticity Model: Formulation and Numerical Implementation. *Pamm* **2015**, *15*, 343-344. [CrossRef]

25. Wulfinghoff, S.; Böhlke, T. Gradient crystal plasticity including dislocation-based work-hardening and dislocation transport. *Int. J. Plast.* **2015**, *69*, 152-169. [CrossRef]

26. Kuroda, M.; Needleman, A. Nonuniform and localized deformation in single crystals under dynamic tensile loading. *J. Mech. Phys. Solids* **2019**, *125*, 347-359. [CrossRef]

27. Berbenni, S.; Taupin, V.; Lebensohn, R.A. A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals. *J. Mech. Phys. Solids* **2020**, *135*, 1-23. [CrossRef]

28. Jebahi, M.; Forest, S. Scalar-based strain gradient plasticity theory to model size-dependent kinematic hardening effects. *Contin. Mech. Thermodyn.* **2021**, 1-23. [CrossRef]

29. Ling, C.; Forest, S.; Besson, J.; Tanguy, B.; Latourte, F. A reduced micromorphic single crystal plasticity model at finite deformations. Application to strain localization and void growth in ductile metals. *Int. J. Solids Struct.* **2018**, *134*, 43-69. [CrossRef]

30. Acharya, A.; Bassani, J.L.; Beaudoin, A. Geometrically necessary dislocations, hardening, and a simple gradient theory of crystal plasticity. *Scr. Mater.* **2003**, *48*, 167-172. [CrossRef]

31. Liu, Z.L.; Zhuang, Z.; Liu, X.M.; Zhao, X.C.; Zhang, Z.H. A dislocation dynamics based higher-order crystal plasticity model and applications on confined thin-film plasticity. *Int. J. Plast.* **2011**, *27*, 201-216. [CrossRef]

32. Bittencourt, E. On the effects of hardening models and lattice rotations in strain gradient crystal plasticity simulations. *Int. J. Plast.* **2018**, *108*, 169-185. [CrossRef]
33. Bittencourt, E. Interpretation of the size effects in micropillar compression by a strain gradient crystal plasticity theory. *Int. J. Plast.* **2019**, *116*, 280–296. [CrossRef]

34. Jebahi, M.; Cai, L.; Abed-Meraim, F. Strain gradient crystal plasticity model based on generalized non-quadratic defect energy and uncoupled dissipation. *Int. J. Plast.* **2020**, *126*, 102617. [CrossRef]

35. Scherer, J.M.; Besson, J.; Forest, S.; Hure, J.; Tanguy, B. Strain gradient crystal plasticity with evolving length scale: Application to voided irradiated materials. *Eur. J. Mech. A/Solids* **2019**, *77*, 103768. [CrossRef]

36. Kröner, E. Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen. *Arch. Ration. Mech. Anal.* **1959**, [CrossRef]

37. Lee, E.H. Elastic-plastic deformation at finite strains. *J. Appl. Mech. Trans. ASME* **1964**, *36*, 1–6. [CrossRef]

38. Gurtin, M.E. The Burgers vector and the flow of screw and edge dislocations in finite-deformation single-crystal plasticity. *J. Mech. Phys. Solids* **2006**, *54*, 1882–1898. [CrossRef]

39. Forest, S.; Aifantis, E.C. Some links between recent gradient thermo-elasto-plasticity theories and the thermomechanics of generalized continua. *Int. J. Solids Struct.* **2010**, *47*, 3367–3376. [CrossRef]

40. Forest, S. Strain Localization Phenomena in Generalized Single Crystal Plasticity. *J. Mech. Behav. Mater.* **2011**, *11*, 45–50. [CrossRef]

41. Germain, P. La méthode des puissances virtuelles en mécanique des milieux continus, première partie: Théorie du second gradient. *J. Mec.* **1973**, *12*, 235–274.

42. Forest, S. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. *J. Eng. Mech.* **2009**, *135*, 117–131. [CrossRef]

43. Gurtin, M.E.; Anand, L.; Lele, S.P. Gradient single-crystal plasticity with free energy dependent on dislocation densities. *J. Mech. Phys. Solids* **2007**, *55*, 1853–1878. [CrossRef]

44. Voyiadjis, G.Z.; Deliktas, B. Mechanics of strain gradient plasticity with particular reference to decomposition of the state variables into energetic and dissipative components. *Int. J. Eng. Sci.* **2009**, *47*, 1405–1423. [CrossRef]

45. Fleck, N.A.; Willis, J.R. A mathematical basis for strain-gradient plasticity theory-Part I: Scalar plastic multiplier. *J. Mech. Phys. Solids* **2009**, *57*, 161–177. [CrossRef]

46. Gurtin, M.E. A finite-deformation, gradient theory of single-crystal plasticity with free energy dependent on the accumulation of geometrically necessary dislocations. *Int. J. Plast.* **2010**, *26*, 1073–1096. [CrossRef]

47. Gurtin, M.E.; Ohno, N. A gradient theory of small-deformation, single-crystal plasticity that accounts for GND-induced interactions between slip systems. *J. Mech. Phys. Solids* **2011**, *59*, 320–343. [CrossRef]

48. Scherer, J.M.; Phalke, V.; Besson, J.; Forest, S.; Hure, J.; Tanguy, B. Lagrange multiplier based vs micromorphic gradient-enhanced rate-(in)dependent crystal plasticity modelling and simulation. *Comput. Methods Appl. Mech. Eng.* **2020**, *372*, 113426. [CrossRef]

49. Kuroda, M.; Tvergaard, V. A finite deformation theory of higher-order gradient crystal plasticity. *J. Mech. Phys. Solids* **2008**, *56*, 2573–2584. [CrossRef]

50. Klusemann, B.; Bargmann, S.; Svendsen, B. Two models for gradient inelasticity based on non-convex energy. *Comput. Mater. Sci.* **2012**, *64*, 96–100. [CrossRef]

51. Pourayevali, H.; Xu, B.X. A study of gradient strengthening based on a finite-deformation gradient crystal-plasticity model. *Contin. Mech. Thermodyn.* **2017**, *29*, 1–24. [CrossRef]

52. Peng, X.L.; Husser, E.; Huang, G.Y.; Bargmann, S. Modeling of surface effects in crystalline materials within the framework of gradient crystal plasticity. *J. Mech. Phys. Solids* **2018**, *112*, 508–522. [CrossRef]

53. Borg, U. Strain gradient crystal plasticity effects on flow localization. *Int. J. Plast.* **2007**, *23*, 1400–1416. [CrossRef]

54. Asaro, R.J. Elastic-plastic memory and kinematic-type hardening. *Acta Metall.* **1975**, *23*, 1255–1265. [CrossRef]

55. Wulfinghoff, S.; Forest, S.; Böhlke, T. Strain gradient plasticity modeling of the cyclic behavior of laminate microstructures. *J. Mech. Phys. Solids* **2015**, *79*, 1–20. [CrossRef]

56. El-Naaman, S.A.; Nielsen, K.L.; Niordson, C.F. An investigation of back stress formulations under cyclic loading. *Mech. Mater.* **2019**, *130*, 76–87. [CrossRef]

57. Forest, S.; Boubidi, P.; Sievert, R. Strain localization patterns at a crack tip in generalized single crystal plasticity. *Scr. Mater.* **2001**, *44*, 953–958. [CrossRef]

58. Sluys, L.J.; Estrin, Y. Analysis of shear banding with a dislocation based gradient plasticity model. *Int. J. Solids Struct.* **2000**, *37*, 7127–7142. [CrossRef]

59. Niordson, C.F.; Redanz, P. Size-effects in plane strain sheet-necking. *J. Mech. Phys. Solids* **2004**, *52*, 2431–2454. [CrossRef]