Mohamed Jebahi
email: mohamed.jebahi@ensam.eu

Samuel Forest
email: samuel.forest@mines-paristech.fr

Scalar-based strain gradient plasticity theory to model size-dependent kinematic hardening effects

Keywords: Strain gradient plasticity, Size effects, Size-dependent isotropic hardening, Size-dependent kinematic hardening, Equivalent plastic strain, Cumulative plastic strain

A common belief in phenomenological strain gradient plasticity modeling is that including the gradient of scalar variables in the constitutive setting leads to size-dependent isotropic hardening, whereas the gradient of second-order tensors induces size-dependent kinematic hardening. The present paper shows that it is also possible to produce size-dependent kinematic hardening using scalar-based gradient theory. For this purpose, a new model involving the gradient of the equivalent plastic strain is developed and compared with two reference scalar-based and tensor-based theories. Theoretical investigations using simple monotonic loading conditions are first presented to assess the ability of the proposed model to solve some issues related to existing scalar-based gradient theories. Simulations under cyclic loading conditions are then provided to investigate the nature of the resulting hardening. These simulations show that the proposed model is capable of producing size-dependent kinematic hardening effects at more affordable costs, compared to existing tensor-based strain gradient plasticity theories.

Introduction

The size-dependent responses of metallic materials have been observed in numerous small-scale experiments typically ranging in size from a few hundreds of nanometers to a few tens of micrometers [START_REF] Hayashi | Strain hardening in bent copper foils[END_REF][START_REF] Liu | Toward a further understanding of size effects in the torsion of thin metal wires: An experimental and theoretical assessment[END_REF][START_REF] Sarac | Spatial distribution of the net burgers vector density in a deformed single crystal[END_REF][START_REF] Dahlberg | Geometrically necessary dislocation density measurements at a grain boundary due to wedge indentation into an aluminum bicrystal[END_REF]. At small scales, the response of metallic materials is characterized by the statistically stored dislocations (SSDs), the density of which increases with the plastic strain, and the geometrically necessary dislocations (GNDs), which are produced by the heterogeneous plastic flow attributable to the gradient of plastic strains, along with the interactions between them. At these scales, GNDs are generally necessary to accommodate the crystal lattice curvature that arises from a nonuniform plastic deformation. It is commonly believed that material hardening arises from SSDs, and size effects are caused by GNDs which, acting as obstacles to SSDs, can lead to further strengthening and/or hardening. Conventional plasticity theories, which rely on the assumption that the stress at a given point is a function of strain and internal variables at this point only and the GNDs are negligible with respect to the SSDs, cannot predict size effects. Due to lacking internal length scale(s), these theories cannot consider the influence of strain gradients.

To overcome limitations of conventional theories, [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF] has proposed in a pioneering work the first gradient theory of plasticity with a single internal length scale embedded within the conventional J 2 plasticity theory. Since then, a considerable number of enhanced phenomenological and physically-based gradient theories, which are typically termed strain gradient plasticity (SGP) theories, have been developed for single-and poly-crystal structures [START_REF] Gurtin | Gradient single-crystal plasticity with free energy dependent on dislocation densities[END_REF][START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF][START_REF] Dahlberg | Evolution of the length scale in strain gradient plasticity[END_REF][START_REF] El-Naaman | An investigation of back stress formulations under cyclic loading[END_REF][START_REF] Forest | Micromorphic crystal plasticity[END_REF][START_REF] Jebahi | Strain gradient crystal plasticity model based on generalized non-quadratic defect energy and uncoupled dissipation[END_REF][START_REF] Panteghini | A potential for higher-order phenomenological strain gradient plasticity to predict reliable response under non-proportional loading[END_REF][START_REF] Rys | A micromorphic crystal plasticity model with the gradientenhanced incremental hardening law[END_REF]. In the present work, only phenomenological SGP theories are discussed. These theories can be divided into two groups: lower-order and higher-order theories. Lower-order theories preserve the structure of classical plasticity (with conventional stresses, equilibrium equations and boundary conditions), with the difference that the yield condition includes strain gradient terms [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF][START_REF] Mühlhaus | A variational principle for gradient plasticity[END_REF][START_REF] Acharya | Lattice incompatibility and a gradient theory of crystal plasticity[END_REF]. As shown by [START_REF] Niordson | On lower order strain gradient plasticity theories[END_REF], these approaches can lead to likely unacceptable behaviors with an inexplicable localization flow. Higher-order theories substantially deviate from conventional plasticity by the inclusion of new (higher-order) stresses, which are work-conjugate to strain gradients, and additional equilibrium equations and boundary conditions [START_REF] Fleck | A reformulation of strain gradient plasticity[END_REF][START_REF] Forest | Elastoviscoplastic constitutive frameworks for generalized continua[END_REF][START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF][START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF]Fleck et al., 2015;[START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF][START_REF] Panteghini | A potential for higher-order phenomenological strain gradient plasticity to predict reliable response under non-proportional loading[END_REF]. It is now generally accepted that higher-order theories offer important modeling capabilities, making them the most commonly used in the literature. Depending on the nature of the gradient terms, these theories can be classified into scalar-based and tensor-based gradient theories. Some representative phenomenological models in these classes are reviewed hereafter.

For more details, the reader is referred to the interesting review of [START_REF] Voyiadjis | Strain gradient continuum plasticity theories: Theoretical, numerical and experimental investigations[END_REF].

In spirit of the [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF] theory, de Borst and [START_REF] De Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF] have proposed a SGP model based on the gradient of the cumulative plastic strain to investigate localization phenomena in solids.

Implementation aspects of this model have been discussed in de Borst and [START_REF] De Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF] and reviewed in de [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF]. With the purpose of generalizing the classical J 2 flow theory to account for gradient effects at small scales, [START_REF] Fleck | A reformulation of strain gradient plasticity[END_REF] have developed a SGP model based on a generalized definition of a monotonically increasing measure of plastic strain. As formulated, the proposed model does not respect the thermodynamic requirement of non-negative dissipation [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF][START_REF] Gurtin | Thermodynamics applied to gradient theories involving the accumulated plastic strain: The theories of Aifantis and Fleck and Hutchinson and their generalization[END_REF]. A modified version of this model based on the assumption of no gradient contributions to dissipation has been proposed by [START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF] to ensure its thermodynamic consistency.

Such an assumption has been used in subsequent works to develop thermodynamically-acceptable SGP models (e.g., Fleck et al., 2015). Fleck and Willis (2015) have investigated the energetic and dissipative aspects of gradient contributions in SGP theories, with a focus on the phenomenon of elastic gaps. A comprehensive mathematical study of SGP models based on the gradient of cumulative plastic strain is proposed by [START_REF] Fleck | A mathematical basis for strain-gradient plasticity theory. Part I: Scalar plastic multiplier[END_REF]. Most of the aforementioned models involve the gradient of a scalar field variable that increases monotonically. These models generally predict size-dependent isotropic hardening effects under classical loading conditions. Although not very common, it is worth noting here that size-dependent kinematic-like hardening effects can also be produced by such models in some particular cases. Indeed, the gradient contributions, when taken as energetic (or recoverable in the terminology of [START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF]), could be canceled or inverted by applying complex boundary conditions on the additional independent scalar field variable. However, such conditions are not common and, in most cases, they would not have a clear physical interpretation.

The other class of phenomenological SGP models involves the gradient of a tensor field variable.

Among the earliest contributions in this class are [START_REF] Steinmann | Views on multiplicative elastoplasticity and the continuum theory of dislocations[END_REF], where the curl of the plastic strain tensor is introduced in the free energy density as a measure of plastic incompatibility, and [START_REF] Forest | Elastoviscoplastic constitutive frameworks for generalized continua[END_REF] and [START_REF] Gurtin | On a framework for small-deformation viscoplasticity: Free energy, microforces, strain gradients[END_REF][START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF], in which the gradient of the plastic strain tensor or the plastic distortion tensor are included. These contributions have inspired a sizable number of subsequent works proposing tensor-based SGP models. [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF] have proposed a unified framework for tensor-based SGP models that include the gradient of the plastic strain tensor. A mathematical study of such models can be found in Fleck and Willis (2009a). Gurtin and Anand (2005a,b) have developed tensor plastic strain gradient models for small and large deformations. Based on the [START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF] model, [START_REF] Bardella | Modelling the torsion of thin metal wires by distortion gradient plasticity[END_REF] have applied a phenomenological distortion gradient plasticity model to study the effects of the plastic spin on the torsional response of thin metal wires. [START_REF] Panteghini | On the role of higher-order conditions in distortion gradient plasticity[END_REF] have investigated the modeling capabilities of such a model and brought light to its kinematic hardening effects. The gradient terms that are naturally involved in tensor-based gradient models generally lead to size-dependent kinematic hardening effects. For such models to predict size-dependent isotropic hardening effects, it is possible to introduce higher-order dissipative stresses, which must be carefully defined to avoid thermodynamic inconsistencies [START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF]. The existing thermodynamicallyconsistent definitions of such stresses generally lead, under certain non-proportional loading conditions, to elastic gaps whose real existence is questioned [START_REF] Fleck | Strain gradient plasticity under non-proportional loading[END_REF].

In the context of phenomenological SGP, most of the existing scalar-based gradient theories lead to size-dependent isotropic hardening effects under classical boundary conditions. To model size-dependent kinematic hardening effects, it is generally required to use tensor-based gradient theories, which are more expensive than scalar-based ones in terms of computational effort and CPU time, due to the larger number of additional nodal degrees of freedom. The present work aims to show that it is possible to produce the latter effects with the gradient of a scalar field variable. To this end, a new model based on the gradient of the equivalent plastic strain, i.e. a norm of the plastic strain tensor, is developed, in contrast to existing theories relying on the gradient of the cumulative plastic strain. Since kinematic hardening effects are overwhelming in GNDs-dominated plasticity, it is believed that a simple (i.e., by adding a single scalarvalued independent field variable) SGP model incorporating a size-dependent back-stress will be useful for the community. It should be noted that the term "equivalent plastic strain" is sometimes used in the literature to designate the cumulative plastic variable [START_REF] Wulfinghoff | Equivalent plastic strain gradient enhancement of single crystal plasticity: theory and numerics[END_REF]. This term is used here to designate the norm of the total plastic strain tensor.

Although based on similar additional independent field variables, the model proposed in this work, which belongs to the class of flow theories, conceptually differs from the deformation theory models proposed by [START_REF] Aifantis | Section 4.12 -gradient plasticity[END_REF]; [START_REF] Fleck | A reformulation of strain gradient plasticity[END_REF]; [START_REF] Hutchinson | Generalizing J 2 flow theory: Fundamental issues in strain gradient plasticity[END_REF] and [START_REF] Fleck | Strain gradient plasticity under non-proportional loading[END_REF].

In the context of flow theories, to the knowledge of the authors, the gradient of the equivalent plastic strain has only been used to develop non-local damage approaches [START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF]. In the paper of [START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localisation[END_REF], a variable named "effective plastic strain" is introduced for some SGP versions. It is not clear whether this variable corresponds to the equivalent or the cumulative plastic strain. However, the expression of the associated consistency condition suggests that it corresponds to the cumulative one.

Indeed, as will be seen later, the use of the equivalent plastic strain modifies the expression of such a condition. Although the "equivalent plastic strain" is introduced in the paper Fleck and Willis (2015), the proposed incremental model used to obtain the paper results involves the gradient of the plastic strain tensor.

After the present introduction, the paper is organized as follows. For comparison purposes, section 2 briefly recalls the main features of scalar-based and tensor-based theories. Section 3 presents the new equivalent plastic strain gradient model. A theoretical investigation of the proposed model is given in section 4 to assess its ability to solve some issues related to conventional scalar-based gradient theories. Section 5 discusses the numerical implementation of the developed model. Numerical results obtained under monotonic and cyclic loading conditions are provided in section 6. Comparison between the proposed model and existing scalar-based and tensor-based theories is also given in this section. Finally, section 7 presents some concluding remarks.

Scalar-based versus tensor-based strain gradient plasticity

The present section aims to review the main features of scalar-based and tensor-based SGP theories. For this purpose, the [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF] theory, which implies the gradient of cumulative plastic strain, and a simplified version of the [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF] theory, which relies on the gradient of the plastic strain tensor, are briefly presented.

Scalar-based strain gradient plasticity

The pioneering scalar-based SGP theory of [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF] is derived hereafter using the generalized power density of internal forces [START_REF] Forest | Some links between recent gradient thermo-elasto-plasticity theories and the thermomechanics of generalized continua[END_REF][START_REF] Forest | Formulations of strain gradient plasticity[END_REF]. In the framework of small deformation, this density can be written as:

p (i) = σ : ε + a ṗ + b • ∇ ṗ (1)
where ε is the linearlized strain tensor which is additively decomposed into elastic ε e and plastic ε p parts, p and ∇p are respectively the cumulative plastic strain and its first gradient, σ is the Cauchy stress tensor, and a and b are generalized stresses acting on respectively p and ∇p. Application of the generalized principle of virtual power [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, première partie : théorie du second gradient[END_REF][START_REF] Forest | Micromorphic Approach for Gradient Elasticity, Viscoplasticity, and Damage[END_REF] results in an additional balance equation, complementing the classical balance equations. In static case and neglecting body forces, these equations can be written:

div σ = 0, a = div b (2)
The entropy imbalance is expressed in local form as:

p (i) -ψ ≥ 0 (3)
where ψ is the free energy potential per unit volume, assumed to be a function of ε e , p and ∇p:

ψ = ψ(ε e , p, ∇p) (4)
Using this function, the Clausius-Duhem inequality becomes:

(σ - ∂ψ ∂ε e ) : εe + σ : εp + (a - ∂ψ ∂p ) ṗ + (b - ∂ψ ∂∇p ) • ∇ ṗ ≥ 0 (5)
with the state laws taken as:

σ = ∂ψ ∂ε e , R = ∂ψ ∂p , b = ∂ψ ∂∇p (6)
The residual dissipation is then:

σ : εp + (a -R) ṗ ≥ 0 (7)
Choosing a simple quadratic free energy potential:

ψ(ε e , p, ∇p) = 1 2 ε e : C : ε e + 1 2 h p 2 + 1 2 c ∇p • ∇p (8)
the state laws become:

σ = C : ε e , R = h p, b = c ∇p (9)
where C is the fourth-order tensor of elastic moduli, h is the usual hardening modulus and c is an additional material parameter (unit MPa • mm 2 = N), as in [START_REF] Aifantis | On the Microstructural Origin of Certain Inelastic Models[END_REF]. The latter parameter is assumed to be positive c = H l 2 , with H a generalized hardening modulus and l an internal length scale.

The yield function is taken as:

f (σ) = J 2 (σ) + a -R -R 0 with J 2 (σ) = 3 2 σ : σ (10)
where R 0 is the initial yield stress and σ is the deviatoric part of σ. Under plastic loading, this gives

J 2 (σ) = σ eq = R 0 + R -a = R 0 + R -div b = R 0 + h p -c ∇ 2 p (11)
where the last term is the Laplacian of the cumulative plastic strain field. As for classical plasticity theory, the flow rule writes:

εp = ṗ ∂f ∂σ ( 12
)
with the plastic multiplier ṗ = 2 3 εp : εp . Time integration of this parameter, with a given initial condition, provides the cumulative plastic strain variable p.

Tensor-based strain gradient plasticity

The main features of tensor-based SGP theories are reviewed hereafter based on the contribution of [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF]. Although a quite general theoretical framework has been proposed in [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF], covering a large range of strain gradient plasticity theories, only a simplified version of this framework is presented in this subsection. In particular higher-order dissipation contributions are not considered. For comparison purposes, this version is derived following the same methodology as the one used in the previous subsection. Following the thermodynamic framework proposed in [START_REF] Forest | Formulations of strain gradient plasticity[END_REF], the generalized power density of internal forces can be written as:

p (i) = σ : ε + q : εp + M . . .∇ εp (13)
where q and M are second-order and third-order generalized stresses, work-conjugate to ε p and ∇ε p , respectively. Application of the generalized principle of virtual power [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, première partie : théorie du second gradient[END_REF][START_REF] Forest | Micromorphic Approach for Gradient Elasticity, Viscoplasticity, and Damage[END_REF] results

in the following balance equations (static case, without body forces):

div σ = 0, q = div M (14)
The entropy imbalance in local form is given by:

p (i) -ψ ≥ 0 (15)
with ψ assumed to be a function of ε e , p and ∇ε p :

ψ = ψ(ε e , p, ∇ε p ) (16)
The Clausius-Duhem inequality can be written as:

(σ -∂ψ ∂ε e ) : εe + (σ + q) : εp -

∂ψ ∂p ṗ + (M - ∂ψ ∂∇ε p ) . . .∇ εp ≥ 0 (17)
with the state laws taken as:

σ = ∂ψ ∂ε e , R = ∂ψ ∂p , M = ∂ψ ∂∇ε p (18)
The residual dissipation is:

(σ + q) : εp -R ṗ ≥ 0 ⇔ (σ -X) : εp -R ṗ ≥ 0 (19)
where the back-stress is defined by:

X = -q = -div M (20)
Using a simple quadratic free energy potential:

ψ(ε e , p, ∇ε p ) = 1 2 ε e : C : ε e + 1 2 h p 2 + 1 2 c ∇ε p . . .∇ε p (21)
the state laws become:

σ = C : ε e , R = h p, M = c ∇ε p (22)
Using the above expression of M, the back-stress X can be obtained:

X = -c ∇ 2 ε p (23)
which depends on the Laplacian of the plastic strain tensor. The yield function is taken as:

f (σ, X) = J 2 (σ -X) -R -R 0 with J 2 (σ -X) = 3 2 σ -X : σ -X (24)
The flow rule writes:

εp = ṗ ∂f ∂σ ( 25
)
3 Equivalent plastic strain gradient theory

The present section details a new scalar-based strain gradient plasticity theory, in which the gradient of the cumulative plastic strain is replaced by the gradient of the equivalent plastic strain defined by:

ε eq = 2 3 ε p : ε p (26)
It must be noted that this norm of the plastic strain tensors differs from the cumulative plastic strain p defined in the previous section by the time integration of the relation :

ṗ = 2 3 εp : εp (27)
Considering the thermodynamic framework proposed in [START_REF] Forest | Formulations of strain gradient plasticity[END_REF], the generalized power density of internal forces can be written as:

p (i) = σ : ε + a εeq + b • ∇ εeq (28)
with a and b generalized stresses acting on ε eq and ∇ε eq , respectively. Application of the generalized principle of virtual power results in the same generalized balance of momentum equations as those obtained using the Aifantis theory (static case, without body forces):

div σ = 0, a = div b (29)
In this section, the free energy potential is assumed to be a function of ε e , p and ∇ε eq1 :

ψ = ψ(ε e , p, ∇ε eq ) (30)
This leads to the following Clausius-Duhem inequality:

(σ - ∂ψ ∂ε e ) : εe + σ : εp + a εeq - ∂ψ ∂p ṗ + (b - ∂ψ ∂∇ε eq ) • ∇ εeq ≥ 0 (31)
with the state laws taken as:

σ = ∂ψ ∂ε e , R = ∂ψ ∂p , b = ∂ψ ∂∇ε eq (32)
The residual dissipation is:

d = σ : εp + a εeq -R ṗ ≥ 0 (33)
The time derivative of the equivalent plastic strain is:

εeq = 2 3 ε p ε eq : εp (34)
Using this expression, the residual dissipation can be rewritten:

d = (σ -X) : εp -R ṗ ≥ 0 (35)
where a back-stress X arises:

X = - 2 3 a ε p ε eq = - 2 3 div b ε p ε eq (36)
Consider a simple quadratic free energy potential:

ψ(ε e , p, ∇ε eq ) = 1 2 ε e : C : ε e + 1 2 h p 2 + 1 2 c ∇ε eq • ∇ε eq (37)
from which the state laws are derived:

σ = C : ε e , R = h p, b = c ∇ε eq (38)
In this case, the back-stress X becomes:

X = - 2 3 c ∇ 2 ε eq ε p ε eq (39)
which depends on the Laplacian of the equivalent plastic strain and on the direction of the plastic strain tensor. In the case of zero plastic strain, X is indeterminate since both the numerator and the denominator vanish. This special case is treated numerically as it will be explained later. The yield function writes:

f (σ, X) = J 2 (σ -X) -R -R 0 with J 2 (σ -X) = 3 2 σ -X : σ -X (40)
The same flow rule as in the previous subsections is adopted:

εp = ṗ ∂f ∂σ (41)
The plastic multiplier ṗ = 2 3 εp : εp is linked to εeq by:

εeq -(n ε • n σ ) ṗ = 0 (42)
with A conceptual problem of gradient plasticity theories including gradient of scalar plastic field variables has been noticed by [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF]. As scalar variables generally contain no information on the direction of the plastic flow, [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF] have shown that an arbitrary small perturbation in the boundary conditions can determine the direction of the plastic deformation in many scalar-based gradient formulations. In this sense, the solution is not stable with respect to the boundary conditions, which seems to be physically-unacceptable. The present subsection aims to investigate this problem in the context of the proposed equivalent plastic strain gradient theory. To this end, a simple one-dimensional elasto-visco-plastic model, inspired by that proposed by [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF], is considered. As in [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF], in the present investigation, plasticity is assumed to be governed by a simple power low flow rule.

n ε = 2 3 ε p ε eq and n σ = 3 2 σ -X J 2 (σ -X)
Using Aifantis theory, the plastic response of the one-dimensional model under tension-compression loading is given by:

εp = ε0 sign(σ) f σ d n = ε0 sign(σ) |σ| + c p ,11 -R -R 0 σ d n ( 43
)
with ε0 reference strain rate, n rate sensitivity parameter, σ d drag stress and " • " denoting the Macaulay brackets defined as

• = (• + |•|) /2 = max {•, 0}.
A conceptual difficulty arises when the yield criterion is satisfied (f = 0) at a given point and at the same time the Cauchy stress is zero. Such a situation occurs, for example, when the higher-order term is sufficiently large (σ = 0 and c p ,11 = R + R 0 ). In this case, the plastic strain is triggered, but the flow rule (43) cannot provide the direction of εp . Indeed, using Aifantis theory, the direction of εp is given by the sign of σ. Therefore, as it will be seen later, an arbitrary small perturbation in σ can determine the direction of the total plastic strain. On the contrary, using the present formulation, the flow rule (43) becomes:

εp = ε0 sign(σ + c ε p ,11 ) f σ d n = ε0 sign(σ + c ε p ,11 ) σ + c ε p ,11 -R -R 0 σ d n (44)
For nonzero yield stress, satisfaction of the yield criterion means that σ + c ε p ,11 = 0. The direction of εp is then always well defined.

In order to further investigate the aforementioned conceptual problem, the model of figure 1, which is inspired by that proposed by [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF], is considered in the following. A one-dimensional bar of length a is subjected to initial conditions ε p 0 and an additional small perturbation F on its right end. The direction of the perturbation is determined by the sign of . For simplicity, the yield stress is set to zero (σ y = R + R 0 = 0) and the rate sensitivity parameter to one (n = 1). The model is then reduced to a visco-elastic gradient model, which is sufficient to analyze the aforementioned conceptual problem.

The initial conditions are assumed to be:

ε p (x, t = 0) = p (x, t = 0) = ε p 0 (x) = F 4 A c (x -a) 2 (45)
with A the cross sectional area of the bar. These initial conditions allow for a simple analytical solution of the problem. The boundary conditions applied on the ends of the bar are assumed to be:

ε p ,1 (x = 0, t) = - F a 2 A c and ε p ,1 (x = a, t) = 0 (46)
Using Aifantis theory, the evolution of the plastic strain is given by:

εp = ε0 sign( ) | F | + A c p ,11 A σ d (47)
A possible solution of this equation, in the limit case of → 0, is given by:

ε p = sign( ) F ε0 2 A σ d t + ε p 0 (x) and p = F ε0 2 A σ d t + ε p 0 (x) (48)
For sufficiently large time, the sign of the strain in the bar, which can be approximated by ε p , is determined by the sign of . This means that the prescribed perturbation F , even though negligible, can determine whether the bar is stretched or compressed. On the contrary, using the present theory, the evolution of the plastic strain is governed by:

εp = ε0 F + A c ε p ,11 A σ d (49)
whose solution, in the limit case of → 0, can be written:

ε p = p = F ε0 2 A σ d t + ε p 0 (x) (50)
This solution is independent of the sign of the introduced perturbation. Although the present theory is based on the gradient of a scalar field variable, it takes into account the direction of the plastic strain, making the solution stable with respect to the boundary conditions. Indeed, the back-stress resulting from the higher-order term is a function of the direction of the plastic strain tensor.

Constrained plasticity in a sheared layer

In this subsection, the equivalent plastic strain gradient model is applied to a constrained elasto-plastic strip subjected to monotonic shear loading. Similar problems have been investigated by [START_REF] Forest | Formulations of strain gradient plasticity[END_REF] and [START_REF] Forest | Questioning size effects as predicted by strain gradient plasticity[END_REF] in the case of Aifantis theory. A 2D strip of length 2a and width 2w is considered (figure 2). The origin O of the coordinate system is located at the center of the strip.

To model the infinite length of this strip in e 1 -direction, periodic boundary conditions are applied on its left and right edges:

u i (a, x 2 , t) = u i (-a, x 2 , t) , for i = 1, 2 ε eq (a, x 2 , t) = ε eq (-a, x 2 , t) (51)
The top and bottom edges of the strip are subjected to opposite displacements in e 1 -direction and zero displacements in e 2 -direction:

u 1 (x 1 , w, t) = -u 1 (x 1 , -w, t) = w γ and u 2 (x 1 , w, t) = u 2 (x 1 , -w, t) = 0 (52)
where γ denotes the mean simple shear applied in e 1 -direction. In addition to these classical conditions, these two edges are also subjected to special boundary conditions, namely vanishing plastic strain, to

ε eq (x 1 , w, t) = ε eq (x 1 , -w, t) = 0 (53)
Assuming plane stain conditions, the displacement field in the strip can be written:

u 1 (x 2 ) = γ x 2 + u (x 2 ) and u 2 = 0 (54)
where u, which describes the fluctuation from the homogeneous solution, is the main unknown of the boundary value problem. The strain and stress tensors are written as:

ε = ε 12 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) and σ = σ 12 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) (55)
where ε 12 (x 2 ) = 1 2 [γ + u ,2 (x 2 )] depends only on x 2 , and σ 12 is homogeneous as a consequence of the balance of momentum equations. After yielding, the plastic strain tensor takes the form:

ε p = ε p 12 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) (56)
The equivalent plastic strain and its gradient can be expressed as:

ε eq = 2 √ 3 |ε p 12 | , ∇ε eq = ε eq,2 e 2 (57)
The yield condition reads:

f (σ, X, R) = √ 3 |σ 12 -X 12 | -R 0 -R, with X 12 = - 2 3 c ε p 12,22
In the monotonic case with positive values of the variables, it gives:

σ 12 = R 0 + R √ 3 - 2 3 c ε p 12,22 (58)
This equation can be solved and provides a parabolic plastic profile when R = 0 and a cosh profile for linear hardening R = h p. The same equation can be obtained using Aifantis theory [START_REF] Forest | Formulations of strain gradient plasticity[END_REF]. Under monotonic conditions, the proposed theory is equivalent to that of Aifantis. As will be seen later, differences will arise under cyclic loading conditions. Although solution of ( 58) is given in the literature for specific boundary conditions [START_REF] Forest | Formulations of strain gradient plasticity[END_REF], it will be reproduced hereafter with further specialization considering the boundary conditions of this subsection. The specialized solution will be used later to validate the implementation of the proposed model. Assuming linear hardening R = h p, with p = ε eq under monotonic conditions, equation (58) yields:

√ 3 σ 12 -R 0 h = ε eq - c h ε eq,22 (59)
Since σ 12 is uniform in the strip, differentiation of this equation with respect to x 2 leads to:

ε eq,2 -λ 2 ε eq,222 = 0, with λ = c h (60)
λ can be interpreted as an effective length scale. Solution of this equation has the following form:

ε eq = A cosh x 2 λ + B ( 61
)
where A and B are integration constants. The displacement fluctuation u in the strip can be obtained from the Hooke constitutive law:

σ 12 = µ γ + u ,2 (x 2 ) - √ 3 ε eq (x 2 ) (62)
which, considering (59), yields:

u (x 2 ) = R 0 √ 3 µ -γ + h + 3 µ √ 3 µ B x 2 + √ 3 A λ sinh x 2 λ ( 63
)
where µ represents the shear modulus. The arbitrary translation for u is set to zero. Constants A and B can be determined based on the boundary conditions at top and bottom edges u (±w) = 0 and ε eq (±w) = 0:

A = R 0 √ 3 µ -γ h + 3 µ √ 3 µ cosh w λ - √ 3λ w sinh w λ B = -A cosh w λ (64)

Bending of thin foil

It has been shown in the literature that scalar-based gradient theories can lead to non-smooth distribution of plastic strains under certain loading conditions. Using bending conditions up to the limit load, Poh et al. ( 2011) have pointed out that the plastic profile obtained using cumulative plastic strain gradient theory is not smooth at the neutral axis. To overcome this difficulty, these authors have proposed to use the gradient of the full plastic strain tensor instead of the gradient of the cumulative plastic strain. The present section aims to investigate this point using the proposed equivalent plastic strain gradient theory.

To this end, the bending example considered by [START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF] will be reproduced using this theory.

In this example, a foil of width 2w is bent around e 3 -axis to a curvature κ (figure 3). For simplicity, κ is supposed to increase monotonically as a function of time t and the mechanical state is assumed to evolve under plane strain conditions. Assuming, in addition, that the elastic and plastic responses are incompressible, the displacement field can be approximated, in the case of small deformation, as:

u = κx 1 x 2 e 1 - κ 2 x 2 1 + x 2 2 e 2 (65)
The strain tensor ε is given by: with ε 11 = κx 2 . Due to the symmetry with respect to the neutral axis, the following analysis is restricted to the upper half of the foil x 2 ≥ 0; the solution for the lower half can easily be obtained using symmetry or anti-symmetry arguments. In the regions where the yield condition is fulfilled, the plastic strain tensor can be written as:

ε = ε 11 (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) (66)
ε p = ε p 11 (e 1 ⊗ e 1 -e 2 ⊗ e 2 ) ( 67
)
where ε p 11 depends only on x 2 . The equivalent plastic strain and its gradient are given by:

ε eq = 2 √ 3 |ε p 11 | , ∇ε eq = ε eq,2 e 2 (68)
The yield criterion reads:

f (σ, X, R) = √ 3 |σ 11 -X 11 | -R 0 -R, with X 11 = - 2 3 c ε p 11,22 (69)
In the present bending problem, the stress tensor is not uniform in the studied foil and the plastic region evolves with loading. The problem will then be solved incrementally using the rate form of (69). Assuming positive variables in the upper half of the studied foil and replacing ε p 11 by its expression in terms of ε eq (68), the rate form of equation ( 69) in the case of linear hardening (R = h p) can be written:

σ11 = h √ 3 εeq - c √ 3 εeq,22 (70)
Applying the Hooke constitutive law, σ11 can be written:

σ11 = µ 2 κ x 2 - √ 3 εeq (71)
Substituting ( 71) into (70), the following equation can be obtained:

εeq -λ 2 εeq,22 = 2 √ 3 µ h + 3µ κ x 2 with λ = c h + 3 µ (72)
This equation, which is obtained assuming monotonic bending conditions, is the same as the one obtained using Aifantis theory [START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localisation[END_REF][START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF]. The general solution of such an equation is given by:

ε eq (x 2 , t) = 2 √ 3µ h + 3µ κ (t) x 2 + C 1 (t) cosh( x 2 λ ) + C 2 (t) sinh( x 2 λ ) + C 3 (x 2 ) ( 73
)
where C 1 , C 2 and C 3 are integration functions to be determined based on the boundary conditions.

Considering the general case with elastic and plastic regions, the position x 2 of the elasto-plastic boundary depends of time (x 2 = d(t)). The function d is also to be determined based on the boundary conditions.

To identify C 1 , C 2 , C 3 and d, the following boundary conditions are considered:

-Generalized traction imposed on the top edge: ε eq,2 (w, t) = 0

-Continuity of ε eq at x 2 = d (t): ε eq (d (t) , t) = 0 -Continuity of classical traction at x 2 = d (t): R0 c -2 √ 3µ c κ (t) d (t) = ε eq,22 (d (t))
-Continuity of generalized traction at x 2 = d (t): ε eq,2 (d (t) , t) = 0

These boundary conditions are similar to those applied by [START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localisation[END_REF] and [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF] in the case of cumulative plastic strain gradient theory. Application of the above boundary conditions leads to the following expression of ε eq :

εeq (x 2 , t) = 2 √ 3 µ 3µ + h κ (t) cosh w-d(t) λ x 2 cosh w -d (t) λ -λ sinh x 2 -d (t) λ -d (t) cosh w -x 2 λ - R 0 3µ + h 1 cosh w-d(t) λ cosh w -d (t) λ -cosh w -x 2 λ (74)
with the function d verifying [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF]:

2 √ 3 µ λ κ (t) 3µ + h 1 -cosh w -d (t) λ - d (t) λ sinh w -d (t) λ + R 0 3µ + h sinh w -d (t) λ = 0 (75)
There exists no analytical solution for equation ( 75). However, as demonstrated by [START_REF] Engelen | An evaluation of higher-order plasticity theories for predicting size effects and localisation[END_REF] and [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF], numerical resolution of this equation shows a growth in the plastic region which is faster than when using classical plasticity theory. Consequently, the elasto-plastic boundary meets the neutral axis at a finite time. Beyond this time, the equivalent plastic strain ε eq reduces to:

εeq (x 2 , t) = 2 √ 3 µ 3µ + h κ (t) cosh w λ x 2 cosh w λ -λ sinh x 2 λ - R 0 3µ + h 1 cosh w λ cosh w λ -cosh w -x 2 λ (76)
It can be verified that ε eq is zero at the neutral axis, but not ε eq,2 which has a finite value. This means that generalized traction exists at the neutral axis but it does not expend work, as ε eq is zero at this axis.

Consequently, the proposed equivalent plastic strain gradient model does not overcome the problem of non-smoothness of plastic strain at the neutral axis, which results in nonzero generalized traction at this axis. However, the obtained traction does not contribute to the system work.

To simplify the analytical derivation of the bending solution, the above bending problem is often solved considering half of the studied foil in the limit case for which the elasto-plastic boundary coincides with the neutral axis. In this case, boundary conditions has to be imposed on this axis to ensure well-posedness of the problem. By imposing zero generalized traction on the neutral axis, smooth plastic strain distribution can be obtained [START_REF] Idiart | Size effects in the bending of thin foils[END_REF][START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF]. However, this choice of boundary conditions can lead to appearance of an internal layer originating from the neutral axis in which the stress state presents an opposite sign with respect to the state in the rest of the foil half. This is questionable as nothing special happens at the neutral axis to change the stress sign. Appearance of such likely-unphysical layer seems to be an apparent problem due to the arbitrary imposition of generalized traction-free boundary conditions on the neutral axis. The general solution, taken in the limit case of d(t) = 0, suggests that essential boundary conditions have to be imposed on the neutral axis. Although these conditions lead to non-smooth plastic strain distributions with nonzero generalized traction at the neutral axis, the obtained traction does not contribute to the system work. Several scalar-based theories leading to nonsmooth plastic strain distributions under certain loading conditions can be found in the literature. Using a minimal gradient-enhancement approach, [START_REF] Stupkiewicz | A minimal gradient-enhancement of the classical continuum theory of crystal plasticity. part ii: Size effects[END_REF] have obtained a non-smooth plastic strain profiles within a constrained strip subjected to simple shear loading. This approach is based on the so-called implicit-gradient regularization [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF][START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF], allowing for reproducing correctly non-smooth profiles. Very recently, [START_REF] Rys | A micromorphic crystal plasticity model with the gradientenhanced incremental hardening law[END_REF] have proposed a new approach combining the micromorphic model [START_REF] Forest | Micromorphic Approach for Gradient Elasticity, Viscoplasticity, and Damage[END_REF] and the minimal gradient-enhanced model [START_REF] Stupkiewicz | A minimal gradient-enhancement of the classical continuum theory of crystal plasticity. part ii: Size effects[END_REF]. This approach also predicts non-smooth plastic strain distributions under constrained shear conditions. Further research effort is needed to investigate the physical nature of this kind of distributions. For the sake of clarity, vector notation will be adopted in this section. In the following, bold symbols denote the vector representation of tensor quantities. The finite element implementation of the equivalent plastic strain gradient model involves the computation of the displacement and equivalent plastic strain field variables (u and ε eq ) that verify the macroscopic balance equations ( 29) 1 and the yield condition (40) in a weak sense. It is worth emphasizing that the equivalent plastic strain field ε eq is taken as independent field variable. The weak forms of (29) 1 and ( 40) can be formulated using the principle of virtual power. Let δ u be a virtual velocity field kinematically admissible to zero on the portion S u of the domain boundary where displacement is given, and let δ εeq be a virtual equivalent plastic strain rate field kinematically admissible to zero on the portion S ε of the domain boundary on which the equivalent plastic strain is given. The weak forms associated with (29) 1 and ( 40) can be written:

G u = V (∇δ u) T • σ dv - St δ uT • t ds = 0 G ε = V δ εeq [J 2 (σ -X) -R -R 0 ] dv = 0 (77)
where S t is the boundary portion on which classical traction force t is imposed. To numerically solve these weak forms, a User-ELement (UEL) subroutine is implemented within the commercial finite element package ABAQUS/Standard. A compact description of the UEL algorithms is given hereafter in twodimensional (2D) plane strain case. Extension to three-dimensional case can easily be derived.

In the above integral equations, there appear at most first-order derivatives of the displacements and second-order derivatives of the equivalent plastic strain. Therefore, the discretization of the displacement field u requires C 0continuous interpolation functions and the discretization of the equivalent plastic strain field ε eq requires C 1continuous interpolation functions. In the context of scalar-based gradient theories, the C 1continuity requirement of the additional scalar field variable has been discussed by de [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF]. In the present theory, this requirement is all the more necessary, due to the nonlinear dependence of the yield condition on the Laplacian term (∇ 2 ε eq ). Because of this nonlinear dependence, it is not easy to eliminate the Laplacian term from the associated weak form by a simple application of the divergence theorem. To ensure C 1continuity of ε eq , the C 1continuous element formulation, which has been initiated by de Borst and [START_REF] De Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF], is applied. Figure 4 presents the finite element used in this formulation. This element employs eight-noded (quadratic) interpolation of u , four-noded cubic Hermitian interpolation of ε eq and 2 × 2 Gauss integration. Such an element has been tested in the context of Aifantis theory by de [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF]. According to these authors, it is the most reliable among several tested C 1continuous elements, with a satisfactory convergence behavior and an exact fulfillment of the yield condition at the integration points. However, because of the presence of the cross-derivative degrees of freedom ε eq,xy , this element is required to be parallel to the global reference system, restricting its use to simple geometrical problems.

Within a finite element, the displacements and equivalent plastic strain fields can be approximated by:

u = N u • U e and ε eq = N T ε • E e eq
N u is the interpolation matrix associated with the displacement field, N ε is the interpolation vector associated with the equivalent plastic strain field, U e and E e eq are respectively the element nodal variables of displacement and equivalent plastic strain. The strain vector ε can be written as:

ε = B u • U e
with B u the gradient matrix associated with the displacement field. The gradient and the Laplacian of the equivalent plastic strain can be approximated as:

∇ε eq = B T ε • E e eq and ∇ 2 ε eq = P T ε • E e eq
where B T ε is the gradient matrix associated with the equivalent plastic strain field and P T ε contains the Laplacian terms of the C 1continuous interpolation functions. Using these approximations, the above weak forms can be written within an element as:

G e u = δ Ue T • V e B T u • σ dv - e t N T u • t ds = 0 G e ε = δ Ėe eq T • V e [J 2 (σ -X) -R -R 0 ] N ε dv = 0 (78)
The principle of virtual power implies that G e u and G e ε are zero for any virtual variations of the element nodal variables δ Ue and δ Ėe eq , which results in:

R e u = V e B T u • σ dv - S e t N T u • t ds = 0 R e ε = V e [J 2 (σ -X) -R -R 0 ] N ε dv = 0 (79)
These equations are linearized with respect to the variations of the element nodal variables U e and E e eq , which leads to an elementary system of linear equations: The global system of linear equations can be obtained by assembling all the elementary systems associated with the overall finite elements. This system is solved by means of a Newton-Raphson iterative solution scheme for the overall increments of the displacement and equivalent plastic strain fields. At each iteration, updated values of these increments are obtained and used to numerically solve the constitutive equations at the Gauss points. Using Aifantis theory, where the cumulative plastic strain increments ∆p at the Gauss points can be determined from the nodal values, these equations can easily be solved using, for example, the procedure detailed in de [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF]. However, this procedure cannot directly be applied to the present theory, where only equivalent plastic strain increments ∆ε eq at the Gauss points can be determined from the nodal values. An additional step to calculate ∆p from ∆ε eq is needed. This can be done by solving the nonlinear equation given by ( 42), which can be expressed in discretized form as:

∆ε eq -n T ε • n σ ∆p = 0 (80)
This equation is solved using implicit backward Euler integration scheme. Note that n ε , which provides the direction of the back-stress, is indeterminate at zero plastic strain. This particular case is treated numerically as follows. Assuming ε eq = 0,

if, in addition, ∆ε eq = 0 (elastic regime), there is no need to calculate n ε and the back-stress is given by:

X n+1 = X n , with X 0 = 0
if ∆ε eq = 0, ∆ε p is first updated and used to calculate n ε :

n ε = 2 3 ε p ε eq ≈ 2 3
∆ε p ∆ε eq

The integration procedure of the constitutive equations is detailed in algorithm 1.

Algorithm 1: Integration of the equivalent plastic strain gradient constitutive equations Input: ∆ε, ε p , ε eq , ∆ε eq , ∇ 2 ε eq , p, σ,

X Compute ∆σ trial = C • ∆ε, σ trial = σ + ∆σ trial , f (σ trial , X) = J 2 (σ trial -X) -R(p) -R 0 if f (σ trial , X) < 0 then Set σ = σ trial ; Set h = E; else
Initialize ∆P = 0, dp = 0, Flag = 0; while (dp > or Flag = 0) do Update Flag = 1;

Compute n ε and n σ ;

Compute ∆P + = dp, ε p + = dp n σ ; Compute σ = σ trial -dp C • n σ , X = -c ∇ 2 (ε eq ) n ε ; Compute L = ∆ε eq -n T ε • n σ ∆p; Compute K p = ∂(dL) ∂(dp) ; Compute dp = -L/K p ; end end 6 Numerical results
After validation of the model implementation, the present section aims at investigating the nature of the hardening effects produced by the proposed equivalent plastic strain gradient model. To this end, results obtained by this model under cyclic loading conditions will be compared with those obtained using Aifantis and Gudmundson theories.

Validation of the model implementation

To validate the implementation of the equivalent plastic strain gradient model, it is applied to simulate the simple shear problem of subsection 4.2, using a special choice of material parameters oriented towards 1) and a prescribed mean simple shear of γ = 0.01. The geometrical model involved in this validation study is discretized using 50 × 50 elements. The boundary conditions are the same as in Fig. 2 and as described in Section 4.2. Figure 5 presents the obtained results in terms of shear stress versus imposed shear and distribution of equivalent plastic strain along e 2 -direction.

The associated analytical results are also presented in this figure. The numerical results are in perfect agreement with the analytical ones, which validates the model implementation.

To study the capability of the proposed C 1continuous element implementation to capture the nonsmooth distribution of equivalent plastic strain at the neutral axis, it is applied to simulate the bending problem of subsection 4.3. To numerically reproduce bending loading, four-point flexural conditions are applied on a foil of length 8 a and width 2 w (a and w are given in table 1). Figure 6 presents the associated results in terms of equivalent plastic strain profiles along e 2 -direction at different imposed curvatures.

These results are obtained with a discretization of 80×20 elements and using the material parameters given in table 1, except for c which is set as c = 0.25 N in this example. For κ = 1.17×10 -2 mm -1 , relatively good results, which compare favorably with the theoretical ones, are obtained. The small differences between the theoretical and numerical plastic strain distributions could be due to the numerical approximation of the bending loading conditions. For κ = 3.67 × 10 -2 mm -1 , which corresponds to the curvature at which the theoretical plasticity reaches the neutral axis, far from the foil center, the plastic strain distributions also compare favorably with the theoretical ones. On the contrary, poor results with non-smooth plastic strain profiles are obtained in the vicinity of the neutral axis. Mesh effects on these results are presented in figure 7, which is obtained using two mesh refinements (80 × 20 elements for coarse mesh and 200 × 50 elements for fine mesh). The smoothness of the plastic strain distributions in the vicinity of the neutral axis can be enhanced by using finer mesh. However, an excessively fine mesh is required to faithfully reproduce the smooth distributions as obtained theoretically. Beyond the limit curvature, convergence problems are encountered. Actually, the proposed model implementation, which is based on C 1continuous element formulation, is not adequate to capture non-smooth plastic strain distributions. Indeed, this formulation is based on the assumption of continuous plastic strains and their gradients. An alternative approach to capture non-smooth plastic strain distributions is to use a micromorphic approach for which the plastic strain is non-smooth but the micromorphic variable remains smooth [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF][START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF]. This is however a different model. In SGP, one could use discontinuous Galerkin methods (Di Pietro and Ern, 2012) or hybrid higher-order (HHO) methods [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF] which allow for modeling discontinuities.

x 2 /w Theoretical -κ = 1.17 × 10 -2 mm -1 Numerical -κ = 1.17 × 10 -2 mm -1 Theoretical -κ = 3.67 × 10 -2 mm -1 Numerical -κ = 3.67 × 10 -2 mm -1

Cyclic shear response of a constrained layer: Equivalent SGP versus Cumulative SGP

To investigate the nature of the higher-order hardening (i.e., coming from the Laplacian term in the yield condition), the present subsection provides a comparison between the proposed and the Aifantis models under cyclic simple shear loading. The shear model involved in the present simulations is similar to that presented in figures 2, with the monotonic loading replaced by a complete cycle of loading-unloading. For a better illustration of the higher-order hardening, the first-order one (i.e., coming from plastic strains) is neglected in the following (h = 0). The other material parameters are the same as those presented in table 1, except for c which is taken as variable parameter. The value of this parameter (c = H l 2 ) is adjusted via the internal length scale l, assuming fixed higher-order hardening parameter H = 10 GPa.

To better analyze the nature of the obtained hardening, the results of the present theory are compared with those of the Aifantis theory. The latter theory has previously been implemented within a UEL subroutine in ABAQUS/Standard, using the procedure detailed in de [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF]. This procedure is also based on C 1continuous element formulation. Therefore, the only differences between the two implementations are the choice of the independent plastic variable (equivalent or cumulative plastic strain) and the consequence of this choice on the integration of the constitutive laws. This makes the comparison between the two theories more reliable.

Figure 8 compares the cyclic shear results of these theories for different values of l. For both theories, the internal length scale l contributes to hardening which increases with increasing l. Using Aifantis theory, isotropic hardening is obtained. On the contrary, although based on the gradient of a scalar field variable, the present theory leads to kinematic hardening with classical Bauschinger effects. Figure 9 shows the profiles of the additional independent scalar field variables involved in the present and the Aifantis theories (i.e., equivalent plastic strain and cumulative plastic strain, respectively) at the loading stages depicted in figure 9a. It can be noted that, contrary to what is presented in figure 5b, these profiles display a parabolic shape which is a consequence of ignoring the first-order hardening (h = 0). In figure 9b, the equivalent plastic strain, on which the present theory relies, presents the same distribution at stages 1 and 3 and vanishes at stage 2. This variable increases within the studied strip up to stage 1.

Under reverse shear loading, it decreases until stage 2 where it vanishes before again increasing between stages 2 and 3. The resulting back-stress exhibits a similar non-monotonic evolution, leading to kinematic hardening. On the contrary, the cumulative plastic strain, on which the Aifantis theory relies, increases monotonically regardless the loading direction, leading to isotropic hardening.

Cyclic tension-shear response of a constrained plate: Equivalent SGP versus Tensorial SGP

A constrained layer subjected to combined tension-shear loading is considered in this subsection to compare the present theory with the tensor-based SGP theory of [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF]. The proposed loading allows for activating all components of the plastic strain tensor, making it possible to draw solid conclusions about the capability of the present model to reproduce the size-dependent effects classically produced by tensor-based SGP theories, while considering the gradient of a scalar field variable. The geometrical model involved in this subsection is similar to that presented in figure 2, with the top edge displacement in e 1 -direction replaced by two-component loading-unloading displacement in e 1 -direction and e 2 -direction: where γ and ε are respectively the mean shear and tensile strains. To ensure that the components of the plastic strain tensor have the same level, the imposed tensile strain ε is set twice the imposed shear stain γ (ε = 2 γ). The remaining boundary conditions are those described in Fig. 2 and Section 4.2.

u 1 (x 1 , w, t) = -u 1 (x 1 , -w, t) = w γ(t) and u 2 (x 1 , w, t) = -u 2 (x 1 , -w, t) = w ε(t)
In particular the equivalent plastic strain is set to zero at the top and bottom lines. In the case of the tensor-based SGP model, these latter conditions are replaced by prescribed vanishing components of the plastic strain tensor, excepted ε p 33 which is left free under plane strain conditions (ε 33 = 0). The material parameters are similar to those given in table 1, except for the first-order hardening parameter h (taken as zero h = 0) and c (taken as a variable parameter). The value of this parameter (c = H l 2 ) is adjusted via the internal length scale l, assuming fixed higher-order hardening parameter H = 10 GPa. The same mesh refinement, consisting of 20 × 20 finite elements, is adopted for both the equivalent and tensorial SGP simulations. However, different types of finite elements are used. In the equivalent SGP simulations, C 1continuous finite elements are applied. Application of the same type of finite elements for the tensorial SGP simulations would be extremely costly, as this will result in 64 degrees of freedom per finite element. Eight-noded fully quadratic elements with 2 × 2 Gauss integration are chosen to perform the latter simulations. This choice represents a good compromise between accuracy and computation cost.

Numerical implementation based on this type of elements has been performed in ABAQUS/Standard using UEL subroutine. For more details about the numerical procedure, the reader is referred to the paper of [START_REF] Martínez-Pañeda | The role of plastic strain gradients in the crack growth resistance of metals[END_REF], which is used as basis to perform the present implementation.

Figures 10 and 11 present the simulation results, in terms of distributions of cumulative plastic strain and plastic strain components within the studied plate, obtained using the equivalent and the tensorial SGP models, with imposed tensile strain ε = 2 γ = 0.02 and energetic length scale l/w = 0.4. It can be verified that all the plastic strain components are activated in the considered plate. The computation times to obtain such results are approximately 71 min for the proposed model and 150 min for the tensorial model. This difference between the computation times is explained by the difference between the total number of degrees of freedom involved in the simulations (4180 degrees of freedom for the present model simulations versus 6200 for the tensorial model simulations) and the resolution steps involved in each model implementation. The proposed model has allowed for approximately a twofold reduction of computation time. Qualitatively, the plastic strain distributions from both models are in good agreement. However, differences can be observe from a quantitative point of view. In particular the thicknesses of the boundary layers differ for both models. Figure 12 illustrates the simulation results in terms of overall shear stress as a function of imposed shear strain, obtained using the equivalent and the tensorial SGP models under a tension-shear loading to ε = 2 γ = 0.01 followed by unloading. The kinematic nature of the size-dependent hardening generated by the two models is illustrated in this figure. Using the same energetic length scale l, different levels of size-dependent kinematic hardening effects are obtained. The present model underestimates such effects. This is, in fact, expected as the back-stresses involved in both models are not the same. The back-stress implied in the present model is based on the Laplacian of the equivalent plastic strain (39), whereas that implied in the tensorial model is based on the Laplacian of plastic strain tensor (23). By applying the Laplacian on the equivalent plastic strain, certain second-order derivatives of this variable can overlap or offset each other, which can result in increasing or decreasing the back-stress. Figure 12 shows that the proposed model can produce the same level of size-dependent kinematic hardening effects as for the tensorial model by adjusting the energetic length scale l (figure 12). Interestingly, with this choice of intrinsic length, we find that the profiles of cumulative plastic strain are the same for both models, which indicates that the differences in the models' responses stem from the distribution among plastic components and not from the cumulative strain.

Conclusion

The present paper deals with the development of a new scalar-based strain gradient approach allowing for predicting size-dependent kinematic hardening effects. Unlike conventional scalar-based gradient theories, the present one is based on the gradient of the equivalent plastic strain which does not increase monotonically. The proposed theory was first analyzed theoretically using simple monotonic loading conditions. Then, numerical simulations under cyclic loading conditions were performed to investigate the nature of the resulting hardening. Comparison with the classical Aifantis (scalar-based) and Gudmundson (tensor-based) theories was given. These theories were implemented in the commercial finite element [START_REF] De Borst | Some Novel Developments in Finite Element Procedures for Gradient-Dependent Plasticity[END_REF] was applied to implement the present and the Aifantis theories, making the comparison between them more reliable. This formulation is widely used in the context of scalar-based gradient theories. Using the same formulation to implement the tensor-based theory of Gudmunson would be computationally expensive. An implementation based on eight-noded fully quadratic finite elements is proposed for the latter theory. This choice represents a good compromise between accuracy and computation cost. Main analytical and numerical results obtained in this work are recalled hereafter.

Using the gradient of the equivalent plastic strain instead of the usual gradient of the cumulative plastic strain, the developed theory overcomes the conceptual problem of scalar-based theories regarding the indeterminacy of the plastic flow direction under certain loading conditions [START_REF] Wulfinghoff | Conceptual difficulties in plasticity including the gradient of one scalar plastic field variable[END_REF].

Although based on the gradient of a scalar field variable, the proposed theory takes into account the plastic flow direction, making the solution stable with respect to the boundary conditions. Indeed, the back-stress resulting from the gradient of the equivalent plastic strain is calculated based on the direction of the plastic strain tensor.

Under monotonic loading conditions, it was found that the results of the proposed model mirror those of the classical Aifantis theory. The same results were obtained under monotonic shear and bending conditions. Under bending conditions beyond the limit load, both models predict non-smooth plastic strain distributions at the neutral axis. [START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF] have shown that this non-smooth distribution can be avoided by using the gradient of the plastic strain tensor. However, with the current understanding of small scale plasticity, it is early to confirm whether this kind of non-smooth distributions is physically acceptable. Other works showing non-smooth plastic strain distributions in the context of generalized continua can be found in the literature [START_REF] Stupkiewicz | A minimal gradient-enhancement of the classical continuum theory of crystal plasticity. part ii: Size effects[END_REF][START_REF] Rys | A micromorphic crystal plasticity model with the gradientenhanced incremental hardening law[END_REF]. The C 1continuous element formulation adopted in this work is not able to capture non-smooth distributions. To overcome this difficulty, an alternative approach is to use a micromorphic model for which the plastic strain is non-smooth but the micromorphic variable remains smooth [START_REF] Peerlings | On the role of moving elastic-plastic boundaries in strain gradient plasticity[END_REF][START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF]. In SGP, one could use implementation techniques allowing for discontinuities, such as discontinuous Galerkin methods [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF][START_REF] Eymard | Discontinuous Galerkin gradient discretisations for the approximation of second-order differential operators in divergence form[END_REF] or hybrid higher-order (HHO) methods [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF]. This point will be treated in a future work.

Application of cyclic shear loading conditions raised a strong difference between the proposed and Aifantis theories, regarding the nature of the resulting hardening. Using cumulative plastic strain gradient, the internal length scale leads to size-dependent isotropic hardening, which increases with increasing l. On the contrary, using equivalent plastic strain gradient, size-dependent kinematic hardening with classical Bauschinger effects is obtained. The kinematic hardening modulus increases with l. In the context of phenomenological strain gradient plasticity, the latter type of hardening is generally obtained using tensorbased gradient theories [START_REF] Poh | An implicit tensorial gradient plasticity model -Formulation and comparison with a scalar gradient model[END_REF], which are generally costly compared to scalar-based ones.

Comparison between the proposed theory and the tensor-based theory of [START_REF] Gudmundson | A unified treatment of strain gradient plasticity[END_REF] has shown that the former is more advantageous in terms of computation time. However, the models provide different levels of size-dependent kinematic hardening effects. This is expected as the models employ different expressions of back-stress. The level of the size-dependent kinematic hardening can be adjusted by appropriate choice of the energetic length scale.

The model proposed in this work can be viewed as a scalar-based counterpart of tensor-based theories, allowing for predicting size-dependent kinematic hardening effects at more affordable costs. For the price of two additional independent scalar field variables (equivalent and cumulative plastic strains), a model can now be proposed combining size-dependent kinematic and isotropic hardening effects using two distinct length scales. Although the present work is presented within small deformation framework, it can easily be extended to finite deformation, in particular to study strain localization phenomena in plasticity. Note also that illustrations of gradient extensions of the J 2 flow theory were given, but the approach equally applies to other equivalent stress measures. This includes plastic compressibility effects induced by the presence of the hydrostatic pressure in addition to invariants of the deviatoric part of the stress tensor.

The equivalent plastic strain measure would then incorporate the possibly non-vanishing trace of the plastic strain tensor.

Figure 1 :

 1 Figure1: One-dimensional model subjected to initial conditions ε p 0 and a small perturbation F , of which the direction is determined by the sign of .

Figure 2 :

 2 Figure 2: 2D strip subjected to simple shear loading and vanishing plastic stain on the top and bottom edges

Figure 3 :

 3 Figure 3: Thin foil subjected to simple bending loading and zero generalized traction on the top and bottom edges

Figure 4 :

 4 Figure 4: C1continuous element

Table 1 :

 1 Figure 5: Simple shear results obtained using the present model

Figure 6 :Figure 7 :

 67 Figure 6: Equivalent plastic strain profiles along e 2 -direction at different imposed curvatures

Figure 8 :

 8 Figure 8: Cyclic shear results obtained using the present and the Aifantis theories for different values of l

( a )

 a Cyclic shear response obtained using the present and the Aifantis theories for l/w = 0Equivalent plastic strain profiles at stages 1 to 3 of figure9a, obtained using the proposed model Cumulative plastic strain profiles at stages 1 to 3 of figure9a, obtained using the Aifantis model

Figure 9 :(Figure 10 :

 910 Figure 9: Profiles of the additional scalar field variables involved in the present and the Aifantis theories at different loading stages for l/w = 0.5

Figure 11 :

 11 Figure11: Distributions of plastic strain components within the studied plate, obtained using imposed tensile strain ε = 2 γ = 0.02 and l/w = 0.4

Figure 12 :

 12 Figure 12: Overall shear stress versus imposed shear strain for different values of energetic length scale l and three different models

In general, εeq is not considered as a physically relevant hardening variable.