A micromorphic crystal plasticity model with the gradient-enhanced incremental hardening law

A model of crystal plasticity is developed in which the effects of plastic flow non-uniformity are described through the full dislocation density tensor. The micromorphic approach is used in which the dislocation density tensor is represented by the curl of an independent constitutive variable called microdeformation. The microdeformation tensor is enforced by an energetic penalty term to be close to the actual plastic distortion tensor. The curl of microdeformation tensor enters the constitutive model in two independent but complementary ways. First, it is an argument of the free energy density function, which describes the kinematic-type hardening in cyclic non-uniform deformation. Second, its rate influences the rates of critical resolved shear stresses, which corresponds to additional isotropic hardening caused by incompatibility of the plastic flow rate. The latter effect, missing in the standard slip-system hardening rule, is described in a simple manner that does not require any extra parameter in comparison to the non-gradient theory. In the proposed model there are two independent internal length scales whose interplay is examined by means of 1D and 2D numerical examples of plastic shearing of a single crystal.


Introduction
The aim of this paper is to examine a crystal plasticity model that combines in a novel way two physically related but conceptually and mathematically distinct effects of plastic flow non-uniformity on the material hardening. The first one is due to the influence of the dislocation density tensor, being a measure of accumulated geometrical incompatibility of the plastic deformation, on the free energy density of the material. In the classical terminology, it leads to kinematic hardening since reverse deformation can reduce this effect to zero. The second effect is due to the influence of the current rate of dislocation density tensor, representing incompatibility of the current plastic flow rate, on the rate of critical resolved shear stresses. In particular, the rate of average density of the dislocations induced by the current slip-rate gradients can be assumed non-negative, which leads to isotropic hardening, also for cyclic deformations. The gradient effects are accompanied by the usual anisotropic hardening of a uniformly deformed crystal that encompasses self-hardening (diagonal) and latent-hardening (non-diagonal) terms of the slip-system hardening matrix.
To describe the first effect outlined above, it is proposed to include the full dislocation density tensor, defined by Nye (1953), Bilby et al. (1955), Kröner (1960), Le and Stumpf (1996), Steinmann (1996) and Cermelli and Gurtin (2001), among others, into the constitutive setting. Such strain gradient plasticity models have been developed for instance by Gurtin (2006) and more recently by Kaiser and Menzel (2019). The use of the full dislocation density tensor instead of individual densities of geometrically necessary dislocations as done for instance in (Gurtin, 2000;Bayley et al., 2007;Gurtin et al., 2007;Bargmann et al., 2014) is more efficient from the computational point of view and contains essential features of strain gradient plasticity as discussed in (Mesarovic et al., 2015;Wulfinghoff et al., 2015). Strain gradient plasticity based on the dislocation density tensor can be viewed as a limit case of the microcurl model proposed by Cordero et al. (2010) which results from the application of the micromorphic approach by Forest (2009) to gradient crystal plasticity. According to the micromorphic theory by Mindlin (1964); Eringen and Suhubi (1964), the material point is endowed with new degrees of freedom, namely a second order (generally) non-symmetric tensor called microdeformation, in addition to the usual displacement vector. In the microcurl model, this microdeformation and its curl are introduced as constitutive variables in the free energy density function. A penalization is introduced such that the microdeformation almost coincides with the plastic distortion tensor. As a result the curl of the microdeformation tensor will practically coincide with the usual dislocation density tensor. This approach can be regarded as an efficient method to implement strain gradient crystal plasticity in finite element codes. The microcurl model was applied to the continuum modeling of grain size effects in polycrystals by Cordero et al. (2012) and to Bauschinger effects in laminate microstructures (Forest, 2008;Wulfinghoff et al., 2015). The microcurl model was extended by Alipour et al. (2019) to include grain boundary yielding effects. A simplified version of the microcurl model was proposed in (Wulfinghoff and Böhlke, 2012;Wulfinghoff et al., 2013;Erdle and Böhlke, 2017). Instead of considering a microdeformation tensor, these authors introduce a scalar micromorphic variable related to the accumulated plastic slip, and its gradient into the free energy density function. The advantage of this model is the reduction of complexity from the computational mechanics point of view. The drawback is that this enhancement does not induce size-dependent back-stress effects in contrast to the original microcurl model. The full micromorphic deformation tensor will therefore be used in the present work.
To describe the second effect, i.e. of the current incompatibility of the plastic flow rate on the isotropic part of the hardening rate, the approach proposed recently by Petryk and Stupkiewicz (2016) is used. Incompatibility of the plastic distortion rate is measured by the effective slip-gradient rate,χ, defined as a norm of the rate of the dislocation density tensor. The associated evolving internal length scale, , has been derived from phenomenological relationships established in the dislocation-based theory of plasticity and expressed through standard parameters of a non-gradient hardening law, including the material constants in Taylor's formula (Taylor, 1934) and the current stress and hardening modulus θ. The product θ χ added to the conventional formula for anisotropic hardening rate for a single crystal constitutes its simple gradient-enhancement. It has been shown that this modification alone can provide realistic predictions of the experimentally observed indentation size effect . An evolving length scale was also introduced recently into other models of strain gradient plasticity (Dahlberg and Boåsen, 2019;Scherer et al., 2019).
The present paper extends the earlier works to a novel combination of the 'microcurl' model of micromorphic type (Cordero et al., 2010) and the 'minimal' gradient-enhancement of the incremental hardening law . It is distinct from the related work by Ryś and Petryk (2018) which was limited to the Gurtin-type models of gradient plasticity (Gurtin, 2000;Gurtin et al., 2007;Bargmann et al., 2014) combined with the gradient-enhanced hardening law. In particular, unlike in the previous works, the effective slip-gradient rateχ is expressed here by a norm of the curl of the micromorphic rate-variable. The numerical treatment here is also different since the micromorphic variable acts as an additional degree of freedom of a global type, so that the dual-mixed formulation applied in (Ryś and Petryk, 2018) is not needed here. One advantage of the micromorphic approach to strain gradient plasticity is its numerical efficiency due to the introduction of independent degrees of freedom connected at the constitutive level to standard mechanical variables, as demonstrated in contributions related to strain localization phenomena by Anand et al. (2012);Peerlings et al. (2012). In particular there is no need for front-tracking methods for the interface between plastically loaded zones and elastic domains of the structure, in contrast to some strain gradient implementation methods (Liebe et al., 2003). Numerical efficiency of the micromorphic-type modelling has recently been confirmed in 3D setting, although only for a scalar micromorphic variable, in the study of martensitic phase transformation by Rezaee-Hajidehi et al. (2019).
The paper is organized as follows. In Section 2 the microcurl model in crystal plasticity is presented in the small deformation framework. The field equation and boundary condition for the generalized couple stress tensor are derived in two alternative ways, with and without the use of the method of virtual power. In Section 3, the condition for plastic flow is derived using the compatibility of actual and virtual dissipation rates. The minimal gradient enhancement of the incremental hardening law with an associated natural length scale is presented along with its adjustment to the micromorphic approach. Sections 4 and 5 contain several numerical examples that visualise two different effects of plastic flow incompatibility on the material behaviour that come from the gradient-dependence of the condition for plastic flow and incremental hardening law. In section 4, 1D examples of plastic simple shear of a two-phase laminate and shearing of a constrained strip are analyzed, and numerical results are compared to analytical ones. In section 5 the analysis is extended to 2D finite element examples of cyclic simple shear of a square single crystal and an idealized polycrystal. Throughout this paper, the attention is limited to quasi-static isothermal deformation.

Kinematics
In the small deformation framework adopted in this paper, the spatial displacement gradient, denoted interchangeably as ∇u ≡ u ⊗ ∇, is split additively into elastic (H e ) and plastic (H p ) parts, The second-order tensors H, H e , H p are generally non-symmetric and can be decomposed into their symmetric and skew-symmetric parts, respectively. The rateḢ p of plastic distortion of a single crystal readṡ whereγ α is the slip rate on the α-th slip system defined by slip direction s α and slip-plane normal m α whose scalar product vanishes to ensure plastic incompressibility, s α · m α = 0. As usual, ⊗ stands for a tensor product. Geometrical incompatibility of a spatial field of H p (or −H e , equivalently) is measured by its curl, commonly called Nye's tensor (Nye, 1953), defined here using the following sign convention where e i denotes an orthonormal basis in Euclidean space R 3 and jkl the permutation symbol, with the summation convention for repeated indices. Following the micromorphic approach (Forest, 2009) that extends the original theory by (Mindlin, 1964;Eringen and Suhubi, 1964), as a counterpart to H p a micromorphic variable κ p is introduced and treated in the calculations as a global variable corresponding to local H p .

Free energy density
In the microcurl model of crystal plasticity (Cordero et al., 2010), extended to include local internal variables in analogy to (Steinmann, 1996), (Menzel and Steinmann, 2000) and (Aslan et al., 2011), the Helmholtz free energy density function per unit volume, ψ, at a given temperature is assumed to have the four arguments: the elastic strain ε e , local internal variables ξ = (ξ α ), the relative plastic strain e p := H p − κ p , and Γ p as the curl of micromorphic variable κ p . In the computational version of the model used in this paper, ψ is assumed in the following additive form where, using a central dot to denote full contraction of tensors, In the previous functions, the fourth order tensor C is the standard tensor of the elastic moduli determined by two independent moduli µ, ν in the isotropic case. For the sake of simplicity, the additional fourth order tensors H κ and A are characterized by the generalized moduli H κ and A, respectively, assuming H κ = H κ 1 and A = A 1, where 1 denotes the fourth-order identity tensor operating on non-symmetric second order tensors. The choice of a quadratic potential ψ curl is justified by micromechanical analyses like the consideration of stacked pile-ups in (Baskaran et al., 2010) or direct comparison with discrete dislocation dynamics in (Chang et al., 2016). Other higher order energy potentials have been proposed involving the norm of the dislocation density tensor or the logarithm of the norm (Berdichevsky, 2006;Le and Sembiring, 2008;Forest and Guéninchault, 2013). However, these potentials are associated with singularities which significantly complicate the numerical analysis (Wulfinghoff et al., 2015).
There are two possible interpretations of the H κ material parameter in ψ micro . It can first be seen as a numerical regularization parameter in order to implement strain gradient plasticity in a simple way. It should then take sufficiently large values. It can also be regarded as a true material parameter that requires calibration from appropriate experimental or numerical data. For instance, Cordero et al. (2012) calibrated values for H κ so as to reproduce the Hall-Petch relationship from finite element polycrystal simulations, at least in a certain range of grain sizes. It involves intermediate values of H κ that will be explored in the examples of Section 4.
The term ψ p is needed for thermodynamic consistency of the model that includes statistical storage of dislocations, as it represents mechanically unrecoverable changes in the residual free energy which do not contribute to isothermal dissipation. Its rate is assumed in the formψ which can be given a physical interpretation but may be left unspecified in the calculations, see (Ryś and Petryk, 2018) and section 3.1 for more details. Quantities related to partial derivatives of function ψ are as follows where σ is the (symmetric) Cauchy stress tensor, p α -internal forces, s -the (non-symmetric) microstress tensor, M -the generalized couple stress tensor, see section 2.4. The stress σ and negative microstress −s resolved on the α-th slip system are denoted by where τ α is the standard resolved shear stress and x α plays below the role of a kinematic hardening component (back-stress). The rate of free energy density ψ is calculated straightforwardly as followṡ The quantity that appears in the above expression forψ may be identified with the thermodynamic driving force conjugate toγ α that results from the assumed form of the free energy density function ψ.

Balance equations
We consider a continuous deformable body that in the geometrically linear setting occupies a domain V of boundary ∂V in an Euclidean space R 3 . The assumed basic kinematic fields are (Cordero et al., 2010) {u,κ p } .
The fieldsu andκ p are assumed to be continuous and sufficiently smooth so that their gradients along with the divergence theorem have a mathematical sense. Note that the field of slip ratesγ α is not treated as a global kinematic variable, due to the fact that the proposed theory is based on the full dislocation density tensor (curlκ p ) instead of the individual GND densities.
When neglecting body and inertia forces, the power density of internal forces in V is assumed in the following form, consistent with Eq. (10), and the power density of external contact forces on ∂V as where t is the traction vector and m -the double traction tensor. The respective volume integrals are assumed to be equal, by the virtual power equality (Germain, 1973) that is assumed to hold for arbitrary fields ofu andκ p and for every subset Π of V . The integral P int is transformed with the help of the divergence theorem as follows where n is an external normal to ∂V .
Since equality (16) is assumed to hold for arbitrary fields ofu andκ p then, following (Cordero et al., 2010), we obtain the field equations in V along with the boundary conditions over ∂V where (M ) ikl = M ij jkl .

Dissipation inequality
By the first law of thermodynamics (for isothermal and quasi-static processes) we havė where F = V ψ dV is the Helmholtz free energy andḊ is the dissipation rate which is non-negative by the second law of thermodynamics. Based on this and using Eq. (16), the following frequently used local form of dissipation inequality is obtained On taking the time derivative of Eq. (5) in the form (10) but without making use of Eqs. (8), it follows that Eqs. (8) 1 , (8) 3 and (8) 4 are recovered if the inequality is required to hold identically, by the usual argument assuming thatε e ,ė p andΓ p are unconstrained. On substituting Eqs. (7) and (10) into inequality (22), the intrinsic dissipation rate inequality is expressed in the reduced form

An alternative approach
Instead of assuming the form of both internal and external power densities (13) and (14) in advance, one can calculate the dissipation rate from Eq. (20) without using Eq. (16). In analogy to (Ryś and Petryk, 2018), for every subdomain Π ⊆ V the external power is taken in the form obtained by replacing the last term in Eq. (14) with a yet unspecified micromorphic extension p ext m of the standard mechanical power of contact tractions. It is the external power expression which enables classification of continuum theories in a manner separated from constitutive assumptions (Del Piero, 2014a,b). On substituting Eqs. (10) and (24) into Eq. (20), the total dissipation rate in the state of mechanical equlibrium readṡ By using the divergence theorem to split the total dissipation rate into the bulk and surface terms, it follows thaṫ The last integrand can be of either sign for any given M ij jklκ p ik = 0 since Π and hence the direction of an external normal n to boundary ∂Π of Π are arbitrary. This statement for a point on a short boundary segment ∂Π * can be extended, by employing a known argument, to the sum of the last two integrals for a sufficiently thin disk Π adjacent to ∂Π * wheṅ κ p = 0 on ∂Π \ ∂Π * , so that the contribution of the volume integral becomes negligible. To ensure thatḊ ≥ 0 for all Π ⊆ V and for every field ofκ p , the term P ext m must cancel out the above surface dissipation term over ∂Π . This is so if which is in agreement with Eq. (19) 2 but has been obtained on another route than in (Cordero et al., 2010). On substituting Eq. (27) into Eq. (26), we arrive at the following expression for the dissipation rateḊ Finally, equality s + curlM = 0 is obtained if an arbitrary field ofκ p is assumed not to affect the dissipation rate.
Remark 1. The right-hand equations in Eqs. (18) and (19) have been derived above without assuming the specific form (13) of the internal power density, and consequently without the need to appeal to the principle of virtual power (16) once the free energy density rate has been expressed by Eq. (10). Eqs. (18) 2 and (19) 2 have been obtained instead by using thermodynamic arguments: equation m = M · n by requiring the surface dissipation term to be always non-negative (and thus to be cancelled out by the external power term), and s + curlM = 0 by assuming that the micromorphic variable rateκ p is not associated with the dissipation rate. Hence, these basic equations of the microcurl model (Cordero et al., 2010) have been recovered on another route.

The plastic flow and hardening laws
3.1. The condition for plastic flow Following Ryś and Petryk (2018), the condition for plastic flow can be derived using the compatibility of the actual and virtual dissipation rates. They are distinct in general since the value of the actual one follows from the entropy imbalance while the virtual one is defined by a constitutive assumption, cf. (Petryk, 2005). Here, Eq. (23) is treated as an expression for the actual dissipation rate density. The virtual dissipation rate is defined pointwise by introducing a dissipation function D ofγ α , assumed to be independent ofκ p , To encompass both the rate-independent and rate-dependent dissipation, it is assumed that where τ c α is the critical resolved shear stress on the α-th slip system, strain-rate sensitivity coefficient m > 0 corresponds to rate-dependent dissipation and m = 0 to rate-independent dissipation, and p α appears in expression (7) forψ p . Ifψ p describes free energy increase due to statistical multiplication of dislocations then it is natural to identifyξ α with (ρ α ) S discussed in the next section. Accordingly, of compatibility of the actual and virtual dissipation rates, separately for each α and for every Π ⊆ V . In view of Eqs. (23) and (29), this is equivalent to imposing the pointwise condition which on substituting Eq. (11) yields the condition for plastic flow On using Eqs. (30) and (31) 4 , it can be rewritten as follows In the limit as m → 0, the last term becomes negligible, and the resulting condition for plastic flow takes the familiar form Recall that x α is the micromorphic back-stress defined by Eq. (9) 2 . Eq. (36) is complemented with the inequality constraint The evolution equation for τ c α (i.e. the hardening law) is addressed in the next section.
Remark 2. The rate-independent condition (36) for plastic flow has not been assumed a priori by analogy to the Schmid law but has been derived using the thermodynamic condition (32) of compatibility of the actual and virtual dissipation rates under assumption (29). Remarkably, condition (36) is independent of p α , in analogy to the non-gradient case (Petryk and Kursa, 2015), hence p α may be left unspecified in calculations. On the other hand, p α is needed in a consistent thermodynamic framework that includes statistical storage of dislocations, especially if rate-sensitivity is also encompassed.
3.2. The hardening law with slip-rate gradient effect This section follows closely Petryk and Stupkiewicz (2016), with the distinction that the role of the effective slip-rate gradient will be played by a norm of the micromorphic variable curlκ p . According to the (generalized) Taylor formula (Taylor, 1934), which is one of the basic phenomenological laws in the materials science literature on plasticity of metals, an isotropic flow stress τ c is a function of the total dislocation density ρ, where coefficient a, elastic shear modulus µ and Burgers vector modulus b are given material constants. To determine the rate of τ c , the rateρ (and not ρ itself) is decomposed into the sum of the average density rate (ρ) S of statistically generated dislocations, and the average density rate (ρ) G of dislocations induced by the current slip-rate gradients, It means that in distinction to the approach by Ashby (1970), Fleck et al. (1994), Nix and Gao (1998) and others, not the existing dislocations themselves but rather their current sources are split according to their statistical or geometrical character. Hence, a time integral of (ρ) G need not define the current density of existing GNDs (geometrically necessary dislocations). By adopting the known formula for the rate of statistical accumulation of dislocations (Kocks and Mecking, 2003) in its simplest form, it is postulated that where λ denotes the dislocation mean free path, specialized as λ α for the α-th slip system, whose inverse is defined as the incremental mean length of dislocation stored per area swept.
Since the analysis in this paper is limited to small plastic strain with respect to an annealed state, annihilation of dislocations and their transport are not included in formulae (40). However, the formula (44) below for the internal length scale remains valid also if Eq. (40) 1 is extended to include another term that corresponds to dislocation annihilation . The geometrically induced dislocation density rate is postulated in another simple form, whereχ is the effective slip-rate gradient, assumed here to be represented by the Euclidean norm of a curl of the micromorphic rate-variableκ p . In case of cyclic loadings, a time integral ofχ scaled by b cannot be interpreted as the current GND density, rather, as the nondecreasing contribution to a total dislocation density coming from the history of incompatible plastic distortion rate. Combination of Eqs. (39)÷(41) with the time derivative of formula (38) defines the contribution ofχ to the rate of the isotropic flow stress τ c , that is missing in the conventional incremental hardening law. Consistency with the standard formulaτ c = θγ, where θ is the scalar hardening modulus in the non-gradient case, requires thatτ If function τ ρ is specified by the Taylor formula (38) 2 then the evolving internal length scale is expressed by It is pointed out that only the standard quantities of a non-gradient hardening law appear on the right-hand side of Eq. (44), so that no additional assumption is needed to determine . Moreover, has a physical interpretation through its direct link to the dislocation mean free path λ which is a well-known length-scale parameter in the physically-based dislocation theory of plasticity.
The conventional incremental hardening law can be expressed as followṡ Assuming that only the isotropic partτ c is influenced by nonzeroχ and substituting Eq. (43), the gradient-enhanced anisotropic hardening law is obtaineḋ The last term that includes a 'natural' length scale defined by Eq. (44) will be referred to as 'P-S term' . It provides a 'minimal' gradient-enhancement of the incremental hardening law as a consequence of Eqs. (38)-(41). Note that the coefficient θ atχ varies in time inversely proportional to τ c .
Hardening moduli h αβ are frequently taken in the form with Here, q is a latent hardening parameter, and θ α is a self-hardening parameter for the α-th slip system, defined by (Anand and Kothari, 1996) where θ 0 , τ max and p are constant parameters, the same for all slip systems. For isotropic counterparts θ and τ c the index α is omitted, so that

Shearing of a strip
In the following section we consider two typical examples of shearing of a strip. The first example concerns plastic simple shear of a single crystal with the slip direction and slip plane normal collinear with global x-y ⇔ 1-2 coordinate system as in (Cordero et al., 2010) and (Aslan et al., 2011). The second example addresses shearing of a constrained strip with two symmetric slip systems as in Stupkiewicz and Petryk (2016). The sketches of the problems under consideration are shown in Fig. 1. In the following 1D examples, 500 elements of a regular mesh and quadratic shape functions for both fields are used.

Single slip
In the first example the strip consists of soft phase (s) in the middle and hard phase (h) on both sides imitating hard elastic inclusions which form obstacles to dislocation movement (Fig. 1a). Such a two-phase laminate is subjected to plastic simple shear γ(x) by applying a mean shear strainγ in the x direction. The periodic boundary condition along x-axis is imposed on the displacement u 2 in y direction and microdeformation κ p . We consider displacement and microdeformation fields of the form: For so-defined problem there is only one resolved shear stress τ = σ 12 = µ(γ − γ + u 2,x ), one component of the generalized couple stress tensor M 13 = −Aκ p 12,x and two components of the microstress tensor s 12 = −H κ (γ − κ p 12 ) and s 21 = H κ κ p 21 . However, from the balance equation s = −curlM = −A curl (curl κ p ), we obtain and thus s 21 = 0. Hence, the plastic flow criterion is obtained in the following form So far the analysis coincides with that in (Cordero et al., 2010) and (Aslan et al., 2011). Now, we extend it by including the P-S term as in Eq. (46), which in the single slip case readsτ From the above equations, it is clear that size effects in the current model stem from two internal length scales mic and , where the former is related to curl(curlκ p ) and the latter to curlκ p . An analytical solution of the above problem, in the case when the P-S term is absent ( = 0) and for linear hardening (p = 0 in Eq. 49), was derived in Aslan et al. (2011). The main equations of the analytical solution are provided in Appendix A. If the P-S term is included then the rate equation resulting from Eq. (53) involves a variable coefficient θ that varies both in time and space. The authors have not found a way to solve that equation analytically.
In the analysis it is assumed that soft (s) and hard (h) phases have different parameters H κ and A specified by the corresponding superscripts. The 1D results, shown below, are calculated for l = 1µm and the soft phase fraction f s = 0.7. Importantly, the value of penalty parameter H κ was taken such that profiles of microdeformation variable κ p 12 and plastic slip γ nearly coincide in the soft phase.
In Fig. 2, analytical and numerical results are presented for the plastic microdeformation variable κ p 12 for different values of A s = 0.05, 1, 50 GPa·µm 2 (i.e. millinewtons) and the same value of H κ = 500 GPa (thus in the soft phase mic = 0.01, 0.045, 0.32 µm), while the other parameters were fixed and are given in Table 1 (with annotation a ). The results computed without the P-S term coincide with the analytical results of Aslan et al. (2011). For smaller values of parameter A s the parabolic profile of microdeformation is observed while for larger values the profile is almost flat in the soft phase. This is related to the dislocation pile-up at the interfaces. In the case of small difference in the material paramater A in two phases (A s ≈ A h ) the GNDs pile up at the interfaces but they can also be smoothly distributed within the soft phase. On the other hand, when increasing the difference between the parameters (but keeping A s > A h ), GNDs tends to concentrate in the close vicinity of the interfaces. Including P-S term in the model increases the maximum value of microdeformation variable in the soft phase. This is more noticeable for smaller values of A s (Fig. 2b).  The effect of the P-S term (Eq. (46)) depends on the evolving length scale, , which is not adjustable as it is uniquely expressed through standard parameters of a non-gradient hardening law. However, the relative strength of this effect depends on the strip size l and is also influenced by other parameters, cf.  and (Ryś and Petryk, 2018). For this reason, for some range of parameters values, the P-S term may play a significant role in the model or not. In the above results, the hardening parameter θ 0 is quite large and despite the difference in microdeformation profiles that are noticable, the differences in stress-strain plots, which are not shown here, are negligible. To show how the P-S term affects the stress-stain response we compute the same example but with parameters taken from (Ryś and Petryk, 2018) (Table 1 with annotation b ).
In Fig. 3a the macroscopic stress Σ 12|γ=0.05 at 0.05 average plastic strain as a function of the length l of the two-phase microstructure is presented in a log-log diagram. The following results are for a constant value of A s = 0.05 GPa·µm 2 (A h = A s /100) and different values of coupling modulus H κ = 1, 10 and 100 GPa. Analytical results are compared with numerical ones in the case where P-S term was not included, i.e. χ = 0, and perfect coincidence was obtained. In the case where P-S term was included in the model, a significant increase in macroscopic stress is observed for the range of l values between 0.1 and 100 µm, with the maximal difference for l ≈ 2 µm (Fig. 4a). It should be noted, however, that the range and significance of the effect of P-S term may change with the value of A s (Fig. 4b). On the other hand, the difference (γ − κ 12 ) is sensitive to the value of H κ at small l.
The influence of P-S term on the maximum value of microdeformation κ p 12 , which is in the middle of the structure, has been shown in Fig. 3b. The maximum values are not only higher when P-S term is included but also shifted to the range of higher values of l (Fig. 4a).

Double slip
In the second 1D example we consider monotonic shear of a constrained strip of thickness H as shown in Fig. 1b. For this problem, a semi-analytical solution has been derived by Stupkiewicz and Petryk (2016) by using the 'minimal gradient-enhancement' of the classical model, i.e. when in the present formulation the terms ψ micro and ψ curl in the Helmholtz free energy are not included. Then there is no length-scale effect related to the second gradient of plastic slip, but there is another length-scale effect due to the first gradient of plastic slip rate, accounted for in the P-S term. The 'minimal' gradient-enhanced model was originally formulated within the classical continuum theory, in which no additional actions except the standard ones are included in the external power. In  it was also indicated that it required some regularization in order to obtain numerical solutions effectively, and non-local slip-rates were introduced as additional global unknowns that average and smoothen the local slip-rates. An averaging equation inspired by the so-called implicitgradient model was used, which in 1D case is basically analogical to Eq. (52) 2 but with an element scale parameter l h rather than mic = A/H κ . Such an algorithmic approach was also used more recently in (Lewandowski and Stupkiewicz, 2018), where in certain circumstances oscillations in numerical solutions for the wedge indentation problem were reported. Similarly, spurious oscillations can occur when the so-called dual-mixed method is used and the regularization is insufficient in vicinity of a kink in the analytical solution (Ryś and Petryk, 2018). Since the P-S term may have a crucial impact on predicting correctly the experimentally observed indentation size effect at the micron scale , it is of interest to check whether a micromorphic model can be more suitable for regularization purposes. It can be expected that the smaller the value of mic (while keeping H κ sufficiently high) the closer the solution to the semi-analytical one obtained in Petryk and Stupkiewicz (2016).
There is one non-zero component of the generalized couple stress tensor M 23 = Aκ p 21,y and two components of the microstress tensor s 12 = −H κ (H p 12 − κ p 12 ) and s 21 = −H κ (H p 21 − κ p 21 ). The balance equation s = −curlM results in Since H p 12 = κ p 12 , the component s 12 does not appear in the flow condition. Note also that κ p 12 does not contribute to the dislocation density tensor. The plastic flow criterion takes the form |σ 21 cos(2φ) − s 21 sin 2 (φ)| = |σ 21 cos(2φ) − Aκ p 21,yy sin 2 (φ)| = τ c , In the calculations of this example, the penalty parameter H κ = 1000 GPa is adopted which has been found sufficiently high for the fields of H p 21 and κ p 21 to be practically coincident. The value of parameter A = 0.1 GPa·µm 2 has been taken small so that the influence of the microstress s 21 be small, but on the other hand, A cannot be too small to provide the required numerical regularization. Other parameters are as in Table 1, annotation b . Comparison of the analytical and numerical results is presented in Fig. 5, where solid lines correspond to analytical results for A = 0 while dashed lines represent numerical results. In Fig. 5a the overall response it terms of the normalized shear stress,σ xy = σ xy /σ xy,0 (with σ xy,0 = | cos 2φ|τ 0 ) versus the overall shear strain is plotted, and in Fig. 5b the shear strain profiles, γ xy (y) at the overall shear strain γ xy = 0.05, are shown for different values of strip thickness H. The size effect, driven almost solely by the length scale , is clearly visible in both figures. The kink in γ xy profiles, resulting from the analytical solution for γ, is more noticable for the smaller strip thickness. The kink in γ 1 solution results in the jump in γ 1,y . In Figs. 5b and 6 it can be seen that the kink as well as the jump were properly smoothed by the finite element solution. Importantly, a better agreement between numerical and experimental results has been obtained than in the previous work where the dual-mixed method was used (Ryś and Petryk, 2018). In particular, no oscillations were observed on γ 1,y and γ 1,yy profiles, neither for local variable γ 1 nor for its global counterpart κ p 21 /(2 sin 2 φ). It can be concluded that the micromorphic approach provides a better reqularization technique in this particular case. However, it is worth mentioning that if a linear shape function for κ p 21 was used instead of the quadratic one then some oscillations occurred in vicinity of the kink. (b) Figure 6: Profiles of the (a) first and (b) second gradient of (local) plastic slip γ 1 , and global counterparts κ p 21,y /(2 sin 2 φ) and κ p 21,yy /(2 sin 2 φ).

Numerical 2D examples
In this section, we consider shearing of a square under plane strain conditions. First, a single crystal is considered and then a square composed of four differently oriented grains. In both cases, each grain deforms plastically on two slip systems rotated relative to each other by an angle 2φ = 2π/3. In the first case the systems are oriented symmetrically to the direction of shear (Fig. 7a), and in the second case, the slip system pairs are rotated relative to each other as shown in Fig. 7b. In the single crystal, displacement in the vertical direction is forbidden at all nodes on the external edges of the square, and zero displacement in the horizontal direction is also prescribed on the bottom edge, while uniform horizontal displacement is prescribed on the upper edge. Plastic slips are indirectly constrained by κ p = 0 at the external edges of the single crystal. In the case of the square with four grains, periodic boundary conditions for the displacement and micromorphic variable κ p are imposed over the external boundary, the shear is enforced at the corner points, while on the grain interfaces continuity of fields u and κ p 12 is assumed (cf. Sec. 5.3). The mesh of a crystal or grain is built with 100×100 regular square finite elements with 9 Gauss points in each element, and quadratic shape functions are used for both fields. The computations have been carried out using the AceFEM package. The material parameters used here, listed in the Table 2, correspond to properties of Cu single crystals and are taken partially from (Anand and Kothari, 1996;Sauzay and Kubin, 2011;Stupkiewicz and Petryk, 2016) following (Ryś and Petryk, 2018). For material parameters as in Table 2 and p = 1, the internal length scale decreases from an initial value of 8.8 µm to a minimum value of 3.4 µm as τ c increases from τ 0 to 1 2 τ max , and then starts to increase. The rate-dependent version of the model is used, specified by Eq. (30) with p α = 0. The adopted value of parameter A corresponds to an energetic length scale L en = 0.2 µm, see below.

Comparison of energy functions of different GND measures
In this section, we compare numerical results for two free energy density functions in which different GND measures are used. The present model with the energy function as in Eq. (6) is compared with the gradient plasticity model with the following form of the free energy part dependent on GND densities ρ G α (Bargmann et al., 2014), Coupling between slip systems is represented by the sum of the coplanarity moduli κ αβ , Eq. (48), and interaction moduli ι αβ defined via (Gurtin and Reddy, 2014) ι αβ = |s α · s β ||m α × m β | .
In the present model, in the limit case as κ p → H p , the defect free energy can be written as The curl of the tensor H p , with the sign convention as in Eq. (4), can be rewritten in the form In the case of 2D numerical simulations it is sufficient to consider only edge dislocations, thus the second term in the above sum vanishes. Taking A = A1 and substituting Eq. (60) 2 , we obtain In the case of two symmetric slip systems which are specified by φ = π/3, full contraction of the last two terms in the above equation gives 1 for α = β and −0.5 for α = β, thus in this particular case, the free energy becomes where δ αβ is the Kronecker delta and ς αβ = 0 for α = β and ς αβ = −0.5 for α = β. Taking that A = h 0 L 2 en and noting that κ αβ = δ αβ in the present case, the two energies (60) and (62) differ only by the non-diagonal coefficients ς αβ and ι αβ of GND interaction between different systems, where ι αβ = 0.43 for α = β and 0 otherwise. In particular, taking h 0 = µ 8(1−ν) after (Evers et al., 2004) and L en = 0.2 µm, the value A = 0.288 GPa·µm 2 given in Table 2 is obtained. A more general analysis regarding differences between various forms of the defect energy was presented in Mesarovic et al. (2015).
In the example below, the square of side length H = 10 µm is subjected to cyclic shear by prescribed displacement u = ±H/10 on the upper edge. The rate dependent model was used with rate sensitivity parameter m = 0.02, hardening parameters p = 2, q = 1.4 and others parameters listed in Table 2. In Fig. 8, the overall shear stress vs overall shear strain for the two types of defect free energies are plotted. Differences between the results are very small due to the relatively small value of L en = 0.2 µm; it has been checked that for L en = 0.7 µm the differences are no longer negligible. Distributions of the plastic slips, |γ α |, and effective plastic slip, γ = t α |γ α |dt, are presented in Fig. 9. The fields of components of curlκ p as measures of the geometrically necessary dislocation density are presented in Fig. 10. Plastic slips and GND distributions are qualitatively and quantitatively similar to those given in (Ryś and Petryk, 2018). The distribution of the GND measures, however, is smoother and more regular in comparison to the previous results, which confirms the advantages of the micrcromorphic approach in regularizing the numerical problem solved.

Effect of the P-S term and the material sample size
In the present work, the GND density tensor affects the material response in two ways, through Γ p in the free energy function and through the P-S term in the hardening law. If p α = 0 then the P-S term is dissipative in nature, leading to the increase of critical thermodynamic forces due to newly created GNDs whose density rate is proportional to the rate of plastic incompatibility tensor. For illustration, two examples of shearing of a square of side length H = 10 and 1 µm are presented for a cycle up to the final overall shear strain of a small value γ xy = 0.01. Here, linear hardening is assumed with p = 0 and q = 1.4; the remaining parameters are listed in Table 2. In Fig. 11, overall shear stress vs overall shear strain is plotted for the two side lengths of the square. In both cases isotropic hardening is influenced by the P-S term, more significantly for smaller H as expected. On the other hand, the non-local back-stress effect due to the adopted free energy form is also clearly visible. Final distributions of plastic slips |γ α | and of the effective slip γ = t α |γ α |dt are presented in Figs. 12-15. The maximum values of plastic slips and effective slip are higher for the model with the P-S term included, but in the case of H = 1 µm the difference is much smaller.
Distributions of the components of the dislocation density tensor and its norm are presented in Fig. 16-19. Maximum values of GNDs are lower when the P-S term is included and the distributions are smoothed, i.e. the concentration of GNDs near the boundary is less significant.      Figure 16: Distribution of the components and norm of curlκ p as the GND density tensor, in µm −1 , for H = 10 µm, γ xy = 0.01, P-S term included. Figure 17: Distribution of the components and norm of curlκ p as the GND density tensor, in µm −1 , for H = 10 µm, γ xy = 0.01, P-S term disregarded. Figure 18: Distribution of the components and norm of curlκ p as the GND density tensor, in µm −1 , for H = 1 µm, γ xy = 0.01, P-S term included. Figure 19: Distribution of the components and norm of curlκ p as the GND density tensor, in µm −1 , for H = 1 µm, γ xy = 0.01, P-S term disregarded.

Idealized polycrystal
The last example is devoted to the analysis of a square composed of four grains with the same set of two slip systems but rotated in each grain as shown in Fig. 7b. Periodic boundary conditions for the displacement and micromorphic variable are applied. Importantly, we do not use any additional conditions on grain interfaces. The interfacial conditions arise from the balance equations for continuous and piecewise smooth fields u, κ p and related simple and double tractions t and m, i.e. any jumps of these fields across grain interfaces are excluded, cf. (Cordero et al., 2012). This might be a shortcoming since physically motivated interfacial conditions have been discussed in the literature and introduced in the micromorphic framework in the form of interface yield conditions by (Wulfinghoff et al., 2013;Alipour et al., 2019). However, in this work our aim is to analyze the influence of the P-S term on plastic slip distribution and the GND pile up within the volume, and not on the dislocation movement across a grain interface.
The parameters are as in the preceding section. The overall shear stress vs overall shear strain for the square H = 20 µm is shown in Fig. 20. Similarly as previously, isotropic hardening that is significantly influenced by the P-S term and non-local back-stress effects are clearly visible. Differences in distributions of the plastic slips |γ α | and the effective slip γ = t α |γ α |dt in the two cases considered are presented in Figs. 20 and 21. Again, the maximum values of plastic slips and the effective slip are higher for the model with included P-S term. The use of P-S term also results in more pronounced accumulation of the plastic slips in the middle of grains. Distributions of the components of the dislocation density tensor and its norm are presented in Figs. 23 and 24. Interestingly, contrary to the single crystal examples, the maximum values of GNDs are higher and their concentrations are more localized when the P-S term is included.

Conclusion
A new version of the micromorphic model of crystal plasticity has been developed by combining the microcurl model (Cordero et al., 2010) with the minimal gradient-enhancement of the hardening law . The combination was not entirely straightforward because both component models required appropriate adjustments to achieve intrinsic consistency of the compound model. Among the modifications introduced, the condition for plastic flow has been derived, following Ryś and Petryk (2018), using the assumed compatibility of actual and virtual dissipation rates, and the effective slip-rate gradient in the enhanced hardening law has been identified here with the curl of the micromorphic ratevariable.
In result, the curl of the microdeformation tensor as a basic constitutive variable has entered the computational model in two complementary ways. First, it is an argument of the free energy density function, and its curl projected on a slip system defines the back-stress in the condition for plastic flow. Second, its rate enters the expression for the rate of critical resolved shear stresses, accompanied by a natural length scale whose value is evolving in a manner uniquely defined by standard parameters of a non-gradient hardening law. In general, the material behaviour is affected by plastic flow non-uniformity in both ways, although not necessarily to a similar extent.
It has been shown by the analysed examples in which circumstances either one or another effect of the curl of the microdeformation tensor can be predominant. Since the internal length scale in the hardening law is not adjustable, the freely adopted values of the energetic length scale mic and characteristic dimension l of the material sample can decide which effect is more substantial, although other material parameters also influence the results. More specific conclusions have been formulated in the preceding sections.
If all other parameters are fixed then by taking a sufficiently small value of the energy parameter A, which is proportional to the square of the length scale mic , the micromorphic approach can serve as a computational regularization tool. Its effectiveness for regularization purposes has been confirmed in the present paper, for instance, by eliminating spurious numerical oscillations which were present in certain solutions obtained using other regularization methods. Therefore, the micromorphic regularization with A small enough and H κ high enough is well suited for exploring predictive capabilities of the numerical modelling based on the minimal gradient-enhancement of the hardening law in estimating size effects in crystal plasticity at the micron scale. Dependence of hardness on the normalized penetration depth in the 3D spherical indentation problem provides a good example of such capabilities . The progress in experimental measurements of GND density fields (Dahlberg et al., 2014;Sarac et al., 2016) offers other possibilities of validating the model, cf. (Lewandowski and Stupkiewicz, 2018).
The proposed model has been presented and applied within the small strain framework. Extension of the approach to finite deformations is possible along the lines developed for strain gradient plasticity in (Gurtin, 2006;Kaiser and Menzel, 2019) and micromorphic crystal plasticity in (Aslan et al., 2011;Forest, 2016;Ling et al., 2018;Alipour et al., 2019). where constants τ c and H are denoted in the present paper by τ 0 and θ 0 , respectively.