# A FFT-based mesoscale continuum dislocation mechanics with defect energy: Applications to composites and polycrystals

Stéphane Berbenni, Vincent Taupin, Ricardo A Lebensohn

► To cite this version:

Stéphane Berbenni, Vincent Taupin, Ricardo A Lebensohn. A FFT-based mesoscale continuum dislocation mechanics with defect energy: Applications to composites and polycrystals. European Journal of Mechanics - A/Solids, 2025, 111, pp.105548. 10.1016/j.euromechsol.2024.105548. hal-04871972

HAL Id: hal-04871972

https://hal.univ-lorraine.fr/hal-04871972v1

Submitted on 7 Jan 2025

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
Full length article

# A FFT-based mesoscale continuum dislocation mechanics with defect energy: Applications to composites and polycrystals

Stéphane Berbenni<sup>a, c, ∗</sup>, Vincent Taupin<sup>a, c, †</sup>, Ricardo A. Lebensohn<sup>b, ‡</sup>

ARTICLE INFO

Keywords:

Dislocation mechanics

Backstress

Crystal plasticity

Composites

Polycrystals

EVPFFT

A B S T R A C T

A crystal plasticity elastoviscoplastic FFT (fast Fourier transform) formulation with a mesoscale continuum field dislocation mechanics model is presented, which incorporates a defect energy density that depends on GND densities and an associated material length scale. This allows to thermodynamically derive internal length scale dependent intra-crystalline backstress and Peach–Koehler force acting on GND densities. The model considers GND density evolution through a filtered numerical spectral approach, which is coupled with stress equilibrium through the elastoviscoplastic FFT algorithm. The discrete Fourier transform (DFT) method together with finite difference (FD) schemes is applied to solve both the backstress tensor and the Fourier–Green operator. Numerical results are first reported for two-phase laminate composites with plastic single crystal channels and elastic precipitates for shear loadings. Channel size effects are simulated and analyzed on the overall and local hardening behaviors during monotonous loadings. In addition, the evolutions of GND densities and the role of their associated backstress on size effects are examined during reversible shear loading. In a second part, the role of the defect energy internal length scale on polycrystal’s hardening during tension–compression is discussed. The results are compared to those obtained using FFT-based continuum field dislocation mechanics without defect energy.

1. Introduction

Experimental works evidenced two different size effects in engi-
neering materials. The first category is due to loading associated with
non homogeneous boundary conditions like torsion (Fleck et al., 1994),
microbending (Stölken and Evans, 1998) and nanoindentation (Nix and
Gao, 1998; Voyiadjis and Peters, 2010; Kucharski et al., 2024). The
second category, which will be more considered in the present study,
is due to plastic strain gradients developing in the microstructure in
the presence of phase or grain boundaries for monotonic and non-
monotonic loadings (uniaxial tension, pure shear, tension–compression,
etc.), i.e. homogeneous boundary conditions. These plastic strain gra-
dients generally lead to a smaller is stronger size effect. For example,
for dispersion-hardened alloys, it was reported that the flow stress
increases as the average particle size and the average particle spac-
ing decrease (Ashby, 1970; Lloyd, 1994; Nan and Clarke, 1996). For
polycrystalline materials, the flow stress is increased when the grain
size is decreased following the Hall–Petch relationship (Hall, 1951;
Petch, 1953; Cordero et al., 2016). In polycrystalline materials, the

lattice curvatures, i.e. the spatial gradients of crystallographic orienta-
tions, can be measured at the intra-granular level with Electron Back-
Scattered Diffraction (EBSD) (Pantleon, 2008; Allain-Bonasso et al.,
2012; Konijnenberg et al., 2015; Wallis et al., 2016; Genée et al., 2021).

Size dependent material behaviors were predicted by discrete meth-
ods like Discrete Dislocation Dynamics (DDD) first developed by Kubin
et al. (1992) and by Van der Giessen and Needleman (1995). DDD,
as a computational method, can naturally account for internal length
scale effects considering discrete distributions of dislocations and their
dynamics. For example, plasticity size effects were studied in micro-
specimens (Motz et al., 2005, 2009; Kiener et al., 2010) and Bauschinger
effect was predicted (Bittencourt et al., 2003; Nicola et al., 2005;
Waheed et al., 2017; Amouzou-Adoun et al., 2023). However, these
methods are rather limited in terms of computational efficiency when
3D polycrystals are studied. This is why continuum-based mechanics
methods were investigated for structural calculations.

Since the pioneering work of Aifantis (1984, 1987), different strain
gradient plasticity (SGP) formulations, either lower-order or higher-
order, have been developed during the last decades to study size

<footnote><sup>a</sup> Université de Lorraine, CNRS, Arts et Métiers Paris Tech, LEM3, F-57000 Metz, France</footnote>

<footnote><sup>b</sup> Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87845, USA</footnote>

<footnote><sup>c</sup> Laboratory of Excellence on Design of Alloy Metals for low-mAss Structures (DAMAS), Université de Lorraine, France</footnote>

<footnote>* Corresponding author at: Université de Lorraine, CNRS, Arts et Métiers Paris Tech, LEM3, F-57000 Metz, France.
*E-mail address:* stephane.berbenni@univ-lorraine.fr (S. Berbenni).</footnote>
effects, because conventional crystal plasticity theories were not able to predict them. Lower-order crystal plasticity SGP theories, denoted “LO-SGP”, generally assume that plastic strain gradients only modify strain-hardening moduli incorporating both SSD (statistically stored dislocations) and GND (geometrically necessary dislocations) densities (Acharya and Bassani, 2000; Acharya and Beaudoin, 2000; Busso et al., 2000; Evers et al., 2002, 2004; Cheong et al., 2005; Han et al., 2005a,b), while higher-order crystal plasticity SGP theories, denoted “HO-SGP”, consider plastic strain gradients as new kinematic variables with new higher-order boundary conditions (Gurtin, 2002, 2004; Gudmundson, 2004; Gurtin and Anand, 2005, 2009; Kuroda and Tvergaard, 2008; Forest, 2008; Cordero et al., 2010, 2012; Niordson and Kysar, 2014; Mayeur and McDowell, 2014; Wulfinghoff et al., 2015; Fleck and Willis, 2015; Ryś and Petryk, 2018; Jebahi and Forest, 2023). In their paper, Fleck and Willis (2015) discussed in the framework of SGP theories, the role of GNDs in terms of energetic vs. dissipative hardening. For example, it was highlighted by these authors that beam size effects during the cyclic experiments of Motz et al. (2005) on copper single crystals, affect isotropic hardening (forest-type hardening) due to plastic strain gradients. This strengthening type and the associated size effect is attributed to dissipative rather than energetic hardening. Likewise, it is admitted that the Hall-Petch strengthening in polycrystals is not associated with kinematic hardening, which indicates a dissipative hardening process. In contrast, GND stored at grain boundaries in the form of dislocation pile-ups are at the origin of a backstress associated with energetic hardening (Fleck and Willis, 2015). In HO-SGP theory, this stems from the incorporation of a *defect energy* that contains a material internal length scale *l* (Gurtin, 2002), which was further assessed by comparisons with DDD simulation results (Bittencourt et al., 2003; Nicola et al., 2005).

An alternative continuum-based approach, called Phenomenological Mesoscopic Field Dislocation Mechanics (hereafter abbreviated “MFDM”), was developed by Acharya and Roy (2006). Later, Acharya (2010, 2011) developed more thoroughly the FDM and MFDM theories with the thermodynamic aspects. It was seen to be successful in modeling different problems in plasticity at mesoscopic scale and small strains. It was able to describe different microstructural size effects in single crystalline materials or multicrystalline thin films (Acharya and Roy, 2006; Acharya et al., 2006; Roy and Acharya, 2006; Roy et al., 2007; Puri et al., 2009, 2010, 2011; Puri and Roy, 2012). A comparison of the MFDM theory (full and reduced versions) with lower order strain gradient crystal plasticity models were provided in Roy et al. (2007). In this mesoscale theory, the constitutive equations for strain-hardening models, slip rule for SSDs and velocity of GNDs need to be specified phenomenologically, but it allows performing simulations at reasonable strain rates compared to DDD. Transport of GND density was numerically solved by the Finite Element (FE) method for three-dimensional setting (Roy and Acharya, 2005, 2006; Varadhan et al., 2006) or one-dimensional finite difference scheme (Das et al., 2016). In this framework, Acharya (2007) also formulated a jump condition at a material interface, like a phase or a grain boundary, adapted for both the FDM and the MFDM theories based on the conservation of Burgers vector. For a fixed interface, this condition states a tangential continuity of the plastic distortion rate, which has consequences on interfacial GND content and slip. Such continuity constraints at material interfaces were used to study the mechanical responses of bicrystals (Puri et al., 2010), multicrystalline thin films (Puri et al., 2011) and metal matrix composites (Richeton et al., 2011; Taupin et al., 2012). More recently, Arora and Acharya (2020a,b) implemented a FE formulation of the finite deformation MFDM theory to study the formation of deformation patterns in single crystals under simple shear in both J2-Von Mises plastic flow theory and crystal plasticity frameworks. Then, the finite deformation MFDM theory was successfully applied to micropillar confined thin film plasticity in Arora et al. (2022), and to describe kink band formation in Cu-Nb nano-metallic laminates (NMLs) in Arora et al. (2023). In their constitutive formulation, a dislocation core energy

density with material length scale *l* was used, and an augmented dissipative plastic strain rate was developed in both J2-Von Mises plastic flow theory and crystal plasticity frameworks. A GND dependent backstress tensor together with a modification of the Peach-Koehler (PK) driving force for GND motion was proposed (Acharya, 2010, 2011; Das et al., 2016; Arora and Acharya, 2020a,b). Finally, it is noteworthy that other continuum dislocation-based models were developed these last decades (Hochrainer et al., 2007; Sandfeld et al., 2010; Hochrainer, 2015; Xia and El-Azab, 2015; Wulfinghoff and Böhlke, 2015; Sudmanns et al., 2019; Schulz et al., 2019; Starkey et al., 2020; Vivekanandan et al., 2021).

Due to computational efficiency, an alternative numerical implementation of the MFDM theory with the FE method is to use the fast Fourier transforms (FFT) algorithm. FFT-based methods were originally developed and applied to composite materials (Moulinec and Suquet, 1994, 1998), and later to viscoplastic crystal plasticity (CP-VPFFT) (Lebensohn, 2001) and elastoviscoplastic crystal plasticity (CP-EVPFFT) calculations (Lebensohn et al., 2012), for applications to polycrystal plasticity. These original CP-EVPFFT-based implementations showed the feasibility of efficiently solving the micromechanical behavior of complex polycrystalline unit cells. The EVPFFT-based methods were also recently extended to consider constitutive behaviors at small strains, based on HO-SGP with the SG-EVPFFT formulation of Lebensohn and Needleman (2016), based on combined LO-SGP and HO-SGCP by Christodoulou et al. (2022), and based on MFDM with the MFDM-EVPFFT formulation by Djaka et al. (2020), Berbenni et al. (2020), Berbenni and Lebensohn (2021). Recently, a non-local large-strain SG-LS-EVPFFT model was developed by Zecevic et al. (2023), and it was applied to successfully simulate kink band formation in Cu-Nb NMLs, which was not possible with the small strain SG-EVPFFT model.

In this contribution, the small strain EVPFFT formulation coupled with MFDM (MFDM-EVPFFT) as previously developed in Djaka et al. (2020), Berbenni et al. (2020) and Berbenni and Lebensohn (2021), is rethought and extended to include a *defect energy* density, that depends on GND densities and an associated material length scale *l*. Therefore, a backstress formulation is introduced in the three-dimensional setting of MFDM, in the same spirit as Acharya (2010) and Das et al. (2016) for a simplified one-dimensional model. The discrete Fourier transform (DFT) method together with a centered finite difference (FD) scheme (Berbenni et al., 2014; Lebensohn and Needleman, 2016) is used to calculate the evolution of the backstress tensor. The paper is organized as follows. In Section 2, the thermodynamic framework related to the Mesoscale Field Dislocation Mechanics (MFDM) with defect energy is introduced after briefly recalling the field equations of the MFDM theory in small strains. Then, in Section 3, the specific numerical implementation of the elasto-viscoplastic FFT-based method for MFDM (MFDM-EVPFFT) is presented. Section 4 reports the different numerical results for two-phase laminate composites with plastic single crystal channels and elastic precipitates for monotonic and reversible shear loadings, and for a polycrystal during tension-compression loading. The comparisons of the results between the MFDM-EVPFFT model with different hardening types (“energetic”, “dissipative”, “combined”) and the CP-EVPFFT model are discussed in the same section. Section 5 concludes and draws perspectives to this work.

## Notation

A bold symbol denotes a tensor or a vector. The symmetric part of tensor A is denoted $A^{sym}$. Its skew-symmetric part is $A^{skew}$ and its transpose is denoted by $A^T$. The tensor $A \cdot B$, with rectangular Cartesian components $A_{ik}B_{kj}$, results from the dot product of tensors A and B, and $A \otimes B$ is their tensorial product, with components $A_{ij}B_{kl}$. The vector $A \cdot V$, with rectangular Cartesian components $A_{ij}V_j$, results from the dot product of tensor A and vector V. The symbol “:” represents the trace inner product of the two second order tensors $A: B = A_{ij}B_{ij}$,
in rectangular Cartesian components, or the product of a higher rank with a second rank tensor, e.g., $\mathbf{A} : \mathbf{B} = A_{ijkl}B_{kl}$. The cross product of a second rank tensor $\mathbf{A}$ and a vector $\mathbf{V}$, the div and curl operations for second rank tensors are defined row by row, in analogy with the vectorial case. For any base vector $\mathbf{e}_i$ of the reference frame:

$$(\mathbf{A} \times \mathbf{V})^t \cdot \mathbf{e}_i = (\mathbf{A}^t \cdot \mathbf{e}_i) \times \mathbf{V} \qquad (1)$$

$$(\mathrm{div} \, \mathbf{A})^t \cdot \mathbf{e}_i = \mathrm{div}(\mathbf{A}^t \cdot \mathbf{e}_i) \qquad (2)$$

$$(\mathrm{curl} \, \mathbf{A})^t \cdot \mathbf{e}_i = \mathrm{curl}(\mathbf{A}^t \cdot \mathbf{e}_i) \qquad (3)$$

In rectangular Cartesian components:

$$(\mathbf{A} \times \mathbf{V})_{ij} = e_{jkl} A_{ik} V_l \qquad (4)$$

$$(\mathrm{div} \, \mathbf{A})_i = A_{ij,j} \qquad (5)$$

$$(\mathrm{curl} \, \mathbf{A})_{ij} = e_{jkl} A_{il,k} = -(\mathrm{grad} \, \mathbf{A} : \mathbf{X})_{ij} \qquad (6)$$

where $e_{jkl}$ is a component of the third rank alternating Levi-Civita tensor $\mathbf{X}$ and the spatial derivative with respect to a Cartesian coordinate is indicated by a comma followed by the component index. The second rank $\delta$ tensor denotes the Kronecker delta operator. The notation $\hat{\mathbf{A}}(\xi)$ will be used for the Fourier transform of $\mathbf{A}(\mathbf{x})$.

In Fourier space, $\xi$ is the Fourier vector of magnitude $\xi = \sqrt{\xi^2 \cdot \xi^2}$ and components $\xi_i$. The imaginary unit is $i = \sqrt{-1}$.

## 2. Mesoscale field dislocation mechanics with defect energy

### 2.1. Mesoscale field equations (reduced form)

A small strain continuum mechanics framework is assumed, where the displacement vector field $\mathbf{u}$ is defined continuously at any point of an elasto-viscoplastic body. The total curl-free distortion tensor $\mathbf{U}$ is the gradient of the displacement vector $\mathbf{U} = \mathrm{grad} \, \mathbf{u}$. At the mesoscopic scale (Acharya and Roy, 2006), both the average plastic distortion $\mathbf{U}^p$, which results from dislocation motion, and the average elastic (or lattice) distortion $\mathbf{U}^e$ are incompatible fields. In classical mechanics, mechanical fields like the Cauchy stress field $\sigma$ and the displacement field $\mathbf{u}$, are derived from the balance of linear momentum together with the balance of angular momentum (see Eq. (7)). In addition, periodic boundary conditions over a unit cell of period denoted $H$ for the FFT-based method are assumed. Field equations are the following:

$$\mathrm{div} \, \mathbf{\sigma} = 0$$

$$\mathbf{\sigma} = \mathbf{\sigma}^T$$

$$\mathbf{U} = \mathrm{grad} \, \mathbf{u} = \mathbf{U}^e + \mathbf{U}^p \qquad (7)$$

$$\mathbf{\sigma} = \mathbf{C} : \varepsilon^e = \mathbf{C} : (\varepsilon - \varepsilon^p)$$

$$\mathbf{u} - \langle \varepsilon \rangle \times \text{periodic, } \mathbf{\sigma} \cdot \mathbf{n} \text{ anti-periodic}$$

where $\mathbf{C}$ is the classic fourth order elastic stiffness tensor satisfying minor and major symmetries, $\varepsilon^e = (\mathbf{U}^e)^{sym}$ is the elastic strain tensor, $\varepsilon^p = (\mathbf{U}^p)^{sym}$ is the plastic strain tensor, and $\varepsilon = (\mathbf{U})^{sym}$ is the total (small) compatible strain tensor. $\langle A \rangle$ denotes the spatial average of $A$ over $H$.

In the Mesoscale Field Dislocation Mechanics (MFDM) and considering here a reduced version (Roy et al., 2007; Acharya and Roy, 2006), the time-evolution of the plastic strain distortion is calculated using the averaging procedure described in Acharya and Roy (2006):

$$\mathbf{U}^p = \mathbf{\alpha} \times \mathbf{v} + \mathbf{L}^p \qquad (8)$$

This relationship describes the Orowan-type law at the mesoscale for Geometrically Necessary Dislocations (GND) through the term $\mathbf{\alpha} \times \mathbf{v}$ (where $\mathbf{\alpha}$ is the average Nye-Kröner dislocation density tensor (Nye, 1953; Kröner, 1958) and $\mathbf{v}$ is the average GND velocity vector), and for Statistically Stored Dislocations (SSD) through the $\mathbf{L}^p$ term. Only the latter is generally prescribed in conventional crystal plasticity, see e.g. Asaro (1983) and Roters et al. (2010).

Then, the space-time evolution of the average dislocation density tensor $\mathbf{\alpha}$ is obtained from the conservation of dislocation flux across

an elementary surface (Mura, 1963, 1987; Acharya, 2001, 2003) and adapted here at the mesoscale:

$$\mathbf{\alpha} = -\mathrm{curl} \, \mathbf{U}^p \qquad (9)$$

This equation will be solved in Section 3 using a spectral approach (Djaka et al., 2015) that was developed to solve in a fast, stable and accurate fashion, the hyperbolic transport equation governing the spatio-temporal evolution of the dislocation density tensor in the MFDM. In addition, the evolution equations for MFDM impose a jump condition on plastic distortion rate across a fixed (i.e. non moving) material interface (Acharya, 2007):

$$[\mathbf{U}^p] \times \mathbf{n} = \mathbf{0} \qquad (10)$$

where $\mathbf{n}$ is the unit outward normal vector to the interface and $[\mathbf{U}^p]$ is the jump of $\mathbf{U}^p$ between both media across the interface. This equation was numerically solved in the FE framework for both impenetrable and penetrable interfaces in the case of thin films by Puri et al. (2011). Here, we will consider only the impenetrable interface condition in the case of two-phase composites (see Section 4.1 for details).

Constitutive specifications on the GND velocity $\mathbf{v}$, and the plastic distortion rate $\mathbf{L}^p$ are required for the theory. This is now discussed using thermodynamic considerations with defect energy, and a crystal plasticity framework is described with both energetic and dissipative hardening.

### 2.2. Thermodynamics-based framework with defect energy

Here, we follow a procedure initiated by Coleman and Noll (1963) and Coleman and Gurtin (1967) for local constitutive equations, which was more recently used by Acharya (2010) for small strain microscopic FDM and by Arora and Acharya (2020a) for finite deformation MFDM. The Helmholtz free energy density function is introduced as follows:

$$\phi = \hat{\phi}(\varepsilon^e) + \phi^d(\mathbf{\alpha}) \qquad (11)$$

where $\hat{\phi}$ is the classic strain energy density, and $\phi^d$ is an augmented phenomenological energy density depending explicitly on the dislocation density tensor $\mathbf{\alpha}$. Here, we assume linear elasticity for $\hat{\phi}$ and a simple quadratic form for $\phi^d$, which will be called defect energy following the same terminology as Gurtin (2002), although the present model does not belong to HO-SGP theories. Therefore, we conjecture that:

$$\phi = \frac{1}{2} \varepsilon^e : \mathbf{C} : \varepsilon^e + \frac{1}{2} \mathbf{A} \mathbf{\alpha} : \mathbf{\alpha} \qquad (12)$$

where $\mathbf{A}$ is a constant material length dependent parameter, which will be discussed later. It is noteworthy that the free energy density depends only on GND densities (and not on SSD densities) in this formulation. It is essentially because GND density defines an observed state from the dislocation density tensor $\mathbf{\alpha}$ which can be measured for example by Electron Back-Scattered Diffraction (EBSD) (Pantleon, 2008; Allain-Bonasso et al., 2012; Wallis et al., 2016; Genée et al., 2021). It follows the statement of the continuum theory of dislocations where GND density is a state variable (Kröner, 1958, 1981, 1996). Therefore, the GND densities present in the free energy density enhances the effect of local incompatible lattice distortion on the constitutive behavior of the material. This is at the origin of new hardening contained in the backstress tensor (which is not due to SSD). In the present MFDM theory, the SSD are only considered in the evolution of plastic distortion through $\mathbf{L}^p$ and in the dissipative (isotropic) hardening law (see Section 2.3).

The present formulation follows the same guidelines as the MFDM theory reported by Acharya (2010, 2011), Das et al. (2016), Arora and Acharya (2020a,b) and Arora et al. (2022), but is restricted to small deformation. Indeed, Arora and Acharya (2020a,b) introduced $\phi^d$ and called it dislocation core energy for finite deformation crystal plasticity and J2 flow theories at the mesoscale. Here, a GND dependent
backstress tensor together with a modification of the Peach-Koehler (PK) driving force for GND motion will be derived as in the equations of Acharya (2010, 2011), Arora and Acharya (2020a,b) and Das et al. (2016).

The constitutive equations for $L^p$ and $v$ are based on the expression of the dissipation rate $D$ with the requirement to be positive, i.e. $D \ge 0$. For a body $V$, neglecting body forces, $D$ reads:

$$D = P_{\text{ext}} - \Phi = \int_{\partial V} \sigma \cdot n \cdot \dot{U} \, dS - \int_{V} \dot{\Phi} \, dV \quad (13)$$

where $n$ is the outward unit normal vector to the boundary surface $\partial V$. In Eq. (13), the time-derivative of the Helmholtz free energy density $\dot{\Phi}$ is given by:

$$\dot{\Phi} = \frac{\partial \Phi}{\partial \varepsilon^e} : \dot{\varepsilon}^e + \frac{\partial \Phi}{\partial \alpha} : \dot{\alpha} = \sigma : \dot{\varepsilon}^e + A\alpha : \dot{\alpha} \quad (14)$$

Then, using the divergence theorem together with Eq. (7) and the symmetries of $C$, the dissipation rate $D$ yields:

$$D = \int_V \sigma : \dot{U}^p dV - A \int_V \alpha : \dot{\alpha} \ge 0 \quad (15)$$

Using simultaneously both equations Eqs. (8) and (9) in Eq. (15), the dissipation rate $D$ becomes:

$$D = \int_V \sigma : L^p dV + \int_V \sigma : (\alpha \times v) dV + A \int_V \alpha : \text{curl } L^p dV + A \int_V \alpha : \text{curl } (\alpha \times v) dV \quad (16)$$

The last expression can be simplified by applying integration by parts and using the third order Levi-Civita tensor $X$ (see Section 1), such that finally Eq. (16) yields:

$$D = \int_V (\sigma - B) : L^p dV + \int_V F^a \cdot v dV - A \int_{\partial V} \alpha : (\dot{U}^p \times n) dS \quad (17)$$

where $B$ is a non-symmetric backstress tensor defined as:

$$B = -A \text{curl } \alpha \quad (18)$$

and, $F^a$ is the Peach-Koehler (PK) driving force for GND motion defined as:

$$F^a = X : (\sigma - B)^T \cdot \alpha \quad (19)$$

or:

$$F^a = X : \sigma \cdot \alpha + AX : (\text{curl } \alpha)^T \cdot \alpha \quad (20)$$

Hence, in comparison with our recent contributions on MFDM-EVPFFT (Djaka et al., 2020; Berbenni et al., 2020; Berbenni and Lebensohn, 2021), a second (length scale dependent) term is present in Eq. (20) in the present MFDM formulation. Following Acharya (2010), we assume that the last term in Eq. (17) (boundary term) is negligible compared to volume terms and will not be considered in the present EVPFFT-based model. Therefore, from Eq. (17), it is seen that $\sigma - B$ is the thermodynamic conjugate of $L^p$ and $F^a$ is the thermodynamic conjugate of $v$. These results exactly match the small deformation specialization of Arora and Acharya (2020a) in terms of free energy. They are also consistent with the thermodynamic derivations of Acharya (2011) in terms of specific internal energy and entropy, with an equivalence in terms of free energy and absolute temperature, see Appendix A in Acharya (2011).

## 2.3. Crystal plasticity with energetic (kinematic) and dissipative (isotropic) hardening specifications

Considering a crystal plasticity framework with dislocation glide (no dislocation climb), the plastic distortion rate tensor $L^p$ due to crystallographic slip reads:

$$L^p = \sum_{s=1}^{N} \dot{\gamma}^s b^s \otimes n^s, \quad (21)$$

where $s$ denotes a slip system, $N$ is the total number of slip systems ($N = 12$ for FCC crystals), $\dot{\gamma}^s$ is the slip rate, $b^s$ denotes the unit slip direction vector, and $n^s$ is the slip plane unit normal. Plastic strain incompressibility holds such that both tensors $L^p$ and $\alpha \times v$ are traceless. From Eq. (21), it means that the crystallographic orientation tensor $m^s = b^s \otimes n^s$ is also traceless. Defining the dilatant direction vector $d = X : \alpha$, then plastic incompressibility imposes that $d \perp v$.

Using Eq. (21) together with Eq. (17) (neglecting the last boundary term in Eq. (17)) yields:

$$D = \int_V \sum_{s=1}^{N} \tau_{eff}^s \dot{\gamma}^s dV + \int_V F^a \cdot v \, dV \ge 0 \quad (22)$$

where $\tau_{eff}^s$ is an effective resolved shear stress on $s$ given by:

$$\tau_{eff}^s = m^s : (\sigma - B) = \tau^s - \tau_{kin}^s \quad (23)$$

where $\tau^s = m^s : \sigma$ and $\tau_{kin}^s = m^s : B$. The term $\tau_{kin}^s$ is a resolved shear stress due to the backstress tensor, which is at the origin of intracrystalline kinematical-type hardening. In the rest of the paper, we will denominate it as energetic hardening, since it stems from the defect energy formulation and it is thus strictly energetic.

Then, the slip rate on each slip system $s$ is taken as a power function of $\tau_{eff}^s$, which fulfills the thermodynamics requirement: $\tau_{eff}^s \dot{\gamma}^s \ge 0$:

$$\dot{\gamma}^s = \dot{\gamma}^0 \left( \frac{|\tau_{eff}^s|}{\tau_c} \right)^n \operatorname{sgn}(\tau_{eff}^s) \quad (24)$$

where $\dot{\gamma}^0$ is a reference slip rate, $n$ is a power law exponent, and $\tau_c$ is the material's shear strength considered identical for all slip systems and evolving with strain (as it will be reported later).

Furthermore, to fulfill the thermodynamics requirement $F^a \cdot v \ge 0$, the average GND velocity $v$ at the mesoscale is constitutively prescribed as follows:

$$v = \frac{g}{|g|} \bar{v} \quad \text{with} \quad \bar{v} \ge 0, \quad (25)$$

where $\bar{v}$ is the magnitude of $v$, and $g$ is the glide force parallel to $v$ defined by:

$$g = F^a - \left( F^a \cdot \frac{d}{|d|} \right) \frac{d}{|d|}, \quad (26)$$

with $F^a$ given by Eq. (19) including the new backstress tensor $B$. Here, to compute $g$ for the present model, we apply the same treatment as in Djaka et al. (2020) and Berbenni et al. (2020) (see details therein). Then, using Eqs. (20) and (22), the component form of $g$ is now given by:

$$g_r = e_{rik} S_{ij} \alpha_{jk} - e_{rik} \alpha_{ik} \frac{S_{mn} \alpha_{np} (\alpha_{mp} - \alpha_{pm})}{|d|^2} + A e_{rik} (\text{curl } \alpha)^T_{mn} \alpha_{np} (\alpha_{mp} - \alpha_{pm}) \frac{(\text{curl } \alpha)^T_{mn} \alpha_{np} (\alpha_{mp} - \alpha_{pm})}{|d|^2} \quad (27)$$

Thus, this expression for $g$ extends our previous works (Djaka et al., 2020; Berbenni et al., 2020; Berbenni and Lebensohn, 2021) by considering a second term which is strictly energetic and derived from the defect energy.

The constitutive specification for the average GND velocity magnitude $\bar{v}$ introduced in Eq. (25) follows the mechanistic formulation given by Puri et al. (2010, 2011), which is based on a combination of Orowan law and Bailey-Hirsch strengthening relationship:

$$\bar{v} = \frac{\eta^2 b}{N} \left( \frac{\mu}{\tau_c} \right)^2 \sum_{s=1}^{N} |\dot{\gamma}^s| \quad (28)$$

where $\eta$ is a material constant close to $1/3$, $b$ is the magnitude of the Burgers vector, and $\mu$ is the isotropic elastic shear modulus of the material.

The evolution law for the shear strength $\tau_c$ is based on the isotropic strain-hardening model used in the earlier MFDM-EVPFFT model (Djaka et al., 2020; Berbenni et al., 2020; Berbenni and Lebensohn, 2021), which was an extension of the pioneering works of Mecking
and Kocks (1981) and Acharya and Beaudoin (2000):

$$ \dot{t}_c = \theta_0 \frac{\tau_s - \tau_c}{\tau_s - \tau_0} \dot{\Gamma} + k_0 \frac{\eta^2 \mu^2 b}{2(\tau_c - \tau_0)} \left( \sum_{s=1}^{N} |\boldsymbol{\alpha} \cdot \mathbf{n}^s| |\dot{\gamma}^s| + \sum_{s=1}^{N} |\boldsymbol{\alpha} \cdot \mathbf{n}^s| |\boldsymbol{\alpha} \times \mathbf{v}| \right) \quad (29) $$

where $\tau_0$ is the initial strength due to lattice friction for FCC materials, $\tau_s$ is the saturation stress, $\theta_0$ is the stage II hardening rate, $\dot{\Gamma}$ is the cumulative slip rate on all slip systems:

$$ \dot{\Gamma} = |\boldsymbol{\alpha} \times \mathbf{v}| + \sum_{s=1}^{N} |\dot{\gamma}^s| \quad (30) $$

Furthermore, the material parameter $k_0$ is related to a geometric mean free path $L_g$ due to GND forest on slip system s (Acharya and Beaudoin, 2000):

$$ L_g = \frac{1}{k_0 |\boldsymbol{\alpha} \cdot \mathbf{n}^s|} \quad (31) $$

Eq. (29) represents isotropic hardening, and it will be called in the rest of the paper dissipative hardening following the same terminology as Gurtin (2002) and Fleck and Willis (2015). The hardening rule used in MFDM (Acharya and Roy, 2006) is strictly dissipative without defect energy formulation (as opposed to backstress and energetic hardening).

Let us finally note that both energetic hardening and dissipative hardening both lead to size dependent hardening. In the case of energetic hardening, the length scale dependent modulus A will be considered as $\frac{1}{2}\tau_0 l^2$, where l is a material length parameter (which has a dimension of length) and $\tau_0$ is a material strength parameter (which has a dimension of stress) (Gurtin, 2002; Bittencourt et al., 2003; Lebensohn and Needleman, 2016; Christodoulou et al., 2022). For dissipative hardening, the material length scale is given by $L_g$ in Eqs. (29) and (31), which is calibrated by material parameter $k_0$. The effect of $k_0$ on cyclic plasticity in two-phase laminate composites was recently studied by Berbenni and Lebensohn (2021) using MFDM-EVPFFT. In the next numerical simulations reported in Section 4, $k_0$ will be fixed to 20 for FCC metals following previous contributions (Acharya and Beaudoin, 2000; Puri et al., 2011; Berbenni et al., 2020; Arora and Acharya, 2020a; Arora et al., 2022).

# 3. Elasto-viscoplastic FFT-based implementation (MFDM-EVPFFT)

## 3.1. EVPFFT method with augmented Lagrangian (AL) scheme

An elasto-viscoplastic crystal plasticity formulation is adopted starting from the small strain EVPFFT formulation earlier developed by Lebensohn et al. (2012). This small strain formulation was further extended to SG-EVPPFT by Lebensohn and Needleman (2016) and Christodoulou et al. (2022) and to MFDM-EVPFFT by Djaka et al. (2020), Berbenni et al. (2020) and Berbenni and Lebensohn (2021). In this formulation, time discretization is handled using a backward Euler implicit scheme. Therefore, the Cauchy stress tensor $\boldsymbol{\sigma}$ at $t + \Delta t$ is given by:

$$ \boldsymbol{\sigma}^{t+\Delta t} = \mathbf{C} : \boldsymbol{\varepsilon}^{e,t+\Delta t} = \mathbf{C} : (\boldsymbol{\varepsilon}^{t+\Delta t} - \boldsymbol{\varepsilon}^{p,t} - \dot{\boldsymbol{\varepsilon}}^{p,t+\Delta t}(\boldsymbol{\sigma}^{t+\Delta t})\Delta t), \quad (32) $$

Now, the supra-indices $t + \Delta t$ are omitted, and the fields corresponding to the previous time step $t$ will be explicitly indicated.

Introducing the Green's function method (Mura, 1987; Moulinec and Suquet, 1994, 1998), the strain tensor $\boldsymbol{\varepsilon}$ defined in Eqs. (7) and (32) is numerically solved iteratively using the fast Fourier transform (FFT)-based Augmented Lagrangian (AL) scheme introduced by Michel et al. (2001) and adapted by Lebensohn et al. (2012) for crystal plasticity-based EVPFFT.

Let $\hat{\boldsymbol{\varepsilon}}(\xi)$ and $\hat{\boldsymbol{\Gamma}}^0(\xi)$ be, respectively, the Fourier transform of $\boldsymbol{\varepsilon}(\mathbf{x})$ and of the strain-based Green operator $\boldsymbol{\Gamma}^0(\mathbf{x})$. As in Djaka et al. (2020)

and Berbenni et al. (2020), the numerical computation of the components of the strain-based Green operator in Fourier space $\hat{\boldsymbol{\Gamma}}^0(\xi)$ as part of the solution of the Lippmann-Schwinger equation is performed using a centered finite difference scheme on a rotated grid as introduced by Willot (2015). The homogeneous reference elastic stiffness tensor $\mathbf{C}^0$ associated with $\hat{\boldsymbol{\Gamma}}^0(\xi)$ is taken as the average of $\mathbf{C}$ over the unit cell.

The iterative procedure performed at $t + \Delta t$ is briefly described now. Let us assume that $\boldsymbol{\lambda}^{(n)}$ and $\boldsymbol{\varepsilon}^{(n)}$ are, respectively, the auxiliary guess stress and strain fields at iteration (n). Then, their Fourier transforms are computed with the FFT algorithm: $\hat{\boldsymbol{\varepsilon}}^{(n)}(\xi)$ and $\hat{\boldsymbol{\lambda}}^{(n)}$. From the Lippmann-Schwinger equation in the Fourier space and using the property of the strain-based Green operator (Michel et al., 2001), the auxiliary strain field $\boldsymbol{\varepsilon}^{(n+1)}(\mathbf{x})$ is solved at (n+1):

$$ \begin{aligned} \hat{\boldsymbol{\varepsilon}}^{(n+1)}(\xi) &= \hat{\boldsymbol{\varepsilon}}^{(n)}(\xi) - \hat{\boldsymbol{\Gamma}}^0(\xi) : \hat{\boldsymbol{\lambda}}^{(n)}(\xi) \quad \forall \xi \neq \mathbf{0} \\ \hat{\boldsymbol{\varepsilon}}^{(n+1)}(\mathbf{0}) &= \langle \boldsymbol{\varepsilon}^{(n)} \rangle \\ \boldsymbol{\varepsilon}^{(n+1)}(\mathbf{x}) &= \text{FFT}^{-1}(\hat{\boldsymbol{\varepsilon}}^{(n+1)}(\xi)) \end{aligned} \quad (33) $$

where FFT$^{-1}$ is the inverse Fourier transform. Once $\boldsymbol{\varepsilon}^{(n+1)}$ is obtained, the nullification of the residual $\mathbf{R}$ is performed, which depends on the stress and strain tensors $\boldsymbol{\sigma}^{(n+1)}$ and $\boldsymbol{\varepsilon}^{(n+1)}$:

$$ \mathbf{R}(\boldsymbol{\sigma}^{(n+1)}) = \boldsymbol{\sigma}^{(n+1)} + \mathbf{C}^0 : \boldsymbol{\varepsilon}^{(n+1)}(\boldsymbol{\sigma}^{(n+1)}) - \boldsymbol{\lambda}^{(n)} - \mathbf{C}^0 : \boldsymbol{\varepsilon}^{(n+1)} = \mathbf{0} \quad (34) $$

This equation is non-linear and is solved with a Newton-Raphson type scheme. The (p + 1)-guess for the stress field $\boldsymbol{\sigma}^{(n+1)}$ is given by:

$$ \boldsymbol{\sigma}^{(n+1,p+1)} = \boldsymbol{\sigma}^{(n+1,p)} - \left( \left( \frac{\partial \mathbf{R}}{\partial \boldsymbol{\sigma}} \right)_{\boldsymbol{\sigma}^{(n+1,p)}} \right)^{-1} : \mathbf{R} \left( \boldsymbol{\sigma}^{(n+1,p)} \right) \quad (35) $$

where the Jacobian $\frac{\partial \mathbf{R}}{\partial \boldsymbol{\sigma}}$ must be calculated:

$$ \left( \frac{\partial \mathbf{R}}{\partial \boldsymbol{\sigma}} \right)_{\boldsymbol{\sigma}^{(n+1,p)}} = \delta \otimes \delta + \mathbf{C}^0 : \mathbf{C}^{-1} + \Delta t \mathbf{C}^0 : \left( \frac{\partial \boldsymbol{\varepsilon}^p}{\partial \boldsymbol{\sigma}} \right)_{\boldsymbol{\sigma}^{(n+1,p)}} \quad (36) $$

Once the convergence is achieved on $\boldsymbol{\sigma}^{(n+1)}$ and $\boldsymbol{\varepsilon}^{(n+1)}$, the new guess for the auxiliary stress field $\boldsymbol{\lambda}^{(n+1)}$ is given using the Uzawa's descent algorithm:

$$ \boldsymbol{\lambda}^{(n+1)} = \boldsymbol{\lambda}^{(n)} + \mathbf{C}^0 : (\boldsymbol{\varepsilon}^{(n+1)} - \boldsymbol{\varepsilon}^{(n+1)}) \quad (37) $$

and the algorithm is stopped when the normalized average differences between the stress fields $\boldsymbol{\sigma}$ and $\boldsymbol{\lambda}$, and the strain fields $\boldsymbol{\varepsilon}$ and $\boldsymbol{\varepsilon}$, are smaller than a given threshold error (typically $10^{-5}$), which fulfill both stress equilibrium and strain compatibility up to sufficient accuracy. In the algorithm described above, an overall macroscopic strain $\mathbf{E} = \langle \boldsymbol{\varepsilon}^{(n)} \rangle$ is applied to the periodic unit cell: $\mathbf{E} = \mathbf{E}^t + \hat{\boldsymbol{\mathcal{E}}} \Delta t$. For mixed boundary conditions with applied strain rate $\hat{\boldsymbol{\mathcal{E}}}$ and applied stress $\boldsymbol{\Sigma}$, the (n + 1)-guess of the macroscopic strain $\mathbf{E}^{(n+1)}$ was detailed in Michel et al. (2001).

In this iterative procedure, the direct and the inverse Fourier transforms are computed here by using fast Fourier transform (FFT) algorithm. The spatial periods of the unit cell are $H_1, H_2$ and $H_3$ in the $x_1, x_2$ and $x_3$ directions, respectively, and discretized by a regular rectangular grid with $N_1 \times N_2 \times N_3$ voxels with position vector $\mathbf{x} = (i_1 \delta_1, i_2 \delta_2, i_3 \delta_3)$, where $i_1 = 0 \rightarrow N_1 - 1, i_2 = 0 \rightarrow N_2 - 1, i_3 = 0 \rightarrow N_3 - 1$ and $\delta_1, \delta_2, \delta_3$ are the voxel sizes in the $x_1, x_2$ and $x_3$ directions (here $\delta_1 = \delta_2 = \delta_3 = \delta$). The computational grid is constituted of a total of $N_{tot} = N_1 \times N_2 \times N_3$ voxels. Then, Fourier vectors $\xi$ are discretized by $\xi_j = 2\pi m_j / H_j$, where $m_j$ ($j = 1 \rightarrow 3$) are given by Moulinec and Suquet (1998):

$$ m_j = \left( \left( -\frac{N_j}{2} + 1 \right), \left( -\frac{N_j}{2} + 2 \right), \dots, -1, 0, 1, \dots, \left( \frac{N_j}{2} - 1 \right), \left( \frac{N_j}{2} \right) \right) \quad (38) $$

if $N_j$ is even, and

$$ m_j = \left( \left( -\frac{N_j - 1}{2} \right), \dots, -1, 0, 1, \dots, \left( \frac{N_j - 1}{2} \right) \right) \quad (39) $$

if $N_j$ is odd.
## 3.2. Numerical details of the present MFDM-EVPFFT with backstress implementation

First of all, the expression of $\partial\dot{\epsilon}^p/\partial\sigma$ in the Jacobian (Eq. (36)) is modified considering the constitutive equations of the MFDM theory (see Section 2):

$$ \left(\frac{\partial \dot{\epsilon}^p}{\partial \sigma}\right)_{\sigma^{(n+1,p)}} = \left(\frac{\partial L^p}{\partial \sigma}\right)_{\sigma^{(n+1,p)}}^{\text{sym}} + \left(\frac{\partial (\alpha \times \mathbf{v})}{\partial \sigma}\right)_{\sigma^{(n+1,p)}}^{\text{sym}} \quad (40) $$

In component form, $(\partial L^p/\partial\sigma)_{\sigma^{(n+1,p)}}^{\text{sym}}$ is calculated using Eqs. (21), (23) and (24):

$$ \left(\frac{\partial L_{kl}^p}{\partial \sigma_{mn}}\right)_{\sigma^{(n+1,p)}} \approx n\gamma^0 \sum_{s=1}^{N} m_{kl}^s m_{mn}^s \frac{(|\tau_{eff}^s|)^{n-1}}{(\tau_c)^n} \quad (41) $$

The calculation of $(\partial(\alpha \times \mathbf{v})/\partial\sigma)_{\sigma^{(n+1,p)}}$ is computed as follows in component form:

$$ \left(\frac{\partial(\alpha \times \mathbf{v})_{kl}}{\partial \sigma_{mn}}\right)_{\sigma^{(n+1,p)}} = e_{lqr} \alpha_{kq} \left( \frac{\partial (g_r/|\mathbf{g}|)}{\partial \sigma_{mn}} \bar{v} + \frac{g_r}{|\mathbf{g}|} \frac{\partial \bar{v}}{\partial \sigma_{mn}} \right)_{\sigma^{(n+1,p)}} \quad (42) $$

Using the new expression of $\mathbf{g}$ (Eq. (27)), we obtain in component form:

$$ \frac{\partial (g_r / |\mathbf{g}|)}{\partial \sigma_{mn}} = \left( \frac{\delta_{rs} |\mathbf{g}|^2 - g_r g_s}{|\mathbf{g}|^3} \right) \left( e_{oks} \alpha_{qk} - e_{iks} \alpha_{ik} \frac{\alpha_{qp} (\alpha_{op} - \alpha_{po})}{\alpha_{ij} (\alpha_{ij} - \alpha_{jl})} \right) \left( \delta_{om} \delta_{qn} - \frac{1}{3} \delta_{mn} \delta_{oq} \right) \quad (43) $$

and using Eq. (28), we obtain:

$$ \left(\frac{\partial \bar{v}}{\partial \sigma_{mn}}\right)_{\sigma^{(n+1,p)}} \approx n\gamma^0 \frac{\eta^2 b}{N} \left(\frac{\mu}{\tau_c}\right)^2 \sum_{s=1}^{N} m_{mn}^s \frac{(|\tau_{eff}^s|)^{n-1}}{(\tau_c)^n} \quad (44) $$

In Eqs. (41) and (44), the approximation lies in the fact that the derivatives $\partial\tau_c/\partial\sigma$ and $\partial\mathbf{m}^2/\partial\sigma$ are neglected. In Eq. (43), the partial derivatives are calculated at fixed $\alpha(\mathbf{x})$.

To solve Eq. (9), let $\hat{\alpha}(\xi)$ be the Fourier transform of $\alpha(\mathbf{x})$. It is explicitly time-integrated using a spectral solver (Gottlieb and Hesthaven, 2001) with an exponential second order spectral low-pass filter, denoted $\kappa(\eta)$, following Djaka et al. (2015, 2020), Berbenni et al. (2020) and Berbenni and Lebensohn (2021):

$$ \hat{\alpha}_{ij}^{t+\Delta t} = \kappa(\eta) \left( \hat{\alpha}_{ij}^t - \Delta t \mathbf{i} \xi_k \left( \widehat{(\alpha_{ij} v_k)}^t - \widehat{(\alpha_{ik} v_j)}^t \right) \right) - \Delta t \mathbf{i} \xi_k e_{ijkl} \widehat{(L_{il}^p)}^t \quad (45) $$

The exponential filter, which is set up to eliminate high frequencies responsible for spurious oscillations (Gottlieb and Hesthaven, 2001; Djaka et al., 2015), is defined as function of discrete frequencies $\eta_j = m_j/N_j$:

$$ \kappa \left( \frac{m_1}{N_1}, \frac{m_2}{N_2}, \frac{m_3}{N_3} \right) = \exp \left( -\beta \left( \left( \frac{m_1}{N_1} \right)^{2p} + \left( \frac{m_2}{N_2} \right)^{2p} + \left( \frac{m_3}{N_3} \right)^{2p} \right) \right), \quad (46) $$

The damping parameter $\beta$ is defined as $\beta = -\log \epsilon_M$, where $\epsilon_M$ is a low value parameter that was optimized by Djaka et al. (2015). For the applications in Section 4, $\epsilon_M = 0.2$ and $p = 1$. In Eq. (45), the time step $\Delta t$ must satisfy the Courant–Friedrichs–Lewy (CFL) limit:

$$ \Delta t_{CFL} = c \frac{\delta}{v_{max}} \quad (47) $$

where c is a user-specified fraction (here c = 0.25), $\delta$ is the voxel size and $v_{max}$ is the maximal GND velocity. Finally, the time step is given by $\Delta t = \min(\Delta t_{CFL}, \Delta t_e)$ where $\Delta t_{CFL}$ is defined in Eq. (35) and the time step $\Delta t_e$ is the classic time step used in the EVPFFT formulation.

Finally, to implement the MFDM-EVPFFT formulation, we need to compute the backstress tensor $\mathbf{B}$ introduced in Eq. (18). This is

performed here by computing its rate $\dot{\mathbf{B}}$, using Eqs. (18) and (9):

$$ \dot{\mathbf{B}} = -A \text{curl } \hat{\alpha} = A \text{curl } \mathbf{curl} \, \dot{\mathbf{U}}^p \quad (48) $$

Then, it is further calculated with the identity $\text{curl curl} = \text{grad div} - \text{div grad}$, then:

$$ \dot{\mathbf{B}} = A (\text{grad div } \dot{\mathbf{U}}^p - \text{div grad } \dot{\mathbf{U}}^p) \quad (49) $$

The last equation, written in Cartesian component form, is solved in the Fourier space:

$$ \hat{\mathbf{B}}_{ij} = -A \left( \hat{\mathbf{U}}_{ik}^p \xi_k \xi_j - \hat{\mathbf{U}}_{ij}^p \xi_k^2 \right) \quad (50) $$

The need for better accuracy of local fields, avoiding Gibbs phenomenon or aliasing, motivates the development of finite difference schemes to calculate the corresponding multipliers in the Fourier space between continuous and discrete Fourier transforms for partial derivatives of second order present in Eq. (50). Therefore, a successful numerical strategy already developed by Berbenni et al. (2014), for the spectral resolution of elasto-static equations of field dislocation and g-disclination mechanics, and by Lebensohn and Needleman (2016) for the spectral resolution of the Gurtin (2002)'s theory through the SG-EVPFFT formulation, consists in approximating the second derivatives in Cartesian space using a centered finite difference (FD) scheme (27-pixel approximation in 3D), and taking discrete Fourier transforms to these FD expressions. Then, following Berbenni et al. (2014), the expressions of $-\xi_k \xi_j$ and $-\xi_k^2$ are approximated by:

$$ \begin{aligned} & -\xi_k \xi_j \leftrightarrow \frac{1}{2\delta_k \delta_j} \left[ \cos \left( 2\pi \left( \frac{m_k}{N_k} + \frac{m_j}{N_j} \right) \right) - \cos \left( 2\pi \left( \frac{m_k}{N_k} - \frac{m_j}{N_j} \right) \right) \right] \\ &= \frac{1}{\delta_k \delta_j} \sin \left( \frac{2\pi m_k}{N_k} \right) \sin \left( \frac{2\pi m_j}{N_j} \right) \\ & - \xi_k^2 \leftrightarrow \frac{2}{\delta_k^2} \left( \cos \left( \frac{2\pi m_k}{N_k} \right) - 1 \right) \end{aligned} \quad (51) $$

where, $m_j$ were defined in Eq. (38) or Eq. (39).

Similar to the numerical procedure described in Lebensohn and Needleman (2016), the backstress tensor $\mathbf{B}$ known at time $t$ is used in the above iterative procedure based on AL scheme at $t + \Delta t$. Then, it is updated according to $\mathbf{B} = \mathbf{B}^t + \mathbf{B}\Delta t$. This means that a small time step is used to ensure numerical convergence (see Section 4.1).

# 4. Applications to two-phase laminate composites and polycrystalline aggregates

## 4.1. Simulations set-up: microstructures, materials parameters and loading conditions

A series of simulations on different unit cells were conducted with our MFDM-EVPFFT with backstress implementation. Two different microstructures were investigated (Fig. 1): a laminate composite microstructure following the works of Djaka et al. (2020) and Berbenni and Lebensohn (2021), and, a polycrystalline microstructure following the work of Berbenni et al. (2020).

In Fig. 1(a), a two-phase periodic laminate composite is considered to have a purely elastic phase (called second phase) and an elasto-viscoplastic phase with face-centered cubic (FCC) crystalline structure (called plastic channel). The period $H$ is set up in the y-direction (corresponding to $e_2$, $H = H_2$) with $H = d + 2h$, where d is the plastic channel size and 2h is the total size of elastic second phase along the y-direction.

Two representative volume elements (RVEs) will be defined for this microstructure. The first one is constituted of a volume fraction of second phase of 0.0625, with a spatial resolution of $N_2 = 128$ voxels along the $e_2$-direction for monotonic shear loading. Therefore, the unit cell is submitted to mixed strain-stress boundary conditions
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Simulation periodic unit cells used for MFDM-EVPFFT: (a) Two-phase laminate composite with period $H = H_2$ constituted of a plastic channel (red color) of size $d$ and elastic plates (second phase in blue color), subjected to pure shear loading, (b) Polycrystalline aggregate with $64 \times 64 \times 64$ voxels (27 grains) and period $H = H_1 = H_2 = H_3$ subjected to tensile loading with average grain size $d = H/\sqrt{27} = H/3$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

and applied shear strain rate $\dot{E}_{13} = \dot{E}_{31} = 0.001 \text{ s}^{-1}$, with $\Delta t_{\epsilon} = 0.01 \text{ s}$. The second one is constituted of a volume fraction of second phase of 0.125, with a spatial resolution of $N_2 = 32$ voxels along the $e_2$-direction for reversible shear loading. The applied shear strain rates are $\dot{E}_{13} = \dot{E}_{31} = \pm 0.001 \text{ s}^{-1}$, with $\Delta t_{\epsilon} = 0.01 \text{ s}$. Due to infinite planar boundaries along the $e_1$- and $e_3$-directions, only 2 voxels are sufficient along these directions due to unit cell periodicity.

The single crystallographic orientation of the plastic channel is given by the three following Bunge–Euler angles: $\phi_1 = 300^{\circ}$, $\phi = 54.7358^{\circ}$, $\phi_2 = 45^{\circ}$. Hence, a predominant single slip mode corresponding to the slip system $(111)\langle 101 \rangle$ is activated in the plastic channel as in Djaka et al. (2020) and Berbenni and Lebensohn (2021). As a result, a Schmid factor of 0.5 is obtained for this single crystalline orientation under this specific shear loading. This active slip system is reported on Fig. 1(a) by two arrows in the plastic channel, representing the unit normal $n^s$, and, the unit vector in the slip direction $b^s$. They are parallel to $e_3$ and $e_1$, respectively. In this configuration, the main GND density created during this forward–reverse shear loading is due to the screw component of the Nye tensor $\alpha_{11}$, with GND density defined by $\alpha_{11}/b$ (in $\text{m}^{-2}$).

In Fig. 1(b), a polycrystalline aggregate made of 27 randomly oriented grains is considered, having FCC crystalline structure within the periodic unit cell. The average grain size $d$ of the polycrystal is related to the period $H$ of the unit cell as $d = H/\sqrt{27} = H/3$ (Lebensohn and Needleman, 2016). The unit cell contains $64 \times 64 \times 64$ voxels. A tension–compression loading is applied in the $e_3$-direction at a macroscopic strain rate $\dot{E}_{33} = 10^{-3} \text{ s}^{-1}$ in tension (resp. $\dot{E}_{33} = -10^{-3} \text{ s}^{-1}$ in compression) with mixed macroscopic strain–stress boundary conditions.

For both microstructures, elasticity is assumed to be isotropic and homogeneous with Young's modulus $E = 69 \text{ GPa}$ and Poisson ratio $\nu = 0.33$ for Aluminum (Al) (elastic shear modulus $\mu = 25.9 \text{ GPa}$). For plastic channel and grains in Fig. 1(a)(b), the material parameters related to slip rule, GND velocity: $\dot{\gamma}^0$, $n$, $b$ (Burgers vector magnitude), $\eta$, and hardening model: $\tau_0$, $\tau_s$, $\theta_0$, $k_0$, are reported in Table 1 and are consistent with Al. Among these parameters, $k_0$ describes the magnitude of dissipative hardening. It is here chosen following the value fitted by Acharya and Beaudoin (2000) for the case of FCC metals, i.e. $k_0 = 20$, which was also taken in Djaka et al. (2020), Berbenni et al.

<table><caption>**Table 1**
List of material parameters used for simulations.</caption><thead><tr><th>γ̇<sup>0</sup> (s<sup>-1</sup>)</th><th>n</th><th>η</th><th>b (m)</th><th>τ<sub>0</sub> (MPa)</th><th>τ<sub>s</sub> (MPa)</th><th>θ<sub>0</sub> (MPa)</th><th>k<sub>0</sub></th></tr></thead><tbody><tr><td>1</td><td>20</td><td>0.33</td><td>2.86 × 10<sup>-10</sup></td><td>3</td><td>12</td><td>150</td><td>20 or 0</td></tr></tbody></table>

(2020) and Berbenni and Lebensohn (2021). An increase of $k_0$ means a decrease of the geometric mean-free path due to forest GND.

Here, the jump condition derived by Acharya (2007) to the case of an impenetrable interface (i.e. impenetrable phase boundaries to dislocations between plastic channel and elastic second phase) has been implemented in the MFDM-EVPFFT formulation only for the two-phase laminate microstructure (Fig. 1(a)) following the previous works of the present authors (Djaka et al., 2020; Berbenni and Lebensohn, 2021). Therefore, the interfacial condition (see Eq. (10)) (Acharya, 2007; Puri et al., 2011) yields in this case: $\dot{U}^p \times \mathbf{n} = 0$. Using the Cartesian coordinates in Fig. 1(a), the unit vectors along the y-, z-, x-directions are defined as $\mathbf{e}_2 = \mathbf{n}$, $\mathbf{e}_3 = \mathbf{l}$ and $\mathbf{e}_1 = \mathbf{n} \times \mathbf{l}$, respectively. Thus, the constraint on $\dot{U}^p$ is only applied to the voxels of the plastic channel adjacent to the interface, which leads to $\dot{U}_{12}^p = 0$ and $\dot{U}_{32}^p = 0$. This corresponds to a single slip situation with slip plane parallel to the interface and non zero slip rate ($\dot{\gamma}^s \neq 0$). Such condition was also discussed by Gurtin and Needleman (2005), which they called micro-hard condition. For polycrystals, a more difficult investigation is needed because more possibilities are present (impenetrable or penetrable GBs or intermediate regimes). Furthermore, the unit normal is dependent on grain boundary normal distributions in the polycrystalline aggregate. In the SG-EVPFFT context, Zecevic et al. (2023) and Zecevic et al. (2024) have recently implemented boundary voxels to solve interfacial problems in nanolaminates and polycrystals. Here, no specific interface condition will be implemented at grain boundaries of polycrystals and the use of the boundary voxels is left for future study.

Finally, the internal material length scale $l$ present in the defect energy through the parameter $A = \tau_0 l^2$ will be normalized by the unit cell period $H$ following different previous contributions in strain-gradient plasticity (Bittencourt et al., 2003; Christodoulou et al., 2022). Therefore, we will study the internal length effect through the normalized length scale $\bar{l} = \frac{l}{H}$ in the following results where $l = \sqrt{\frac{A}{\tau_0}}$.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Monotonic overall shear responses $<\sigma_{13}>$ vs. $<\epsilon_{13}>$ predicted by the non-local MFDM-EVPFFT model: Effect of the normalized internal material length scale $\bar{l} = l/H$ on the overall shear response for a two-phase laminate composite with a fixed volume fraction of second phase of 0.0625.</figcaption></figure>

We will set $k_0 = 0$ when considering energetic hardening only in the model. We will set $\bar{l} = 0$ when considering dissipative hardening only in the model. We will set $k_0 = 20$ and $\bar{l} \neq 0$ when considering both energetic hardening and dissipative hardening in the model, which is called combined hardening.

## 4.2. Results for two-phase composites: monotonic shear loading responses

First of all, the effect of the normalized length scale $\bar{l} = l/H$ on the overall shear responses of the two-phase composite is investigated until 1% shear strain in Fig. 2, assuming energetic hardening only. Here, the unit cell period $H$ along $e_2$ is fixed to $1 \mu m$, and the volume fraction of second phase is 0.0625. It is seen that the overall flow stress is increasing when $\bar{l}$ varies from 0.05 to 0.5, which is consistent with a smaller is stronger effect (at a given $l$). This result is similar to those of Lebensohn and Needleman (2016) for bi-crystals.

In Fig. 3, the results of four different models are studied for the same composite: (i) the local model corresponding to classic CP-EVPFFT (Lebensohn et al., 2012), (ii) the present non-local model corresponding to MFDM-EVPFFT with only energetic hardening ($\bar{l} = \frac{l}{H} \approx 0.5$, $k_0 = 0$), (iii) the non-local model corresponding to MFDM-EVPFFT with only dissipative hardening (Djaka et al., 2020; Berbenni et al., 2020; Berbenni and Lebensohn, 2021) ($\bar{l} = 0$, $k_0 = 20$), (iv) the non-local model corresponding to MFDM-EVPFFT with combined energetic hardening and dissipative hardening ($\bar{l} = \frac{l}{H} \approx 0.5$, $k_0 = 20$). For classic CP-EVPFFT, the strong overall hardening is due to the classic composite effect (i.e. kinematical hardening due to stress accommodation between both phases), see e.g. Berbenni and Lebensohn (2021) for a detailed study of this hardening type. Non-local models enhance hardening and it is interesting to note that for this material length parameter ($\bar{l} \approx 0.5$), the purely energetic and purely dissipative models give approximately the same hardening rate. The combined hardening model gives the most strain hardening.

Following Fig. 3 with the same color code, Fig. 4 show the local and the non-local predictions after 1% shear strain for shear stress component $\sigma_{13}$ (Fig. 4(a)), screw dislocation density $\alpha_{11}/b$ (Fig. 4(b)), backstress field defined by $\tau_{kin}^{(1)} = m^{(1)} : B$ (Fig. 4(c)), slip $\gamma^{(1)}$ (Fig. 4(d)), where (1) denotes the active slip system in the channel. In contrast with the local model (classic CP-EVPFFT) where stress and slip profiles are piece-wise uniform (see Fig. 4(a)(d)), the non-local models predict the formation of shear stress $\sigma_{13}$ and slip gradients $\gamma^{(1)}$ in the plastic channel with a boundary layer region in which these fields smoothly vary between channel center and phase boundaries. It is interesting to note that the new MFDM-EVPFFT model with energetic hardening predicts

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Monotonic overall shear responses $\langle\sigma_{13}\rangle$ vs. $\langle\epsilon_{13}\rangle$ for four different models: (i) local model (CP-EVPFFT), (ii) non-local MFDM-EVPFFT models with energetic hardening ($\bar{l} = \frac{l}{H} \approx 0.5$, $k_0 = 0$), (iii) non-local MFDM-EVPFFT model with dissipative hardening ($\bar{l} = 0$, $k_0 = 20$), (iv) non-local MFDM-EVPFFT with combined energetic and dissipative hardening ($\bar{l} = \frac{l}{H} \approx 0.5$, $k_0 = 20$).</figcaption></figure>

a smoother variation of both fields as compared to the MFDM-EVPFFT model, which gives less slip gradient. In Fig. 4(c), the backstress field, only present in the MFDM-EVPFFT model with energetic or combined hardening, exhibits peaks close to phase boundaries with opposite signs. These peaks satisfy self-equilibrium: a positive sign is observed in the plastic channel with a smooth decay and a negative sign is observed in the elastic second phase with a sharper decay with more intense values. Here, it is interesting to make a comparison with the results of Lebensohn and Needleman (2016) (see figure 2 therein), where such backstress peaks at interfaces between soft and hard grains were reported. However, let us note that their convention for backstress field was derived from the gradients of micro-stress following the Gurtin's theory (Gurtin, 2002; Bittencourt et al., 2003), not from $\tau_{kin}^{(1)}$ as here. Fig. 4(b) first shows that the CP-EVPFFT model gives very low values of screw GND density with a sharp profile only concentrated at interfaces, i.e. interfacial dislocations. In contrast, both MFDM-EVPFFT formulations without or with backstress are able to predict double ended pile-ups of screw GND density $\alpha_{11}/b$ (note the different polarities at both sides). In the case of the MFDM-EVPFFT model with energetic or combined hardening, larger values of GND densities with smoother profiles are obtained as compared to the MFDM-EVPFFT model without backstress.

## 4.3. Results for two-phase composites: reverse shear loading responses

Let us consider now the reversible shear response loading of a two-phase composite, where the volume fraction of elastic second phase is now set to 0.125 with a unit cell period of $H = 0.25 \mu m$ (which means a smaller channel size $d$ as before). The role of energetic hardening on the Bauschinger effect is here studied considering the material length scale parameter $\bar{l} = \frac{l}{H} \approx 0.9$ (which is equivalent to $\frac{l}{d} = 1$). The predictions of the four different models are also studied for this configuration. Fig. 5 shows the overall shear strain vs. shear strain response after a forward loading up to 0.2% shear strain followed by a reverse loading up to -0.2%. Here, the MFDM-EVPFFT model with dissipative hardening model leads to higher strain-hardening rate during the forward loading stage as compared to the MFDM-EVPFFT model with energetic hardening. Among all models, the MFDM-EVPFFT model with combined hardening displays the largest flow stress.

The Bauschinger effect can be characterized by the quantity $\Sigma_b = (\Sigma_f - \Sigma_r)/2$ (Atkinson et al., 1974; Corbin et al., 1996), where $\Sigma_f$ is here the forward shear stress $\langle\sigma_{13}\rangle$ at $\langle\epsilon_{13}\rangle = 0.2\%$ (in the plastic regime) and $\Sigma_r$ is the reverse shear stress $\langle\sigma_{13}\rangle$ at the onset of plasticity during
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Predictions of the four models after 1% shear strain for: shear stress component $\sigma_{13}$ (a), screw dislocation density $\alpha_{11}/b$ (b), backstress field defined by $\tau_{kin}^{(1)} = m^{(1)} : B$ (c), slip $\gamma^{(1)}$ (d), where (1) denotes the active slip system in the channel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

reverse loading. Fig. 5 first shows that a strong Bauschinger effect is already predicted with the CP-EVPFFT model due to the composite effect, and a sharp elastic plastic transition is observed during reverse loading. This result is consistent with the model of Mughrabi (1988). In contrast, the MFDM-EVPFFT model with dissipative hardening predicts a slower elastic-plastic transition as already observed in Berbenni and Lebensohn (2021). Fig. 5 shows that the MFDM-EVPFFT model with only energetic hardening predicts an increase of the Bauschinger effect as compared to the CP-EVPFFT model. Considering a combined hardening with the MFDM-EVPFFT model, it predicts a further increase of the Bauschinger effect. This increase of the linear kinematic hardening observed with the MFDM-EVPFFT model including a backstress due to $\tau_{kin}^{(1)}$, can be interpreted in comparison with the dislocation model of Tanaka and Mura (1981), based on the static equilibrium of a double-ended edge dislocation pile-up (Forest, 2008).

The evolutions of backstress field and screw GND density, are studied at the four different states denoted A, B, C, D in Fig. 5. Fig. 6 shows the spatial variation along the e₂-direction (in voxels) of the backstress field, while Fig. 7 shows the spatial variation of the screw GND density component $\alpha_{11}/b$. Regarding the MFDM-EVPFFT model with both energetic or combined hardening, Fig. 6 shows a change of sign of backstress peaks with similar maximal intensities in the plastic channel ($\approx 6.2$ MPa) between states A and D. In between, a decrease of backstress peak intensities is observed at state B, i.e during the first stage of reverse loading, while it is canceled at state C for $\langle\epsilon_{13}\rangle = 0$. This change of sign can be attributed to the evolution of the GND density. For both MFDM-EVPFFT without or with backstress formulation, Fig. 7 shows the formation of a continuous screw double-ended GND pile up at the end of the first forward shear loading (state A), with maximal magnitudes for $\alpha_{11}/b$ near phase boundaries reaching $8.2 \times 10^{13} \text{ m}^{-2}$. At states B and C, an unpiling-up mechanism operates with the progressive annihilation of the screw GND density first formed at state A, and the built-up of a GND pile up with opposite polarities reaching maximal GND densities at state D.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Prediction of the reversible behavior in shear deformation with four different models for a two-phase laminate composite. The second phase volume fraction is set to 0.125, with a unit cell period $H = 0.25$ µm: (i) local model (CP-EVPFFT), (ii) non-local MFDM-EVPFFT models with energetic hardening ($\bar{l} \approx 0.9, k_0 = 0$), (iii) non-local MFDM-EVPFFT model with dissipative hardening ($\bar{l} = 0, k_0 = 20$), (iv) non-local MFDM-EVPFFT with combined hardening ($\bar{l} \approx 0.9, k_0 = 20$).</figcaption></figure>

### 4.4. Results for polycrystals: reversible tension–compression responses

Let us finally study the reversible tension–compression loading response along the e₃-direction of a polycrystal constituted of 27 randomly oriented grains. For this example, $H = 3$ µm, which means an average grain size $d \approx 1$ µm. Again, the role of energetic hardening on the Bauschinger effect is compared to the results obtained with the MFDM-EVPFFT model developed in Berbenni et al. (2020) for FCC polycrystal plasticity. The normalized material length scale parameter $\bar{l} = l/H$ is set to $\bar{l} \approx 0.6$.

Fig. 8 shows some comparisons in terms of stress–strain responses between CP-EVPFFT, MFDM-EVPFFT with dissipative hardening, energetic hardening or combined hardening, for the FCC polycrystal. As for
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Spatial variations along the $e_2$-direction (in voxels) of the backstress field at four different states: A, B, C, D (see Fig. 5) for three models: MFDM-EVPFFT model with dissipative hardening, MFDM-EVPFFT model with energetic hardening, MFDM-EVPFFT model with combined hardening.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Spatial variations along the $e_2$-direction (in voxels) of the GND density $\alpha_{11}/b$ at four different states: A, B, C, D (see Fig. 5) for three models: MFDM-EVPFFT model with dissipative hardening, MFDM-EVPFFT model with energetic hardening, MFDM-EVPFFT model with combined hardening.</figcaption></figure>

the composite material, the Bauschinger effect is characterized by the quantity $\Sigma_b = (\Sigma_f - \Sigma_r)/2$, where $\Sigma_f$ is now the forward tensile stress $\langle\sigma_{33}\rangle$ at $\langle\epsilon_{33}\rangle = 0.2\%$ (in the plastic regime) and $\Sigma_r$ is the compressive stress $\langle\sigma_{33}\rangle$ at the onset of plasticity during compressive loading. Then, the mechanical test was stopped at $\langle\epsilon_{33}\rangle = -0.08\%$ during compressive loading. The predicted tension–compression responses are reported on Fig. 8 for the four investigated models. The lowest Bauschinger effect is obtained with conventional CP-EVPFFT and corresponds to classic inter-granular hardening in polycrystal plasticity. Fig. 8 shows that the MFDM-EVPFFT model with only energetic hardening predicts an increase of the Bauschinger effect as compared to the CP-EVPFFT model. Considering a combined hardening with the MFDM-EVPFFT model, it predicts a further increase of the Bauschinger effect. Therefore, conclusions similar to the case of the two-phase composite can be drawn. However, all the stress–strain responses exhibit non-linear hardening due to polycrystalline plasticity, which is different from the results obtained with the composite. Under reversible tension–compression loading, the Bauschinger stress is thus increased by adding energetic hardening, which is in agreement with the results of Christodoulou

et al. (2022) who used HO-SGP plasticity theory following Gurtin’s theory (Gurtin, 2002).

Two states denoted A and B are also indicated in figure 8, to compare different local mechanical fields for MFDM-EVPFFT model with dissipative hardening and energetic hardening. Figs. 9 and 10 report the evolutions of three local fields for both models. The fields are the longitudinal stress component along the $e_3$ direction ($\sigma_{33}$), the plastic strain component along the $e_3$ direction ($\epsilon_{33}^p$) and the scalar GND density $\rho_{GND}$ taken as the $L^2$ norm of the Nye tensor $\alpha$ normalized by the Burgers vector magnitude: $\rho_{GND} = \frac{\sqrt{\alpha : \alpha}}{b}$. In both figures, it is observed at state A that the regions near grain boundaries are decorated by large GND density pile-ups and plastic strain gradients from grain center to grain boundaries with both dissipative and energetic hardening. This is in contrast to local CP-EVPFFT for which the results were reported in Berbenni et al. (2020), where no long-range GND pile-ups are observed. These GND density pile ups are reduced in magnitude at state B, due to a transient regime between both states A and B, with annihilation of first GNDs formed during forward loading, followed
<figure><img src="image_10.png" /><figcaption>Fig. 8. Prediction of tension–compression tests obtained with four different models for a polycrystal made of 27 grains: (i) local model (CP-EVPFFT), (ii) non-local MFDM-EVPFFT model with energetic hardening ($\bar{I} \approx 0.6, k_0 = 0$), (iii) non-local MFDM-EVPFFT model with dissipative hardening ($\bar{I} = 0, k_0 = 20$), (iv) non-local MFDM-EVPFFT with combined energetic and dissipative hardening ($\bar{I} \approx 0.6, k_0 = 20$). The states denoted A and B are also indicated to compare different local mechanical fields in Figs. 9 and 10.</figcaption></figure>

by a rebuilt of new GNDs during compression. It is also observed that the longitudinal stress hot spots are related to the GND density ones. In comparison with the MFDM-EVPFFT model with dissipative hardening (Fig. 9), a smoothing effect is obtained in Fig. 10 with the MFDM-EVPFFT model with energetic hardening, mainly on plastic strain and GND density distributions. This effect is clearly attributed to energetic hardening, which spreads more plastic strain and GND density as compared to dissipative hardening. It is noteworthy that a similar smoothing effect was also reported in the SG-EVPFFT simulation results for polycrystals in Lebensohn and Needleman (2016).

## 5. Concluding remarks

A new version of the MFDM-EVPFFT formulation was proposed to extend the works of Djaka et al. (2020), Berbenni et al. (2020) and Berbenni and Lebensohn (2021), by incorporating a *defect energy* in the thermodynamic framework of the MFDM. This internal length scale dependent defect energy has been assumed to have a quadratic form, which is at the origin of energetic hardening due to GND density. This

leads to a new backstress tensor B, which in contrast with other MFDM formulations using a core energy (Arora and Acharya, 2020a,b), is used in the expressions of the slip flow rule, with an effective shear stress on each slip system, and the expression of the Peach–Koehler glide force. As compared to the previous version of MFDM-EVPFFT (Djaka et al., 2020; Berbenni et al., 2020; Berbenni and Lebensohn, 2021), two new features are presented: (i) a new expression of the Jacobian for the augmented Lagrangian scheme considering plastic strain rate due to the backstress tensor, (ii) a spectral discrete method based on centered finite difference (FD) schemes to calculate the backstress tensor evolution.

As illustrations, two types of microstructures were studied: a two-phase laminate composite with plastic channel sandwiched between elastic plates, and a polycrystalline aggregate with randomly oriented grains. Furthermore, two types of loadings were simulated: monotonic and reversible loadings. For the two-phase laminate, increasing the internal length $l$ for fixed size of the unit cell $H$ was shown to increase strain hardening. Furthermore, it is shown that energetic hardening due to long-range backstress produces size dependent hardening. When this energetic hardening is combined with size dependent dissipative hardening (calibrated by the parameter $k_0 = 20$), the overall hardening rate is amplified. By analyzing the different intracrystalline mechanical fields like stress, GND density, backstress and slip profiles in the plastic channel, it is shown that a long-range intracrystalline backstress superimposes to the classic backstress due to the composite effect, which leads to larger slip and stress gradients. In the case of reversible shear loading, an enhanced Bauschinger effect is produced by the non-local MFDM-EVPFFT formulation with both energetic and dissipative hardening, as compared to the local CP-EVPFFT and previous MFDM-EVPFFT formulations. A piling up/unpiling up mechanism of screw GND density is observed as in Berbenni and Lebensohn (2021), and the evolution of associated backstress peaks with different polarities were observed. It is shown that energetic hardening and combined energetic/dissipative hardenings increase the spreading width of double-ended screw GND density pile ups.

Numerical simulations with a polycrystalline RVE obtained with periodic Voronoi tessellation with a grain size of $d = 1$ µm were performed with 64 × 64 × 64 voxels (27 grains). Under reversible tension–compression loading, a larger Bauschinger stress is obtained with energetic hardening, which is in agreement with the results of Christodoulou et al. (2022), even though the present MFDM-EVPFFT formulation is not based on Higher Order strain gradient (HO-SGP)

<figure><img src="image_11.png" /><figcaption>Fig. 9. From left to right: evolutions of $\sigma_{33}$ (in MPa), $\epsilon_{33}^p$, and $\rho_{GND}$ (in m$^{-2}$), during a tension–compression test using the MFDM-EVPFFT model with dissipative hardening. Both states A (top line) and B (bottom line) are indicated in Fig. 8.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** From left to right: evolutions of $\sigma_{33}$ (in MPa), $\epsilon^{p}_{33}$, and $\rho_{GND}$ (in m$^{-2}$), during a tension–compression test using the MFDM-EVPFFT model with energetic hardening. Both states A (top line) and B (bottom line) are indicated in Fig. 8.</figcaption></figure>

plasticity theory as in Gurtin's theory (Gurtin, 2002). As opposed to local CP-EVPFFT, large GND density pile ups are observed from grain center to grain boundaries with the present MFDM-EVPFFT formulation. From these full-field polycrystalline simulations, it is shown that, as compared to the previous MFDM-EVPFFT formulation with dissipative hardening, the incorporation of energetic hardening in MFDM-EVPFFT formulation spreads more GND densities over the grains. Furthermore, the smoothing of the fields with energetic hardening, is also consistent with the SG-EVPFFT formulation of Lebensohn and Needleman (2016).

As a perspective to this work, the present MFDM-EVPFFT model with combined energetic and dissipative hardening can be applied to polycrystals with different slip assumptions at grain boundaries. Since new robust FFT-based solvers for crystal elasto-viscoplasticity at finite strains for non-local plasticity are now available (Zecevic et al., 2023, 2024), the MFDM-EVPFFT formulation can also be extended to a large deformation framework using the constitutive finite deformation MFDM equations (Acharya, 2004; Arora and Acharya, 2020a; Arora et al., 2022). Finally, only a quadratic form of defect energy was used. However, non-quadratic gradient-related defect energy forms (e.g. power laws with $n \neq 2$ or hyperbolic trigonometric functions) give better results in the context of Strain Gradient Plasticity (SGP), see e.g. Bayerschen and Böhlke (2016), Jebahi et al. (2020), Zecevic et al. (2023). Therefore, other forms of defect energy can be studied to explore backstress profiles and GND density patterns during reversible plasticity in both small and finite deformations settings with the MFDM-EVPFFT formulation.

### CRediT authorship contribution statement

**Stéphane Berbenni:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**Vincent Taupin:** Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Conceptualization.
**Ricardo A. Lebensohn:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

SB, VT thank the French State (ANR) through the program "Investment in the future" (LabEx "DAMAS" referenced as ANR-11-LABX-0008-01) for financial support. RAL work was funded by Los Alamos National Laboratory under Laboratory-Directed Research and Development (LDRD) program.

### Data availability

Data will be made available on request.

### References

Acharya, A., 2001. A model of crystal plasticity based on the theory of continuously distributed dislocations. J. Mech. Phys. Solids 49, 761–784.

Acharya, A., 2003. Driving forces and boundary conditions in continuum dislocation mechanics. Proc. R. Soc. A 459, 1343–1363.

Acharya, A., 2004. Constitutive analysis of finite deformation field dislocation mechanics. J. Mech. Phys. Solids 52, 301–316.

Acharya, A., 2007. Jump condition for GND evolution as a constraint on slip transmission at grain boundaries. Phil. Mag. 87, 1349–1359.

Acharya, A., 2010. New inroads in an old subject: Plasticity, from around the atomic to the macroscopic scale. J. Mech. Phys. Solids 58 (5), 766–778.

Acharya, A., 2011. Microcanonical entropy and mesoscale dislocation mechanics and plasticity. J. Elasticity 104, 23–44.

Acharya, A., Bassani, J.L., 2000. Lattice incompatibility and a gradient theory of crystal plasticity. J. Mech. Phys. Solids 48 (8), 1565–1595.

Acharya, A., Beaudoin, A.J., 2000. Grain size effect in viscoplastic polycrystals at moderate strains. J. Mech. Phys. Solids 48, 2213–2230.

Acharya, A., Roy, A., 2006. Size effects and idealized dislocation microstructure at small scales : Predictions of a phenomenological model of mesoscopic field dislocation mechanics : Part I. J. Mech. Phys. Solids 54, 1687–1710.

Acharya, A., Roy, A., Sawant, A., 2006. Continuum theory and methods for coarse-grained plasticity. Scr. Mater. 54, 705–710.

Aifantis, E.C., 1984. On the microstructural origin of certain inelastic models. Trans. ASME, J. Eng. Mater. Technol. 106, 326–330.

Aifantis, E.C., 1987. The physics of plastic deformation. Int. J. Plast. 3, 211–247.

Allain-Bonasso, N., Wagner, F., Berbenni, S., Field, D.P., 2012. A study of the heterogeneity of plastic deformation in IF steel by EBSD. Mater. Sci. Eng. A 548, 56–63.

Amouzou-Adoun, Y.A., Jebahi, M., Fivel, M., Forest, S., Lecomte, J.-S., Schuman, C., Abed-Meraim, F., 2023. On elastic gaps in strain gradient plasticity: 3D discrete dislocation dynamics investigation. Acta Mater. 252, 118920.

Arora, R., Acharya, A., 2020a. Dislocation pattern formation in finite deformation crystal plasticity. Int. J. Solids Struct. 184, 114–135.
Arora, R., Acharya, A., 2020b. A unification of finite deformation J2 Von-Mises plasticity and quantitative dislocation mechanics. J. Mech. Phys. Solids 143, 104050.

Arora, A., Arora, R., Acharya, A., 2022. Mechanics of micropillar confined thin film plasticity. Acta Mater. 238, 118192.

Arora, A., Arora, R., Acharya, A., 2023. Interface-dominated plasticity and kink bands in metallic nanolaminates. Crystals 13 (5).

Asaro, R., 1983. Crystal plasticity. ASME J. Appl. Mech. 50, 921-934.

Ashby, M.F., 1970. Deformation of plastically non-homogeneous materials. Phil. Mag. 21, 399-424.

Atkinson, J.D., Brown, L.M., Stobbs, W.M., 1974. The work-hardening of copper-silica. IV The Bauschinger effect and plastic relaxation. Phil. Mag. 30, 1247-1280.

Bayerschen, E., Böhlke, 2016. Power-law defect energy in a single-crystal gradient plasticity framework - A computational study. Comput. Mech. 31, 13-27.

Berbenni, S., Lebensohn, R.A., 2021. A numerical study of reversible plasticity using continuum dislocation mechanics. C. R. Physique 22 (S3), 295-312.

Berbenni, S., Taupin, V., Djaka, K.S., Fressengeas, C., 2014. A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics. Int. J. Solids Struct. 51, 4157-4175.

Berbenni, S., Taupin, V., Lebensohn, R.A., 2020. A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals. J. Mech. Phys. Solids 135, 103808.

Bittencourt, E., Needleman, A., Gurtin, M.E., VanderGiessen, E., 2003. A comparison of non local continuum and discrete dislocation plasticity predictions. J. Mech. Phys. Solids 51, 281-310.

Busso, E., Meissonnier, F., O'Dowd, N., 2000. Gradient-dependent deformation of two-phase single crystals. J. Mech. Phys. Solids 48 (11), 2333-2361.

Cheong, K., Busso, E., Arsenlis, A., 2005. A study of microstructural length scale effects on the behaviour of FCC polycrystals using strain gradient concepts. Int. J. Plast. 21 (9), 1797-1814.

Christodoulou, P.G., Lebensohn, R.A., Beyerlein, J.L., 2022. Study of the interplay between lower-order and higher-order energetic strain-gradient effects in polycrystal plasticity. J. Mech. Phys. Solids 164, 104906.

Coleman, B.D., Gurtin, M.E., 1967. Thermodynamics with internal state variables. J. Chem. Phys. 47, 597-613.

Coleman, B.D., Noll, W., 1963. The thermodynamics of elastic materials with heat conduction and viscosity. Arch. Ration. Mech. Anal. 13, 167-178.

Corbin, S.F., Wilkinson, D.S., Embury, J.D., 1996. The Bauschinger effect in a particulate reinforced Al alloy. Mater. Sci. Eng. A 207, 1-11.

Cordero, N.M., Forest, S., Busso, E.P., Berbenni, S., Cherkaoui, M., 2012. Grain size effects on plastic strain and dislocation density tensor fields in metal polycrystals. Comput. Mater. Sci. 52, 7-13.

Cordero, N.M., Gaubert, A., Forest, S., Busso, E.P., Galerneau, F., Kruch, S., 2010. Size effects in generalised continuum crystal plasticity for two-phase laminates. J. Mech. Phys. Solids 58, 1963-1994.

Cordero, Z.C., Knight, B.E., Schuh, C.A., 2016. Six decades of the Hall-Petch effect: a survey of grain size strengthening studies on pure metals. Int. Mater. Rev. 61 (8), 495-512.

Das, A., Acharya, A., Suquet, P., 2016. Microstructure in plasticity without nonconvexity. Comput. Mech. 57, 387-403.

Djaka, K.S., Berbenni, S., Taupin, V., Lebensohn, R.A., 2020. A FFT-based numerical implementation of mesoscale field dislocation mechanics: Application to two-phase laminates. Int. J. Solids Struct. 184, 136-152.

Djaka, K.S., Taupin, V., Berbenni, S., Fressengeas, C., 2015. A numerical spectral approach to solve the dislocation density transport equation. Modelling Simul. Mater. Sci. Eng. 23, 065008, (27pp).

Evers, L.P., Brekelmans, W.A.M., Geers, M.G.D., 2004. Non-local crystal plasticity model with intrinsic SSD and GND effects. J. Mech. Phys. Solids 52, 2379-2401.

Evers, L.P., Parks, D.M., Brekelmans, W.A.M., Geers, M.G.D., 2002. Crystal plasticity model with enhanced hardening by geometrically necessary dislocation accumulation. J. Mech. Phys. Solids 50, 2403-2424.

Fleck, N., Muller, G., Ashby, M., Hutchinson, J., 1994. Strain gradient plasticity: Theory and experiment. Acta Metall. Mater. 42 (2), 475-487.

Fleck, N.A., Willis, J.R., 2015. Strain gradient plasticity: energetic or dissipative? Acta Mech. Sin. 31, 465-472.

Forest, S., 2008. Some links between Cosserat, strain gradient crystal plasticity and the statistical theory of dislocations. Phil. Mag. 88, 30-32.

Genée, J., Gey, N., Bonnet, F., Lebensohn, R.A., Berbenni, S., 2021. Experimental and numerical investigation of key microstructural features influencing the localization of plastic deformation in Fe-TiB2 metal matrix composite. J. Mater. Sci. 56 (19), 11278-11297.

Gottlieb, D., Hesthaven, J.S., 2001. Spectral methods for hyperbolic problems. J. Comput. Appl. Math. 128, 83-131.

Gudmundson, P., 2004. A unified treatment of strain gradient plasticity. J. Mech. Phys. Solids 52, 1379-1406.

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

Gurtin, M., 2004. A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin. J. Mech. Phys. Solids 52, 2545-2568.

Gurtin, M.E., Anand, L., 2005. A theory of strain-gradient plasticity for isotropic, plastically irrotational materials, Part I: Small deformations. J. Mech. Phys. Solids 53, 2545-2568.

Gurtin, M.E., Anand, L., 2009. Thermodynamics applied to gradient theories involving accumulated plastic strain: The theories of Aifantis and Fleck and Hutchinson and their generalization. J. Mech. Phys. Solids 57, 405-421.

Gurtin, M.E., Needleman, A., 2005. Boundary conditions in small-deformation, single crystal plasticity that account for the Burgers vector. J. Mech. Phys. Solids 53, 1-31.

Hall, E.O., 1951. The deformation and ageing of mild steels. Proc. R. Soc. Lond. [Biol] 64, 747-753.

Han, C.-S., Gao, H., Huang, Y., Nix, W.D., 2005a. Mechanism-based strain gradient crystal plasticity—I. Theory. J. Mech. Phys. Solids 53 (5), 1188-1203.

Han, C.-S., Gao, H., Huang, Y., Nix, W.D., 2005b. Mechanism-based strain gradient crystal plasticity—II. Analysis. J. Mech. Phys. Solids 53 (5), 1204-1222.

Hochrainer, T., 2015. Multipole expansion of continuum dislocations dynamics in terms of alignment tensors. Phylosophical Mag. 95 (12), 1321-1367.

Hochrainer, T., Zaiser, M., Gumbsch, P., 2007. A three-dimensional continuum theory of dislocation systems: kinematics and mean-field formulation. Phil. Mag. 87 (8-9), 1261-1282.

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

Jebahi, M., Forest, S., 2023. An alternative way to describe thermodynamically-consistent higher-order dissipation within strain gradient plasticity. J. Mech. Phys. Solids 170, 105103.

Kiener, D., Motz, C., Grosinger, W., Weygand, D., Pippan, R., 2010. Cyclic response of copper single crystal micro-beams. Scr. Mater. 63 (5), 500-503.

Konijnenberg, P.J., Zaefferer, S., Raabe, D., 2015. Assessment of geometrically necessary dislocation levels derived by 3D EBSD. Acta Mater. 99, 402-414.

Kröner, E., 1958. Kontinuumstheorie der versetzungen und eigenspannungen. In: Collatz, L., Loesch, F. (Eds.), Ergebnisse der Angewandte Mathematik 5. Springer Verlag, Berlin.

Kröner, E., 1981. Continuum theory of defects. In: Balian, R., et al. (Eds.), Physics of Defects Les Houches Session 35. North Holland, New York, pp. 215-315.

Kröner, E., 1996. Dislocation theory as a physical field theory. Meccanica 31, 577-587.

Kubin, L.P., Canova, G., Condat, M., Devincre, B., Pontikis, V., Bréchet, Y., 1992. Dislocation microstructure and plastic flow : a 3-D simulation. Solid State Phenomena 23-24, 455-472.

Kucharski, S., Maj, M., Ryś, M., Petryk, H., 2024. Size effects in spherical indentation of single crystal copper. Int. J. Mech. Sci. 272, 109138.

Kuroda, M., Tvergaard, V., 2008. On the formulations of higher-order strain gradient crystal plasticity models. J. Mech. Phys. Solids 56 (4), 1591-1608.

Lebensohn, R., 2001. N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform. Acta Mater. 49, 2723-2737.

Lebensohn, R.A., Kanjarla, A.K., Eisenlohr, P., 2012. An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials. Int. J. Plast. 32-33, 59-69.

Lebensohn, R.A., Needleman, A., 2016. Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms. J. Mech. Phys. Solids 97, 333-351.

Lloyd, D.J., 1994. Particle reinforced aluminium an magnesium matrix composites. Int. Metall. Rev. 39, 1-23.

Mayeur, J., McDowell, D., 2014. A comparison of Gurtin type and micropolar theories of generalized single crystal plasticity. Int. J. Plast. 57, 29-51.

Mecking, H., Kocks, U.F., 1981. Kinetics of flow and strain-hardening. Acta Metall. 29, 1865-1875.

Michel, J.C., Moulinec, H., Suquet, P., 2001. A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int. J. Numer. Methods Eng. 52, 139-160.

Motz, C., Schöberl, T., Pippan, R., 2005. Mechanical properties of micro-sized copper bending beams machined by the focused ion beam technique. Acta Mater. 53 (15), 4269-4279.

Motz, C., Weygand, D., Senger, J., Gumbsch, P., 2009. Initial dislocation structures in 3-D discrete dislocation dynamics and their influence on microscale plasticity. Acta Mater. 57 (6), 1744-1754.

Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and non linear properties of composites. C. R. Acad. Sci. Paris II 318, 1417-1423.

Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Engrg. 157, 69-94.

Mughrabi, H., 1988. Dislocation clustering and long-range internal stresses in monotonically and cyclically deformed metal crystals. Rev. Phys. Appl. 23, 367-379.

Mura, T., 1963. Continuous distribution of moving dislocations. Phil. Mag. 89, 843-857.

Mura, T., 1987. Micromechanics of Defects in Solids, dordrecht ed. Kluwer Academic Plublishers.

Nan, C.W., Clarke, D.R., 1996. The influence of particle size and particle fracture on the elastic/plastic deformation of metal matrix composites. Acta Mater. 44, 3801-3811.
Nicola, L., Van der Giessen, E., Gurtin, M.E., 2005. Effect of defect energy on strain-gradient predictions of confined single-crystal plasticity. J. Mech. Phys. Solids 53 (6), 1280-1294.

Niordson, C.F., Kysar, J.W., 2014. Computational strain gradient crystal plasticity. J. Mech. Phys. Solids 62, 31-47.

Nix, W.D., Gao, H., 1998. Indentation size effects in crystalline materials: A law for strain gradient plasticity. J. Mech. Phys. Solids 46 (3), 411-425.

Nye, J.F., 1953. Some geometrical relations in dislocated crystals. Acta Mater. 1, 153-162.

Pantleon, W., 2008. Resolving the geometrically necessary dislocation content by conventional electron backscattering diffraction. Scr. Mater. 58 (11), 994-997.

Petch, N.J., 1953. The cleavage strength of polycrystals. J. Iron Steel Inst. 174, 25-28.

Puri, S., Acharya, A., Rollett, A.D., 2010. Controlling plastic flow across grain boundaries in a continuum model. Metall. Mater. Trans. A 42, 669-675.

Puri, S., Das, A., Acharya, A., 2011. Mechanical response of multicrystalline thin films in mesoscale field dislocation mechanics. J. Mech. Phys. Solids 59, 2400-2417.

Puri, S., Roy, A., 2012. Plastic deformation of multicrystalline thin films: Grain size distribution vs. grain orientation. Comput. Mater. Sci. 52, 20-24.

Puri, S., Roy, A., Acharya, A., Dimiduk, D., 2009. Modeling dislocation sources and size effects at initial yield in continuum plasticity. J. Mech. Mater. Struct. 4 (9), 1603-1618.

Richeton, T., Wang, G.F., Fressengeas, C., 2011. Continuity constraints at the interfaces and their consequences on the work hardening of metal-matrix composites. J. Mech. Phys. Solids 59, 2023-2043.

Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D.D., Raabe, D., 2010. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications. Acta Mater. 58, 1152-1211.

Roy, A., Acharya, A., 2005. Finite element approximation of field dislocation mechanics. J. Mech. Phys. Solids 53, 143-170.

Roy, A., Acharya, A., 2006. Size effects and idealized dislocation microstructure at small scales : Predictions of a phenomenological model of mesoscopic field dislocation mechanics : Part II. J. Mech. Phys. Solids 54, 1711-1743.

Roy, A., Puri, S., Acharya, A., 2007. Phenomenological mesoscopic field dislocation mechanics, lower-order gradient plasticity, and transport of mean excess dislocation density. Modelling Simul. Mater. Sci. Eng. 15, 167-180.

Ryś, M., Petryk, H., 2018. Gradient crystal plasticity models with a natural length scale in the hardening law. Int. J. Plast. 111, 168-187.

Sandfeld, S., Hochrainer, T., Gumbsch, P., Zaiser, M., 2010. Numerical implementation of a 3D continuum theory of dislocation dynamics and application to micro-bending. Phil. Mag. 90 (27-28), 3697-3728.

Schulz, K., Wagner, L., Wieners, C., 2019. A mesoscale continuum approach of dislocation dynamics and the approximation by a Runge-Kutta discontinuous Galerkin method. Int. J. Plast. 120, 248-261.

Starkey, K., Winther, G., El-Azab, A., 2020. Theoretical development of continuum dislocation dynamics for finite-deformation crystal plasticity at the mesoscale. J. Mech. Phys. Solids 139, 103926.

Stölken, J., Evans, A., 1998. A microbend test method for measuring the plasticity length scale. Acta Mater. 46 (14), 5109-5115.

Sudmanns, M., Stricker, M., Weygand, D., Hochrainer, T., Schulz, K., 2019. Dislocation multiplication by cross-slip and glissile reaction in a dislocation based continuum formulation of crystal plasticity. J. Mech. Phys. Solids 132, 103695.

Tanaka, K., Mura, T., 1981. A dislocation model for fatigue crack initiation. ASME J. Appl. Mech. 48, 97-103.

Taupin, V., Berbenni, S., Fressengeas, C., 2012. Size effects on the hardening of channel-type microstructures: a field dislocation mechanics-based approach. Acta Mater. 60, 664-673.

Van der Giessen, E., Needleman, A., 1995. Discrete dislocation plasticity: A simple planar approach. Modelling Simul. Mater. Sci. Eng. 3, 689-735.

Varadhan, S., Baudoin, A.J., Acharya, A., Fressengeas, C., 2006. Dislocation transport using Galerkin/least squares formulation. Modelling Simul. Mater. Sci. Eng. 14, 1245-1270.

Vivekanandan, V., Lin, P., Winther, G., El-Azab, A., 2021. On the implementation of dislocation reactions in continuum dislocation dynamics modeling of mesoscale plasticity. J. Mech. Phys. Solids 149, 104327.

Voyiadjis, G., Peters, R., 2010. Size effects in nanoindentation: an experimental and analytical study. Acta Mech. 211 (19), 131-153.

Waheed, S., Hao, R., Bhowmik, A., Balint, D.S., Giuliani, F., 2017. A unifying scaling for the Bauschinger effect in highly confined thin films: a discrete dislocation plasticity study. Modelling Simul. Mater. Sci. Eng. 25 (5), 054003.

Wallis, D., Hansen, L.N., Britton, T.B., Wilkinson, A.J., 2016. Geometrically necessary dislocation densities in olivine obtained using high-angular resolution electron backscatter diffraction. Ultramicroscopy 168, 34-45.

Willot, F., 2015. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. C. R. Mec. 343, 232-245.

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

Wulfinghoff, S., Forest, S., Böhlke, T., 2015. Strain gradient plasticity modeling of cyclic behavior of laminate structures. J. Mech. Phys. Solids 79, 1-20.

Xia, S., El-Azab, A., 2015. Computational modelling of mesoscale dislocation patterning and plastic deformation of single crystals. Modelling Simul. Mater. Sci. Eng. 23 (5), 055009.

Zecevic, M., Lebensohn, R.A., Capolungo, L., 2023. Non-local large-strain FFT-based formulation and its application to interface-dominated plasticity of nano-metallic laminates. J. Mech. Phys. Solids 173, 105187.

Zecevic, M., Lebensohn, R.A., Capolungo, L., 2024. Generalized grain boundary constitutive description implemented in a strain-gradient large-strain FFT-based formulation: Application to nano-metallic laminates. J. Mech. Phys. Solids 193, 105859.