## Abstract

We propose a visco-plastic strain gradient plasticity theory for single crystals. The gradient enhancement is based on an equivalent plastic strain measure. Two physically equivalent variational settings for the problem are discussed: a direct formulation and an alternative version with an additional micromorphic-like field variable, which is coupled to the equivalent plastic strain by a Lagrange multiplier. The alternative formulation implies a significant reduction of nodal degrees of freedom. The local algorithm and element stiffness matrices of the finite-element discretization are discussed. Numerical examples illustrate the advantages of the alternative formulation in three-dimensional simulations of oligo-crystals. By means of the suggested formulation, complex boundary value problems of the proposed plastic strain gradient theory can be solved numerically very efficiently.

## 1. Introduction

The continuum mechanics of single crystal models (Hill 1966; Lee 1969; Rice 1971; Ortiz & Stainier 1999) are based on a well-established physical basis. The discovery of dislocations as the fundamental carriers of plasticity justifies the mesoscopic kinematical assumptions concerning the plastic slip rates (Orowan 1934) and the slip system geometry associated with these models. Experiments (Schmid & Boas 1935) provide evidence of the correlation between the projected shear stresses and the slip system activity. Also non-Schmid effects can be taken into account within this setting (Kocks 1987; Yalcinkaya *et al.* 2008).

In contrast, single crystal hardening models and the prediction of the closely related dislocation micro-structure are still subjects of current research. Besides purely phenomenological hardening laws (Taylor 1938; Koiter 1953; Hill 1966) local dislocation density evolution or balance equations (Gillis & Gilman 1965; Mecking & Kocks 1981; Franciosi & Zaoui 1982; Estrin 1996) allow to estimate the evolution of the total dislocation line length per unit volume and to model the hardening behaviour of crystals more physically. These models are successfully applied at scales and specimen dimensions, where no size effect is observed. However, the inability to predict the size-dependent behaviour at smaller scales observed in experiments, e.g. the Hall–Petch effect (Hall 1951; Petch 1953) or the size-dependent strength of micro-specimens (Fleck *et al.* 1994; Stölken & Evans 1998; Xiang & Vlassak 2006; Dimiduk *et al.* 2007; Gruber *et al.* 2008) motivated many authors to enrich the classical local theory, for example, based on Nye's dislocation density tensor (Nye 1953; Bilby *et al.* 1955; Kröner 1958) or related measures associated with the gradients of the plastic slips. Nye's tensor can be computed from the spatial gradients of the plastic variables and introduces an internal length into the theory. It can be interpreted as the total Burgers vector per unit area, i.e. it has a clear physical meaning and must be distinguished from the above-mentioned dislocation density (the total line length per unit volume). Nye's tensor or familiar gradient-type quantities (Fleck & Hutchinson 1993; Liebe & Steinmann 2001; Liebe 2004; Gurtin *et al.* 2007) provide the physical motivation for size-dependent hardening models. However, a generally accepted correlation between the gradients of the variables and the hardening behaviour of the crystal could not be established, yet. Several thermodynamic gradient theories for polycrystals have been proposed associating the gradients of plastic quantities with focus on an increase of the free energy (Fleck *et al.* 1994; Steinmann 1996; Menzel & Steinmann 2000) or dissipation (Gurtin & Anand 2005; Fleck & Willis 2009). In the case of single crystals, a refined kinematical theory allows to compute geometrically necessary screw and edge dislocation densities for the different slip systems, closely related to Nye's tensor (Gurtin 2002), which has been generalized to the context of large deformations by Gurtin (2006). Associated thermodynamic theories have been proposed by e.g. Cermelli & Gurtin (2002), Berdichevsky (2006), Ohno & Okumura (2007) or Ekh *et al.* (2007). The dislocation field theory of Evers *et al.* (2004), see also Geers *et al.* (2006) and its generalization (Bayley *et al.* 2006) have been compared with thermodynamical approaches by Ertürk *et al.* (2009) and Bargmann *et al.* (2010). Similar non-local theories are, for example, based on incompatibility-dependent hardening moduli (Acharya & Bassani 2000) or other non-work-conjugate formulations (Kuroda & Tvergaard 2008). In the micromorphic approach of Forest (2009), a local inelastic variable is energetically coupled to an additional degree of freedom. A strong energetical coupling can formally be interpreted as penalty approximation of, e.g. a strain gradient plasticity theory. However, the coupling parameter allows for a specific adjustment of the scaling behaviour of the model, which was investigated by Cordero *et al.* (2010) and Aslan *et al.* (2011) for single crystal laminates and by Cordero *et al.* (2012) for periodic single crystals. Based on a closely related large-deformation theory for isotropic elastic–plastic materials, Anand *et al.* (2011) regularized strain-softening phenomena by means of the gradient of an additional variable related to the equivalent plastic strain. Hochrainer (2006) generalized Nye's concept to a higher-dimensional continuum dislocation theory, which was implemented by Sandfeld *et al.* (2010).

The fact that size-dependent theories usually require the computation of additional gradients makes the numerical treatment by the finite-element method more complex and computationally expensive (compared with local theories) since the variables have to be approximated by (at least piecewise) differentiable shape functions. One of the first works on this issue was published by de Borst & Mühlhaus (1992). Shu *et al.* (1999) presented and compared different element formulations for the Toupin–Mindlin framework of strain gradient theory. They introduced the displacement gradients as additional nodal degrees of freedom (associated with *C*^{0}-continuous shape functions) and coupled them to the true gradients by Lagrange multipliers. Liebe & Steinmann (2001) proposed a thermodynamic framework, where the yield condition is treated in a weak sense, leading to a distinction between plastically active and inactive nodes. A multi-field incremental variational framework for gradient-extended standard dissipative solids, including a generalization of this framework, can be found in Miehe (2011). Becker (2006) and Han *et al.* (2007) propose a projection of the plastic variables from the integration points to the nodes which subsequently allows the computation of plastic strain gradients and an implementation close to local theories. Forest (2009) proposes a generalized micromorphic approach with a local variable and a micromorphic variable. In addition, Wieners & Wohlmuth (2011) propose a primal–dual finite-element approximation. Both approaches preserve a local evaluation of the yield function and facilitate the application of the radial return algorithm.

One drawback of gradient plasticity theories is that a unique and generally accepted correlation between the gradients of the variables and the hardening behaviour is missing. Another disadvantage of implementations, including the plastic variables as additional nodal variables, is the significant increase in the computational effort owing to the extended number of degrees of freedom, especially in single crystal computations. Consequently, the application of the theory to three-dimensional problems is numerically extraordinarily expensive.

The present work provides a gradient theory for crystals which leads to efficient numerics. The formulation and model behaviour are close to existing gradient theories. We choose an equivalent plastic strain measure (which unites the plastic strain history of all slip systems in one scalar quantity) for the gradient enhancement aiming at a significant decrease of additional nodal degrees of freedom compared with existing theories. One central motivation for this proceeding is that equivalent plastic strain measures have been used successfully in the context of phenomenological local hardening models. Most of the above-mentioned works are based on geometrically necessary dislocations (GND), given by Nye's tensor which has a clearer physical interpretation than the gradient of the equivalent plastic strain. However, the present theory allows the simulation of three-dimensional systems (consisting of several grains) at reduced computational cost, owing to the significant reduction of degrees of freedom. The close relation to other gradient theories will be shown, and as long as the exact dependence of the hardening behaviour on geometrically necessary dislocation density measures (if there is any) remains unclear, even more complex gradient theories remain imprecise.

The outline of the paper is as follows: first, we discuss the basic kinematical and energetic assumptions. Subsequently, two different but physically identical variational settings are introduced. Details of the finite-element implementation are discussed before some numerical results illustrate the performance of the alternative formulation in three-dimensional applications. Finally, we analyse the ability of the theory to reproduce the experimentally observed size effects.

*Notation.* A direct tensor notation is preferred throughout the text. Vectors and second-order tensors are denoted by bold letters, e.g. ** a** or

**. A linear mapping of second-order tensors by a fourth-order tensor is written as . The scalar product and the dyadic product are denoted, e.g. by**

*A***⋅**

*A***and**

*B***⊗**

*A***, respectively. The composition of two second-order tensors is formulated by**

*B***. Completely symmetric and traceless tensors are designated by a prime, e.g.**

*AB***′. Matrices are denoted by a hat, e.g. . The transpose and the inverse of a matrix are indicated by and , respectively.**

*A*## 2. Equivalent plastic strain gradient enhancement of single crystal plasticity

### (a) Constitutive equations

#### (i) Kinematics

We describe the motion of a body by the displacement field ** u**(

**,**

*x**t*), which maps the positions of the particles to their associated displacements

**at time**

*u**t*. Here, is open with boundary and closure . We restrict ourselves to the small deformation context, i.e. ∥

**∥≪1, where**

*H***:=∇**

*H***=∂**

*u**u*

_{i}/∂

*x*

_{j}

*e*_{i}⊗

*e*_{j}. Furthermore, we define the strain tensor by

**:=sym(**

*ε***)=∇**

*H*^{S}

**and assume the classical relation 2.1for the plastic part of the displacement gradient, where the scalars**

*u**λ*

_{α}are slip parameters and

*N*is the number of slip systems which are characterized by slip plane normals

*n*_{α}and slip directions

*d*_{α}(orthogonal to

*n*_{α}) of unit length. The symmetric part of

*H*^{p}is the plastic strain 2.2The Schmid-tensor

*M*^{S}

_{α}of slip system

*α*is defined by . In order to capture the actual single crystals kinematics, the slip systems are defined pairwise with common slip plane normal and opposite slip directions. This implies for octahedral slip systems, e.g.

*N*=24. The slip parameters are constrained to increase monotonously, i.e. with initial values equal to zero. The effective plastic slips (in each slip system) are given by the difference between two pairwise defined slip parameters. It should be noted that alternatively 12 slip systems with possibly positive and negative slip increments for face-centred cubic (FCC)-crystals can be introduced. For most gradient theories, this leads to a reduced number of degrees of freedom in the finite-element implementation. The elastic strain is defined by

*ε*^{e}:=

**−**

*ε*

*ε*^{p}. We denote the plastic domain of a slip system

*α*by (where ).

Additionally, the equivalent plastic strain measure is given by the expression
2.3which maps the slip parameters to the equivalent plastic strain *γ*_{eq}. Recall that denotes the use of matrix notation. Note that, we choose definition (2.3) for simplicity and that other definitions of are possible. The theory can easily be generalized to more complex (e.g. tensorial) quantities like the plastic part of the displacement gradient . This generalization allows the application of the theory presented in the following to more complex models, which assume the free energy to depend, for example, on Nye's dislocation density tensor.

#### (ii) Free energy and dissipation

We assume the free energy per unit volume to be of the additive form (formally similar to, e.g. Steinmann 1996)
2.4with , where is the elastic stiffness tensor. The second part phenomenologically accounts for the isotopic hardening behaviour of crystals. Physically, this observation is primarily attributed to the multiplication of dislocations, which are mutually acting as obstacles leading to an increase of the macroscopically observed yield stress with proceeding plastic deformation. The third part *W*_{g}(∇*γ*_{eq}) introduces an internal length into the theory, i.e. it leads to a size-dependent mechanical response of the model. In the case of metals, size effects are observed experimentally for, e.g. micro specimens, steels with varying grain sizes or nanoindentation tests. The size effect can partially be explained by geometrically necessary dislocation configurations represented by Nye's tensor, which is closely related to the gradients of the plastic slips. The dissipation 𝒟 (dissipated power per unit volume) is assumed to be a superposition of the dissipation contributions of the individual slip systems
2.5where are dissipative forces which extend power over the plastic slip rates (cf. Cermelli & Gurtin 2002).

### (b) Variational form and strong form of the field equations

The basis for the subsequent theoretical development and for the numerical implementation is the principle of virtual power, which states that for a given state, described by the primary variables , the virtual power of the external forces is equal to the virtual power of the internal forces. The formal treatment is similar to that of Gurtin *et al.* (2007). The external forces are assumed to be given by the traction field on the associated Neumann boundary . The tractions generate power over . The micro-tractions are defined analogically, i.e. generates power over on the Neumann boundary . The primary variables are assumed to be known on the Dirichlet boundary and , respectively. Here, we restrict ourselves to micro-hard Dirichlet boundary conditions .

We choose ** v*** and to be the independent virtual kinematical quantities and require both of them to respect the Dirichlet boundary conditions on and . The dependent virtual rates can be expressed by
2.6The principle of virtual power states
2.7with
2.8and the virtual power of the dissipative forces
2.9Expression (2.8) motivates the introduction of the micro-stress
2.10which is energetically conjugate to ∇

*γ*

_{eq}.

Based on these results, equation (2.7) can be expressed in terms of the primary variables 2.11Using Gauss' theorem and the chain rule, we obtain 2.12 2.13 For vanishing , this expression yields 2.14

where the Cauchy stress has been introduced. Since ** v*** is arbitrary the integrands of both integrals must vanish independently and point-wise. The first integral in equation (2.14) implies
2.15
The second integral in equation (2.14) yields the Neumann boundary conditions
2.16
Choosing now

***=**

*v***0**and , we find 2.17 Introducing the resolved shear stress , equation (2.17) yields the field equations 2.18 Note that the contribution is equal for all slip systems due to the special choice of the scalar equivalent plastic strain

*γ*

_{eq}. For the Neumann boundary conditions, we obtain 2.19 Intending to use a visco-plastic overstress constitutive formulation, we define the yield criteria (cf. Miehe 2011) by 2.20 with the initial yield stress

*τ*

^{c}

_{0}. From equation (2.18), it follows that 2.21 Rate-independent single crystal plasticity theories suffer from the problem of possibly occurring linearly dependent constraints (Kocks 1970). In order to circumvent this problem, we base our formulation on a visco-plastic constitutive flow rule of the Perzyna-type (Perzyna 1971) 2.22 where

*g*(

*x*) is monotone with

*x*≤0⇔

*g*(

*x*)=0 and

*η*is a viscosity parameter (cf. Simo & Hughes (1998) or Miehe (2011) in the context of gradient theories). We define the viscous stress by , i.e. in the case of plastic yielding, and find from equation (2.21) the shear stress equilibrium in the plastic zones 2.23 Note that the negative of can formally be interpreted as a back-stress (cf. Gurtin

*et al.*2007), which is slip system independent, since

*W*

_{g}depends only on the gradient of the equivalent plastic strain.

### Equivalent plastic strain gradient theory

From the equations (2.5), (2.18) and (2.23), the dissipation can be deduced
2.24
This result formally corresponds to the result of the local theory (i.e. *W*_{g}=0). Hence, compared with a purely local theory, the dissipative mechanisms due to plastic slip rates are not altered in the theory at hand. In the case of plastic yielding, the dissipative stress *τ*^{d}_{α}, as obtained in equation (2.18), is
2.25
The introduction of a dissipation potential (here discussed only for linear viscosity), e.g. by
2.26
allows to construct a global incremental potential of the model (cf. Miehe 2011). However, this procedure is not part of the present work. The model equations introduced in this section are summarized in box 1.

### (c) Alternative formulation

In the following, we define a model in terms of four primary variables which is (concerning its mechanical behaviour) physically equivalent to the model introduced in §2*b*. We introduce an additional field variable *ζ* and weakly enforce its equivalence to *γ*_{eq} by a Lagrange parameter . This apparently complex approach leads to a massive reduction of the computational effort when solving the problem by use of the finite-element method. This idea is inspired by the work of Simo *et al.* (1985) (also by Forest 2009), who introduced the weak enforcement of the equality of two (*a priori* different) field variables by a force-like field which is treated like an additional field variable (see also Shu *et al.* 1999).

The free energy is assumed to take the form
2.27
We define ** ξ**:=∂

*W*

_{g}/∂(∇

*ζ*) for the subsequent developments and start from the weak form 2.28 Application of Gauss' theorem and the chain rule yields 2.29 The field equations and Neumann boundary conditions associated with (2.29) are summarized in box 2. Additionally, the yield criteria and flow rule are listed. It should be noted that the slip parameters can alternatively be handled as internal variables. The associated evolution rules can then be derived based on the reduced dissipation inequality (for further details, see Forest (2009)).

Similar to the case before (equation (2.23)), we obtain the relation 2.30 for the plastic zones.

Note that formally, can be interpreted as a back stress again and ∂*W*_{h}/∂*λ*_{α} is a local hardening stress. The back stress is equivalent for all slip systems due to the reduction to only one scalar variable *ζ* and the special choice of .

From the set of equations in box 2, the equations in box 1 can be derived. In this sense, the models (2.11) and (2.28) can be considered to be equivalent.

### Equivalent plastic strain gradient theory, alternative formulation

The original system of equations (box 1) includes three scalar partial differential equations for the linear momentum balance and *N* (e.g. *N*=24/12 in the case of FCC-crystals) partial differential equations for the micro force balances, whereas the total number of scalar partial differential equations of the new set of equations (box 2) is only four. It is emphasized that this massively reduced degree of non-locality makes the new set of equations (box 2) interesting for a numerical implementation via finite elements. In contrast, the discretization of the first model leads to elaborate element stiffness matrices and active set search algorithms on the nodal level which have barely been studied for crystals, yet. Therefore, only the second model is implemented. The numerical implementation is presented in the following.

## 3. Numerical implementation

In this section, the local (integration point) algorithm, the algorithmic tangent moduli and the element stiffness matrices for the alternative formulation are derived. It is shown that the number of degrees of freedom per node increases only by one compared with a local theory. For the numerical implementation, explicit expressions for the free energy are required. We assume quadratic forms for simplicity
3.1
with constant hardening moduli *h*_{αβ} and *K*_{g}.

The starting point for the discretization via finite elements is the time-discretized problem (2.28), evaluated at time step *t*_{n+1} (for convenience, the index ‘*n*+1’ is dropped in the following), given by the residuals associated with the variations of the displacements
3.2
of the variable *ζ*
3.3
of the back stress
3.4
and of the slip parameters *λ*_{α}
3.5
Here,
3.6
is the slip rate which is assumed constant during the time interval [*t*_{n},*t*]. For simplicity, we prescribe a linear-viscous behaviour
3.7

### (a) Local algorithm

The structure of the set of residuals allows for an algorithmically decoupled determination of the quantities (** u**,

*ζ*) and the local quantities , which is close to standard iterative procedures in computational inelasticity, solving the discretized momentum equations and updates of the local variables separately (e.g. Simo & Hughes (1998); Miehe & Schröder (2001), among many others). For the computation of the local quantities, we assume the solution (

**,**

*u**ζ*) to be given. The local problem, given by equations (3.4) and (3.5), is constrained by the requirement Δ

*λ*

_{α}≥0, which will be exploited to identify the local set of active slip systems (containing the indices of the active systems). For the moment, assume an active set to be given. The aim here is the development of an algorithm which allows the computation of 3.8 where

*n*

_{act}is the number of active slip systems. Recall that denotes the use of matrix notation. The subscripts in the parenthesis and the superscript ‘S’ account for the active slip parameters.

To identify the solution, we compute the residuals (3.4) and (3.5) and therefore use the following matrix notation 3.9 based on the state . The solution can then be obtained as the solution of the linear system of equations with the symmetric tangent 3.10 The explicit elements of are given by 3.11 For convenience, the complete local algorithm is summarized in box 3 containing a standard active set search procedure (see Ortiz & Stainier (1999); Miehe & Schröder (2001), and references therein).

### (b) Algorithmic tangent moduli

The local update procedure is coupled to the global algorithm via a generalized form of the classical algorithmic tangent moduli
3.12
In order to compute the algorithmic tangent, consider the definition The basis for the computation of the algorithmic tangent moduli is the requirement, that the residuals remain zero, if the quantities (** u**,

*ζ*) change 3.13 Accordingly, the differential of the local variables can be reconstructed explicitly as a function of the differential of the quantities 3.14 Hence, the tangent operator , containing the partial derivatives 3.15 can be computed as a matrix product and allows to identify two of the required algorithmic tangent moduli. The matrix is given by the expression 3.16 The identification of the remaining algorithmic moduli and can be effectuated based on the upper entries of , which allow the direct computation of the increment of the plastic parameters 3.17 Finally, the computation of the stress-increment 3.18 leads to the identification of the algorithmic tangents and . It should be noted that due to the symmetry of the problem, the equivalence can be shown.

### Local algorithm

In the case of an elastic update, i.e. *f*_{α}<0 ∀ *α*∈{1,2,…,*N*}, we set
3.19
Here, *C*^{pζ}_{el} is a large (negative) number which penalizes a deviation of *ζ* from *γ*_{eq}. Alternatively, the equation , evaluated at the elastic integration point, can be interpreted as kinematical constraint on the nodal values associated with *ζ* of the corresponding element.

### (c) Element stiffness matrices

With the algorithmic moduli at hand, the element stiffness matrices associated with the linearized global residuals *G*^{u} and *G*^{ζ} can be computed. The linearization of *G*^{u} with respect to ** u** yields
3.20
Here, is the weight due to the numerical quadrature scheme associated with integration point p of element e. The matrix is the standard interpolation operator, projecting the nodal displacements to the strain at the integration points.

The linearization with respect to *ζ* is incorporated by
3.21
where interpolates the nodal values of *ζ* to the integration points. The symmetry of the problem was exploited again in order to compute . Finally, the linearization of *G*^{ζ} with respect to *ζ* leads to
3.22
where interpolates the nodal values of *ζ* to its gradient in the integration points. In the finite-element procedure, the matrices , , and represent parts of the (symmetric) element stiffness matrix.

## 4. Numerical examples

We demonstrate the performance of the implementation of the alternative formulation by the simulation results of micro-mechanical tension and torsion tests of micro-components which have been carried out with the in-house finite-element program of the Institute of Engineering Mechanics (Continuum Mechanics).

In both cases, we use an oligo-crystal consisting of eight grains with a simplified grain geometry as shown in figure 1 and random crystal orientations. The hardening parameters are assumed to be given by *h*_{αβ}=*H*[*q*+(1−*q*)*δ*^{αβ}] (Hutchinson 1970) with *q*∈[1;1.4] as suggested by Kocks (1970). We use the following material parameters

C_{11} | C_{12} | C_{12} | K_{g} | H | q | τ_{0}^{c} | η |
---|---|---|---|---|---|---|---|

168GPa | 121GPa | 75GPa | 3 × 10^{−2}N | 200MPa | 1.1 | 35MPa | 0.001MPa s |

The boundary conditions for the tension tests allow for free lateral contraction. We choose micro-free boundary conditions for *ζ*, i.e. . All simulations were performed with standard linear hexahedrons and adaptive time steps. The total simulation time is 1 s. Figure 2 shows the results for the tensile and torsion test simulations for varying cube sizes. The model response for torsion is normalized in order to allow a comparison of the responses of the cubes with different dimensions. The results (figure 2) show that the size-dependence of the mechanical response is stronger in the case of torsion (compared with tension). This difference, which is also observed in experiments (Fleck *et al.* 1994), occurs (here) due to increased strain gradients in the case of torsion, which are penalized by *W*_{g}(∇*ζ*). On the integration point level, this effect is represented by the back stress (cf. equation (2.30)). Figure 1(right) illustrates that can be considered as a hindrance for further plastic deformation in strongly deformed regions while it augments plastic processes in little distorted parts of the body and thereby reduces gradients of *ζ*. For both load cases, the simulations show a strong size effect for . The negligible size effect in the case of larger bodies is due to the fact that the strain gradients decrease with increasing dimensions, i.e. the theory at hand tends to a classical local theory for large bodies with small gradients. An idealized (e.g. misorientation-induced) boundary slip resistance can be introduced into the model by micro-hard boundary conditions. For the simulation of the micro-component (cf. figures 1 and 3), we model the internal grain boundaries (represented by *Γ*) as micro-hard, i.e. *ζ*=0 ∀ ** x**∈

*Γ*. Formally, this case is covered by the presented theory by considering the oligo-crystal as a union of eight bodies. The implementation is effectuated by setting the prescribed nodal values of

*ζ*to zero and eliminating the associated lines and columns from the system matrix and the residual vector. The non-local hardening modulus

*K*

_{g}has not been changed to allow a comparison with the micro-free boundary conditions (figure 2). However, for more realistic applications,

*K*

_{g}should be reduced in the case of micro-hard grain boundaries, as the response shown in figure 3 is much too hard compared with experimental results. The mesh dependence of the results is depicted in figure 4. The study was performed with a cube size of

*L*=30 μm. The choice of different cube sizes resulted in comparable results. The convergence study suggests that, in the case of the oligo-crystal with simplified grain geometry, the mesh with approx. 37 000 d.f. (20×20×20 elements), which has been used in figures 2 and 3, yields reasonable global results close to the converged solution (especially for micro-free boundary conditions). A more realistic (and therefore more complex) grain geometry is expected to require an increased number of degrees of freedom per grain. Additionally, a mesh dependence study with 512 crystals (8×8×8 crystals) was carried out for micro-free boundary conditions and

*L*=30 μm (figure 4) to get a first impression of the mechanical response of a polycrystal.

Table 1 shows the Euclidean norm of the residual for a tensile test simulation with micro-hard boundary conditions, a cube size of *L*=30 μm and a 8×8×8-mesh. For time steps 3 and 4 (elasto-plastic transition with 16/13 iterations) not all values are represented. The total CPU time for the solution was 25.89 *s* (table 1) on a Pentium Dual Core PC with 3.0 *GHz* and 6 *GB* RAM (which limits the maximum number of d.f., in the case of a direct solver, to approx. 300 000).

## 5. Summary

A geometrically linear single crystal model extended by an energetic hardening term based on the gradient of the equivalent plastic strain *γ*_{eq} has been introduced. The direct finite-element implementation of the theory is computationally expensive due to the need of piecewise differentiable shape functions for the plastic slip parameters *λ*_{α}. The introduction of an alternative (but physically equivalent) formulation based on a Lagrange multiplier and a micromorphic-like variable *ζ* drastically reduces the number of partial differential equations and nodal degrees of freedom of the finite-element implementation. The local algorithm and its linearization, represented by the extended algorithmic tangent moduli, have been discussed in detail. Numerical simulations of three-dimensional tension and torsion tests of oligo-crystals with a simplified geometry illustrate the performance of the implementation. The overall model response qualitatively mimics the size-dependent material behaviour observed in experiments (Fleck *et al.* 1994; Stölken & Evans 1998; Xiang & Vlassak 2006; Gruber *et al.* 2008). The experimentally found size-dependence of the yield stress is not captured, a well-known consequence of the simple quadratic ansatz for the strain gradient-dependent energy contribution, which can partially be remedied by more elaborate approaches (Ohno & Okumura 2007).

Overall, the proposed geometrically linear model allows computationally relatively cheap three-dimensional size-dependent simulations of oligo-crystal models.

The theory is not restricted to the special scalar variable *γ*_{eq}. Instead, the presented concepts can be generalized to more complex (e.g. tensorial) quantities like the full plastic part of the displacement gradient and thereby cover theories taking, e.g. Nye's tensor as argument of the free energy.

The boundary conditions for the field *ζ* remain an issue for further research. The physical relevance of micro-free and micro-hard boundary conditions at the grain boundaries is questionable. Both conditions seem to be applicable in order to reproduce the global size-dependent behaviour. However, their validity is restricted to certain ranges of the deformation and size. Therefore, more elaborated grain boundary models should be examined.

## Acknowledgements

The authors acknowledge the support rendered by the German Research Foundation (DFG) under grant BO 1466/5-1. The funded project ‘Dislocation-based gradient plasticity theory’ is part of the DFG Research Group 1650 ‘Dislocation-based plasticity’.

- Received February 7, 2012.
- Accepted March 21, 2012.

- This journal is © 2012 The Royal Society