The Role of Geometrically Necessary Dislocations in Cantilever Beam Bending Experiments of Single Crystals

The mechanical behavior of single crystalline, micro-sized copper is investigated in the context of cantilever beam bending experiments. Particular focus is on the role of geometrically necessary dislocations (GNDs) during bending-dominated load conditions and their impact on the characteristic bending size effect. Three different sample sizes are considered in this work with main variation in thickness. A gradient extended crystal plasticity model is presented and applied in a three-dimensional finite-element (FE) framework considering slip system-based edge and screw components of the dislocation density vector. The underlying mathematical model contains non-standard evolution equations for GNDs, crystal-specific interaction relations, and higher-order boundary conditions. Moreover, two element formulations are examined and compared with respect to size-independent as well as size-dependent bending behavior. The first formulation is based on a linear interpolation of the displacement and the GND density field together with a full integration scheme whereas the second is based on a mixed interpolation scheme. While the GND density fields are treated equivalently, the displacement field is interpolated quadratically in combination with a reduced integration scheme. Computational results indicate that GND storage in small cantilever beams strongly influences the evolution of statistically stored dislocations (SSDs) and, hence, the distribution of the total dislocation density. As a particular example, the mechanical bending behavior in the case of a physically motivated limitation of GND storage is studied. The resulting impact on the mechanical bending response as well as on the predicted size effect is analyzed. Obtained results are discussed and related to experimental findings from the literature.


Introduction
Micromechanical testing of small-scaled single crystals has been excessively practiced in the last two decades to study the mechanical size-dependence of diverse materials [1][2][3][4]. Different intrinsic (microstructural) effects have been found to be triggered by the interplay of physical size limitation such as free surfaces and the underlying microstructure in which the initial density of dislocations plays a crucial role. For example, experiments and 3D-discrete dislocation dynamics (DDD) simulations indicated that the size-dependent response of Ni single crystals decreases with increasing starting dislocation density [5]. Further, stable plastic deformation for crystals of "360 nm size was achieved and no strengthening effect in the range of 360 nm-1500 nm was observed for Mo alloy microcrystals [6]. In both cases, the crystals were machined from relatively strong pre-strained bulk. Conversely, a sufficiently low initial dislocation density provokes a rapid starvation of available dislocations [7]. However, this effect becomes dominant if, for instance, the diameter D of micropillar samples falls below a critical value ("1 µm in [7]). In the study of Shan et al. [8], an initial dislocation density of "10 15 m´2 has fallen very fast to zero for D " 160 nm whereas already larger samples (D ě 250 nm) were less likely to be dislocation-free after testing. Besides the characteristic sample dimension, the overall crystal size may additionally be considered as suggested recently by El-Awady [9].
With increasing sample size, other mechanisms will superimpose and start to affect or even govern the mechanical behavior. Depending on the geometry of deformation induced into the sample in terms of the applied loading, plastic strain gradients and associated geometrically necessary dislocations therewith, may increase significantly the plastic work hardening. Bending represents a typical deformation mode in which plastic strain gradients are induced extrinsically. Fleck and Hutchinson [10] suggest that the stored density of GNDs is proportional to the curvature of the bended beam. Experiments have found a strong inverse correlation between beam thickness and increase in strength of the material. For instance, an increase in the work hardening behavior with decreasing beam thickness has been reported by Stölken and Evans [11] for thin beams with a variation in thickness between 12.5 µm and 100 µm. Microbending experiments of Motz et al. [12] confirmed the strong correlation between flow stress and beam thickness. The investigated beam sizes varied in thickness between 7.5 µm and 1 µm. In the same work, it was pointed out that a regular alignment of GNDs (as would be expected in pure bending) is insufficient to explain the size effect. The responsible mechanism could be dislocation source limitation as a result of a rapid starvation of initially available dislocations as well as back-stress effects induced by dislocation pile-ups along the neutral axis. Numerical results of 3D-DDD beam bending simulations with thicknesses between 0.5 µm and 1.5 µm revealed that a combination of pile-up of GNDs and source size limitation mimics the experimental data quiet well [13]. These findings indicate the complexity of size-dependent strengthening behavior due to superposition of different mechanisms. Accordingly, it is difficult to independently estimate the size-dependent hardening contribution associated with each mechanism based on experimental data. For that reason, a gradient extended crystal plasticity model is used to investigate the size range (characterized by means of the beam thickness) and the extent to which GNDs affect the size-dependent bending behavior of small-scaled cantilever beams.
The aim of this work is to examine the size-dependent bending response for a range of sample sizes. In particular, we focus on three sample sizes with a variation in thickness between 2.5 µm and 5.0 µm. For this purpose, a higher-order gradient crystal plasticity model is implemented in a three-dimensional finite element framework. Highlights of the model are non-standard evolution relations for edge and screw components of the slip system-based dislocation density vector, crystal-specific interaction relations, and higher-order gradient boundary conditions. The gradient effect associated with accumulating GND densities is of primary interest. In order to gain additional insights, the size-dependent impact of GNDs on the evolution of SSDs is investigated as well. The mechanical response of cantilever beams is further addressed for the case when the evolution of GNDs saturates due to a physical limitation. The number of dislocations to be stored locally due to compatibility reasons cannot be arbitrary high. Accordingly, a feasible limit for GND storage is introduced. All results are discussed and related to experimental findings from the literature.

Higher-Order Gradient Crystal Plasticity Theory
The large deformation-based, gradient-enhanced crystal plasticity theory of Bargmann et al. [14] is followed, see also [15]. As part of this, the multiplicative decomposition of the deformation gradient F " F E¨F P into an elastic F E and a plastic part F P is adopted. The state space is specified by s " pC E , tγ α u, t∇ i γ α uq, in which C E " F T E¨F E denotes the elastic right Cauchy-Green strain tensor, tγ α u represents a set of plastic slip variables, and t∇ i γ α u is the corresponding set of plastic slip gradients. Both sets scale with the number of admitted slip systems, where α denotes a particular slip system. In this framework, the plastic slip γ α is used as a primary measure for plastic deformation Although the plastic slip rate ν α may not necessarily represent the material-time derivative of γ α , such an approximation is often applied in practice. This assumption is adopted in the following.
We employ a thermodynamic consistent formulation in which the free energy storage with respect to the intermediate configuration is expressed by the following additive split such that ψ e i pC E q represents a hyperelastic contribution, ψ l i ptγ α uq covers the local energy storage due to dislocation glide, and ψ g i pt∇ i γ α uq is the gradient related energy storage due to accumulation of GNDs. The affiliation of a physical quantity ‚ is indicated by the following notation: ‚ r (reference space), ‚ i (intermediate space), ‚ c (current space). This notation is also applied to differential operators, e.g., ∇ i ‚, Div i p‚q, and Curl i p‚q.
Under quasi-static, isothermal conditions, the dissipation inequality with respect to the intermediate configuration reads In this, the stress power density P " P E`PP is decomposed into an elastic part P E " 1{r2ρ i sS E : 9 C E , based on the second Piola-Kirchhoff stress tensor S E , and a plastic part P P " 1{ρ i M E : L P , in which M E " C E¨SE is the Mandel stress tensor and L P " 9 F P¨F´1 P " ř α ν α rs α b n α s denotes the plastic velocity gradient defined as a superposition of individual plastic slip rates on corresponding slip systems. As usual, s α , n α , and t α " n αˆsα are the slip direction, slip plane normal, and transverse direction vectors respectively. The intermediate configuration is assumed to be isoclinic (Any set of orthonormal vectors ts α , t α , n α u which describes the crystallographic geometry of a slip system does not change between the reference and the intermediate configuration [16]. Hence plastic deformation is lattice (volume) preserving, i.e., detpF P q " 1.). Further, Υ i represents the specific energy storage rate Following standard thermodynamic arguments, the constitutive relation for the hyperelastic state results in S E " 2 Bψ e i {BC E . In accordance to Ekh et al. [17], we introduce the scalar-valued micro-stress κ α " Bψ l i {Bγ α and the back-stress vector (or vector-valued micro-stress) κ α " Bψ g i {B∇ i γ α work conjugate to γ α and ∇ i γ α respectively. Substituting these constitutive relations into Equation (3) leads to the reduced dissipation inequality Integration by parts and applying Gauss's theorem gives A more restrictive validity of the dissipation inequality is obtained by its local form. This gives with respect to the bulk dissipation and with respect to the boundary dissipation where τ α " M E : rs α b n α s represents the resolved shear stress and κ i denotes the micro-hardening stress projected on the outward pointing normal vector N (b) i -here associated with an infinitesimal intermediate surface element. The relation between the referential surface normal vector and the intermediate surface normal vector is given by the cofactor of the plastic deformation, i.e., N i dA i " cofpF P qN r dA r . The current rate of dissipating work depends on the history of plastic deformation in terms of γ α but also on the current plastic slip rate ν α , cf. Equations (7) and (8). Moreover, for homogeneous plastic deformation, i.e., κ α " 0, Equation (7) reduces to D L red " ř α rτ α´κα s ν α ě 0 whereas the boundary contribution in Equation (8) completely vanishes.
Then, the onset of plastic yielding is regulated by means of a yield function φ α for each slip system, here defined as where Y α defines the initiation of plastic yield in the absence of hardening, i.e., the initial slip system resistance which is equivalent to the critical resolved shear stress of the system, whereas S α " κ α´D iv i pκ α q`Y α is the current slip system resistance.
In the context of viscoplasticity, overstress states τ α´Sα ą 0 are generally allowed and typically regularized via a power law relation. Here, a Perzyna form is chosen as a viscoplastic regularization: The brackets define a ramp function of the form ă φ α ą" 1{2 rφ α`| φ α |s and ensure that ν α ‰ 0 only for τ α ą S α . Moreover, m denotes the rate sensitivity exponent of the stress ratio, ν 0 represents the reference shear rate, and C 0 is the drag stress.

Governing Equations
The so-called dislocation tensor represents a continuum measurement for GND densities from which the vector-valued GND vector 9 g iα is derived, cf. [18]. In extension to [14,15], the evolution of 9 g iα -here with respect to the intermediate configuration-is taken as where b α denotes the Burgers vector. This relation accounts for dislocation collision processes, for instance, interactions between mobile dislocations of active slip systems mimicked by the plastic slip rates and stored dislocations on latent slip systems. For most applications, it is sufficient to account for the edge and screw component of the GND density vector. It is assumed that there is no development of plastic slip gradients perpendicular to the slip planes, i.e., 9 g iαˆn α " 0 [19]. With respect to fcc crystals, modified interaction relations are proposed for the edge and screw components of the GND density based on moduli for edge-edge ι ee αβ , edge-screw ι es αβ , screw-screw ι ss αβ , and screw-edge ι se αβ dislocation intersections: ι ee αβ "ˇˇs α¨sβˇˇnαˆnβˇ, ι es αβ "ˇˇs α¨tβˇˇnαˆnβˇ, ι ss αβ "ˇˇt α¨tβˇˇnαˆnβˇ, ι se αβ "ˇˇt α¨sβˇˇnαˆnβˇ. (12) In addition, slip coplanarity moduli χ αβ are introduced in accordance to previous works (e.g., [20][21][22]) as In the case of coplanar slip systems, i.e., slip planes of system α and β are parallel to each other, all intersection moduli vanish and slip-system interactions are solely determined by coplanarity moduli. Recall that the GND components are obtained by projecting the plastic slip gradients on the slip directions s α and the transverse slip directions t α (see Arsenlis and Parks [23], Gurtin and Anand [24]). Then, under consideration of the above introduced interaction moduli, Equation (11) is substituted by the following two scalar-valued field equations which read and 9 g s iα " respectively. The first relation in Equation (14) accounts for the impact of stored edge and screw GND densities (with respect to latent slip systems β) on the evolution of the edge GND density of slip system α whereas the second term measures the geometrically necessary edge dislocation density by means of the plastic slip gradient ∇ i ν α . A similar relation is given for the geometrically necessary screw dislocation density. The negative sign of the second term in Equation (15) yields a convention in which the right-handed screw dislocation segment is treated as positive. An equivalent expression is found in the literature where a positive sign is used together with the line direction vector l α "´t α . For the sake of completeness, we recall the balance of linear momentum for quasi-static and isothermal conditions and in the absence of body forces 0 " Div i pF E¨SE q.
where P " F E¨SE¨F´T P is the relation between the first and second Piola-Kirchhoff stress tensor.

Constitutive Relations
The hyperelastic energy contribution is postulated in terms of the Neo-Hookean law where J E " detpF E q " detpC E q 1{2 is the elastic Jacobian determinant, I C E " trpC E q the first invariant of C E , and µ together with λ denote the Lamé parameters. Then, the hyperelastic stress response results in As observed in various investigations (cf. [12,25,26]), the stress-strain response of copper single crystal is characterized by a pronounced saturation behavior. In the study of Kleemola and Nieminen [27], it was found that Voce hardening [28] represents the best choice to describe the hardening behavior of Cu crystals. Therefore, the local energy contribution associated with dislocation glide (see [29]) reads c sat is the saturation rate parameter and ∆H l α " H l 0´Y α is the saturation hardening defined as the difference between local hardening modulus H l 0 and initial yield resistance Y α of the slip system. In this hardening law, latent hardening is addressed in terms of the accumulated plastic slip ř β γ β within the exponential function such that high slip system activity increases the effective local hardening contribution. The corresponding scalar-valued micro-stress κ l α is derived as Recalling that the free energy increases due to storage of GNDs in the material, edge and screw dislocation characters associated with distortion and twisting of the crystal lattice respectively, are introduced into the defect energy by the form H e 0 and H s 0 are the gradient hardening moduli related to the edge and screw component respectively, and l α is a constitutive internal length scale parameter. The relation between the gradient hardening moduli is taken as according to their elastic strain energy ratio (cf. [30]) where ν " λ{r2λ`2µs is the Poisson's ratio. Then, the constitutive relation for the back-stress vector is obtained from the chain rule Finally, due to the dependence of 9 g e iα and 9 g e iα on ν α , the scalar-valued micro-stress κ g α related to interaction processes is introduced via the relation In fact, κ g α describes a non-local hardening contribution based on GND intersection and collision effects which are particularly pronounced upon load reversal. In the end, the micro-stress κ α is comprised of a local as well as a non-local contribution, i.e., κ α " κ l α`κ g α .

Boundary Conditions
As we study the bending behavior of a single crystal with free surfaces, microfree boundaries are chosen. They refuse any dislocation pile-ups at the exterior of the crystal and, hence, boundaries appear to be transparent to dislocation motion This results in zero boundary dissipation according to Equation (8). An alternative approach based on non-idealized boundary conditions with non-zero boundary dissipation involving the effect of boundary yielding is presented in Husser et al. [30].

Numerical Implementation
The solution algorithm for the highly coupled and strongly non-linear multi-field problem is based on the dual-mixed finite element method as proposed by Ekh et al. [17]; see also Bargmann et al. [14,31]. In this, GND densities are introduced as nodal degrees of freedom in addition to the displacement. The basis for implementing the material model into a finite element framework is the variational form of the underlying governing equations. Applying the principle of virtual work to Equation (16) where δu is a vector-valued test function and δC E " rF T E¨∇i pδuq`∇ i pδuq T¨F E s denotes the variation of the right Cauchy-Green strain tensor. The corresponding mechanical boundary conditions read In a similar manner, the variational forms of Equations (14) and (15) are obtained. We further choose an implicit finite-difference method for the time discretization of the global field relations such that 9 g e iα " ∆g e iα {∆t, where ∆t measures the current time increment and ∆g e iα " g e iαpn`1q´g e iαpnq . The same holds for the screw component, i.e., ∆g s iα " g s iαpn`1q´g s iαpnq . As an approximation, plastic slip rates are discretized analogically, i.e., ν α " ∆γ α {∆t, where ∆γ α " γ αpn`1q´γαpnq . With this at hand, the variational forms are written as and 0 " respectively. Here, δg s α and δg e α are arbitrary test functions.
The time-discretized versions of the micro-stresses (Equations (23) and (24)) read Both equations are fully implicit. Therefore, it is not explicitly indicated that the time-dependent quantities are associated with the new time increment t pn`1q " t pnq`∆ t.

Finite Element Model
The presented gradient-based crystal plasticity model is applied to microbending experiments of copper single crystal. Three micron-sized cantilever beams with varying thickness t in the range between 2.5 µm and 5.0 µm are under investigation, see Table 1 for exact sample dimensions. In order to prevent additional influences on the size-dependent hardening, the momentum arm l b as well as the edge length in width direction w are kept constant for all geometries. This allows for a meaningful interpretation of the results for which a strong correlation between the strength of the material and the thickness of the beam is expected (Evans and Hutchinson [32]). Sample geometries and crystallographic orientation were exemplary chosen in accordance to the experimental set-up of Motz et al. [12]. Details regarding sample preparation and fabrication are provided in [33]. In that experimental study, various single crystalline cantilever beam samples were fabricated by the focused-ion beam (FIB) technique and loaded with an indenter tip at the free end. With regard to the finite element model, quadratic serendipity elements (twenty-node) have been used for geometry approximation whereas displacement and the GND density degrees of freedom have been solved with a different number of nodes, cf. Figure 1. Here, the solution of the displacement field is based on a fully quadratic FE-approximation combined with a reduced integration scheme (2ˆ2ˆ2 Gauss points). This approach is known to be well suitable for bending-dominated problems. In contrast, the GND densities are only evaluated at the corner nodes in terms of a linear FE-approximation using the full integration scheme, see Figure 1b. This mixed-element formulation-henceforth denoted as 20RI8FI-was examined by Kuroda [34] for the two-dimensional case with respect to simple shear and compression problems and was revealed to be well suitable for applications in the context of higher-order gradient crystal plasticity as it exhibits a reliable performance. The finite element meshes of the cantilever beams are illustrated in Figure 2 whereas the corresponding geometry and discretization data is provided in Table 1. Figure 2 includes a mesh for which both fields are approximated with trilinear eight-node elements and 2ˆ2ˆ2 Gauss points (full integration scheme), cf. Figure 1a. In the results section, a comparison in performance between the mixed-element (20RI8FI) and the fully linear formulation (8FI8FI) is carried out for selected cases concerning cantilever beam #1.  Table 1.

Crystallography and Material
The crystallographic orientation of the crystals is chosen in accordance to the experiments in [12]: the r110s direction is aligned parallel to the longitudinal beam axis (parallel to X 1 ) and the p111q-plane is oriented parallel to the X 1´X2 plane, cf. Figure 2a. The applied deflection results in the typical tensile and compressive dominated zones which determine the resolved shear stress on the individual slip systems and, thus, dominate their activation (at least for the here relevant bending load regime).
There exist four non-zero Schmid factors f α for the particular crystal orientation as indicated in Table 2 together with the slip system designation of the fcc lattice. Table 2. C.s.s. of copper (fcc): slip direction s α , transversal slip direction t α , slip plane normal n α , and initial Schmid factor f α . Slip activation is mainly determined by the bending stresses as an accommodation of plastic deformation by slip systems with an initially zero Schmid factor is rather unlikely. ?
The elasticity parameters for copper, i.e., Young's modulus E " 126.9 GPa and Poisson's ratio ν " 0.35, are taken from [12], respectively. Further, the initial yield limit is chosen to be Y α " 1.5 MPa in agreement with experiments of single crystals, see for instance [36,37]. The magnitude of the Burgers vector is taken as b " |b α | " a{ ? 2 " 0.2552 nm, based on a lattice constant of a " 0.3609 nm [38]. In order to minimize rate effects, the rate-sensitivity parameter is chosen to be m " 20 and the reference slip rate is put on a level with the macroscopic (quasi-static) load rate, i.e., ν 0 " 9 ε ben . The drag stress parameter is assumed to be C 0 " 10 MPa. Moreover, the same constitutive length-scale parameter l " l α is applied to all slip systems in analogy to the Burgers vector magnitude. In [11], l was found to be 4 µm for highly pure Ni. Since a similar magnitude was obtained for torsion tests of copper wires, cf. [39], this value is adopted here. The values for the hardening moduli and the saturation rate are discussed in Section 5.1. All material parameters are summarized in Table 3.

Element Choice: Eight-Node Hexahedron Element with Full Integration (8FI8FI) vs. Twenty-Node Brick Element with Mixed Integration (20RI8FI)
As a starting point, the size-independent response is studied, i.e., the bending response in the absence of gradient effects which is obtained when l{t « 0. By setting l " 0, the computations mimic the response of bulk samples independent of the actual sample dimensions. A macroscopic bending test from [12] serves as a reference in order to calibrate the local hardening modulus H l 0 as well as the saturation rate c sat . All numerical bending tests were loaded up to 10% normalized deflection (applied deflection/initial momentum arm). Optimal values are identified as H l 0 " 77.5 MPa and c sat " 10 3 . The final results are shown in Figure 3. As seen, this choice results in a good saturation behavior with rapid hardening behavior for all three sample sizes associated with the 20RI8FI-formulation. All curves yield the reference flow stress of «227 MPa, i.e., the bending response is clearly size-independent (independent of the beam thickness).  Simulations with all three sample sizes (20RI8FI) yield the reference data of [12]. The FE-mesh with linear elements (8FI8FI) overestimates the stiffness as well as the strength due to volumetic locking in case of bending-dominated loading for the size-independent as well as the size-dependent case.
By comparing the bending response of cantilever beam #1 of the two element formulations, it is clearly seen that the 8FI8FI-element formulation overestimates the strength and the bending stiffness. A deviation in the stress-strain curve is already obvious after «1% normalized deflection. Further, the saturation level as well as saturation behavior are not captured correctly. This is due to a bending-dominated deformation mode which cannot be properly captured by the linear element formulation due to locking effects.
Next, the size-dependent material behavior is investigated. Exemplarily, cantilever beam #1 is studied as the strongest impact is expected for the beam sample with the smallest thickness. A linear-like hardening behavior is observed (cf. Figure 3) which is associated with a continuously increasing (plastic) strain gradient during bending. The strain gradient scales the higher-order gradient hardening in terms of the back-stress Div i pκ α q. Hence, the gradient hardening contribution constantly increases which in the end prevents a saturation of the overall hardening response. As in the size-independent case, the 8FI8FI-formulation overestimates the bending response. Yet, the difference at the final deformation state appears to be not that pronounced which indicates that the computational accuracy of the GND field is not affected by the 8FI8FI-element formulation.
In the following, a qualitative comparison with regard to the GND density and the plastic slip is carried out. For that reason, the effective GND density g eff i is introduced as a function of the admitted edge and screw GND components where g tot iα is the total GND density associated with slip system α. In analogy, the effective plastic slip serves as a representative variable in order to compare the distribution of plastic deformation. Exemplarily, Figure 4 depicts the distribution of the effective GND density g eff i along the central middle axis for both element formulations. The results are within the same order of magnitude or even identical, i.e., the GND density field is well captured by both formulations. This holds true for both, the qualitative distribution along the middle central axis as well as the quantitative distribution within the X 1´X3 -plane. Consequently, the bending response (but not the GND evolution) is affected by volumetric locking in case of the 8FI8FI-element formulation. The contour plot for the case of 20RI8FI-elements is not presented as only marginal differences are present compared to the contour plot in Figure 4. Nevertheless, the contour plot of g eff i can be found in Figure 7 (Section 5.2) where the pile-up characteristic of GNDs is analyzed in more detail. The distribution of the effective plastic slip γ eff is qualitatively compared along the highest and lowest central path parallel to the beam axis, see Figure 5. In addition, a quantitative comparison is provided on the basis of the effective plastic slip distribution within the central X 1´X3 -plane. As seen, differences are found near the supporting end, i.e., close to a normalized position of « 0. Here, the 20RI8FI-element formulation resolves a higher magnitude of γ eff which is confirmed by the contour plots. Slightly higher values of γ eff are also computed in regions where both formulations are close to each other. Besides, it can be seen that plastic deformation is only accommodated within the first third/half of the beam sample (referred to the fixed side) whereas the rest of the beam finger remains straight. Such a localized plastic deformation is characteristic for cantilever beam bending.

Bending Size Effect-Influence of Sample Thickness
Experiments show that the bending size effect is strongly correlated to the beam thickness. For that reason, the ability of the model to predict a strengthening effect as a result of thickness reduction is studied. This investigation allows in turn to quantify the role of GNDs within cantilever beam bending experiments. All subsequent computations are based on the 20RI8FI-element formulation. All three sample sizes are loaded up to 10% normalized deflection using a gradient hardening modulus of H e 0 " 1 GPa, which corresponds to H s 0 " 650 MPa via Equation (22). A relatively high value is chosen in order to obtain a rather strong size effect such that the impact of the characteristic sample dimension is immediately recognizable. The length scale parameter l is involved in all subsequent computations with the value stated in Table 3. Figure 4. Distribution of the effective GND density g eff i for the cantilever beam sample #1 along the middle beam axis (from left to right) for the particular case H e 0 " 1 GPa. The corresponding contour plot within the central X 1´X3 -plane is additionally shown for the 8FI8FI-element formulation (cf. Figure 7 for the 20RI8FI-element formulation). Both formulations yield very similar results, indicating that the GND density field is not affected by the choice of element formulation.
The results in terms of stress-strain curves for all three cantilever beam samples are presented in Figure 6. A clearly recognizable increase in strength is obtained for all three beam sizes compared to the reference saturation stress, i.e., the size-independent response measured in [12]. The smallest sample (cantilever beam #1, thickness of t " 2.5 µm) exhibits the stiffest response. With increasing beam thickness, the bending response becomes softer. Accordingly, the response of cantilever beam #3 (thickness of t " 5 µm) is the softest and for t " 3.5 µm (cantilever beam #2), the response is located in between. The obtained size effect is well captured by the non-local crystal plasticity model. The slop in the elasticity-dominated regime shows the exact opposite trend. Here, the slope is determined by the second moment of area and, hence, by the beam thickness t. At a later stage of deformation, the back-stress effect resulting from the storage of GNDs becomes dominant in terms of the work-hardening rate.
By considering the distribution of GNDs within the central cross section of the beams, i.e., within the X 1´X3 -plane as displayed in Figure 7, it is seen that high densities of GNDs are accumulated at the supported end of the beams where the deformation is concentrated. In this region, GNDs pile up along the neutral plane (zero stress isoline), which does not necessarily coincide with the middle beam line. In fact, a shift of the stress field towards the bottom is caused to some extent by the supported end of the cantilever beam.
For all three sample sizes, a similar distribution of the effective GND density is found with a maximum value of g eff i « 7.0ˆ10 13 m´2, located at the lower half of the clamped side. In other words, the population of GNDs directly scales with the deformation gradient imposed by the normalized deflection as this is the same for all three beam geometries. Nevertheless, the resulting back-stress contributes differently for the samples leading to higher bending stresses σ ben for thinner beam samples. This can be partially explained by the fact that the computed bending stress. σ ben " M p b {S p depends on the plastic section modulus S p " wt 2 {4 which is a pure geometrical information. It is assumed that plastic deformation dominates which holds true at least in the most relevant strain regime ε ben ą «0.35, i.e., in the regime where the bending stress is saturating (cf. [12]). Further, M p b " Fl b denotes the plastic bending moment. Thus, the resulting flow stress scales inversely proportional to the square of the beam thickness t. Accordingly, there must be a scaling effect resulting from the beam width w, even if this effect is expected to be smaller. As w is kept constant in the present study, a pure dependence on the sample thickness is obtained. In the experiments, however, an influence of the beam width w is likely due to dimensional deviations caused by the fabrication process. Figure 5. Distribution of the effective plastic slip γ eff for the cantilever beam sample #1 along the indicated paths of the central cross section. In addition, the distribution is shown within the central X 1´X3 -plane. Results show some noticeable differences between both element formulations, in particular, close to the supported end, which refer mainly to the magnitude of γ eff .
The role of GNDs is further assessed by means of the statistically stored dislocation density. To do so, we introduce the effective SSD density ζ eff i as a function of γ eff and the effective free path of moving dislocations L eff where i is the total (effective) dislocation density, and K is a material constant. Equation (34) represents a modified version of the originally proposed relation by Essmann and Mughrabi [40,41]. Corresponding interaction processes are intrinsically considered via the calculation of γ α and g iα , respectively. K " 10 ( [42]) and ζ eff 0 " 2.0ˆ10 12 m´2 ( [13]) are given in the literature. Figure 6. Bending response of cantilever beam samples with different thicknesses. The strengthening effect associated with the accumulation of GNDs is the strongest for cantilever beam #1, having the smallest thickness of t " 2.5 µm. The softest response is obtained for cantilever beam #3 with t " 5 µm. The general trend in terms of smaller is stronger is captured well by the underlying model. All computations refer to the 20RI8FI-element formulation.
In Figure 8, the contour plots illustrate that the density of SSDs is naturally concentrated within highly deformed zones. More interestingly, the magnitude of ζ eff i reduces considerably with decreasing thickness of the cantilever beam. This indicates that the impact of GNDs on the evolution of SSDs is increased for smaller samples. In other words, the magnitude of SSDs approaches the magnitude of GNDs with decreasing sample dimensions, leading to a more pronounced influence on the bending response in case of smaller beams coming from accumulated GNDs. This indicates that the location of most prominent dislocation accumulation, for instance in terms of the total dislocation density ρ eff , shifts from the beam surface towards the beam center with decreasing beam thickness.

Bending Size Effect-Impact of a Saturating GND Density
The plastic deformation is strongly localized at the supported end of the beam. In this region, the number of GNDs increases with increasing bending load due to compatibility reasons. Accordingly, the hardening contribution from the higher-order gradient Div i pκ α q continuously increases which is reflected by the slope in the hardening behavior, cf. for instance Figure 6. However, from a physical point of view, the number of dislocations stored locally cannot be arbitrary high [43]. For that reason, the impact of a maximum permissible GND density on the bending response is investigated for cantilever beam #1 by setting up two different saturation values which are applied to each slip system independently: g max iα " 1ˆ10 13 m´2 and g max iα " 2ˆ10 13 m´2, respectively. Hence, the evolution equations for edge and screw GND densities are exposed to the following case differentiations: $ & % Equationp14q, for g e iαpn`1q ă g max iα ; 0, for g e iαpn`1q ě g max iα , and 9 g s iα " $ & % Equationp15q, for g s iαpn`1q ă g max iα ; 0, for g s iαpn`1q ě g max iα .
(37) Figure 7. Distribution of the effective GND density g eff i for all three sample sizes within the central cross section (X 1´X3 -plane) after 10% normalized deflection. GNDs pile up along the neutral plane indicated by the dashed line. Although the distribution of g eff i appears to be similar with respect to all three beam sizes, their impact on the bending response is strongly correlated to the beam thickness via the plastic section modulus S p . In view of the effective SSD density shown in Figure 8, a strong effect on the evolution is found due to the impact of GNDs. Figure 8. Distribution of the effective SSD density ζ eff i for all three sample sizes within the central cross section (X 1´X3 -plane) after 10% normalized deflection. The density of SSDs approaches the density of GNDs with reducing beam thickness, leading to a more pronounced influence on the mechanical bending response coming from accumulated GNDs.
As in the previous cases, all computations are based on the 20RI8FI-element formulation and are performed up to 10% normalized deflection. The resulting stress-strain curves are presented in Figure 9 together with the size-independent reference and the response of an unrestricted GND density evolution. Those cases represent the lower resp. upper bound of the mechanical bending behavior in terms of strength. The bending response in case of a limited GND evolution shows a saturation-like hardening behavior where the increase in strength relative to the size-independent case is related to the applied GND density limit. Hence, the stress saturation level is steered by the magnitude of g max iα in the way that a higher saturation limit causes a delayed deviation from the unrestricted case. In fact, once the saturation limit of g max iα is reached, the size-dependent micro-hardening stresses become decoupled from the gradient of plastic slip. Accordingly the size-dependent hardening contribution coming from Div i pκ α q is limited.
In related experiments, the bending response for a cantilever beam of size 2.5ˆ5.0ˆ16.3 µm (tˆwˆl b ) was found to show a stress saturation already at about 3.5% normalized deflection [12]. Thus, a saturation limit for GNDs around 1ˆ10 13 m´2 appears to be reasonable with respect to the current problem (cf. Figure 9). In comparison to the size-independent case, this yields an increase of about 118 MPa in flow stress which is solely related to the geometrically necessary storage of dislocations. Furthermore, it can be seen from Figure 9 that the bending size effect is conserved when comparing the bending response between the differently sized cantilever beams for the particular case if g max iα " 1ˆ10 13 m´2. The resulting effective back stress τ eff b , defined as is compared in Figure 10 for selected magnitudes of g max iα and sample thicknesses t. As can be seen, the size-dependent hardening contribution is much higher for the higher GND density limit which is consistent with the determined bending response in Figure 9. At the same time, the impact of the thickness is negligible for a particular g max iα value. Figure 9. Impact of a maximum permissible GND density g max iα on the mechanical bending response of Cu single crystals. As seen for cantilever beam # 1, the resulting saturation level is steered by the magnitude of g max iα . Furthermore, the bending size effect is conserved when comparing the response between the differently sized cantilever beams for a GND saturation limit of g max iα " 1ˆ10 13 m´2.

Discussion
The bending size effect is characterized by an increasing strength with decreasing sample thickness. A strong correlation between flow stress and beam thickness was found by Motz et al. [12] for thicknesses in the range 7.5 µm to 1 µm. A similar trend was found by Demir et al. [26] for alternative single crystal geometries with average thicknesses in the range 4.23 µm to 1.02 µm. In this size regime, the bending size effect is associated with a combination of different mechanism. Besides the geometrically necessary storage of dislocations, dislocation starvation and dislocation source limitation are known to affect the plastic deformation behavior. Dislocation source limitation plays a dominant role for very small beam sizes as the statistical distribution of dislocation sources becomes then more and more crucial, cf. [44]. As a consequence of a limited availability of source density within the localized region of plastic deformation (supported beam end), the yield limit of the material may increase significantly. In this respect, the storage of GNDs imposes an additional resistance to dislocation nucleation, yielding an increasing nucleation strength. Last but not least, plasticity is strongly controlled by the initial dislocation density of the crystal [8]. In the particular case of microbending of single crystals, dislocations are able to leave the crystal through the free surfaces at some point [12]. This process of starving or escaping dislocations finally leads to the extinction of initially available dislocations [4]. Hence, with a sufficiently small number of obstacles in the single crystalline sample, the size-dependent flow stress is additionally governed by the applied deformation rate relative to the dislocation nucleation rate which determines the required stress level for continuing operation of individual dislocation sources, see also Balint et al. [45]. Keeping in mind that an accurate fabrication of micro-cantilever samples is a challenging task, there is currently no sufficient data available allowing for a meaningful interpretation of the impact of each of the sample dimensions-thickness, width, and length-independently. For that reason, strengthening effects mentioned above are neglected in the current numerical study. Instead, the focus is solely on the micromechanical role of GNDs and their impact on the mechanical bending response as a function of the beam thickness t.

Figure 10.
Effective back stress τ eff b for selected g max iα values and cantilever beam sizes. The size-dependent resistance to bending deformation increases with increasing GND density limit. The sample thickness has a negligible impact on the back-stress evolution for a fixed g max iα .
The numerical results correctly capture the commonly observed trend 'smaller is stronger'. In fact, the size-dependent strengthening effect due to the back-stress effect induced by the storage of GNDs is predicted very well for two different cases: (i) unrestricted evolution of GND densities and (ii) physically limited GND densities. In case (i), GND densities evolve with increasing plastic slip gradients which, for the particular case of cantilever beam bending, continuously grow within the deformation-localized region. In contrast, the evolution of GND densities is limited in case (ii) and completely vanishes if a certain saturation limit is reached, i.e., if the plastic strain gradient becomes very large. This mimics a more realistic micromechanical behavior as is supported by the characteristics of related experimental force-deflection curves. Moreover, the local number of GNDs to be stored locally cannot be arbitrary high (cf. [43]). The numerically determined flow stresses for the here investigated beam thicknesses and for the particular case of g max iα " 1ˆ10 13 m´2 are illustrated in Figure 11 along with available experimental data from the literature. The simulation data is solely a function of the beam thickness t while experimentally determined data sets might (undesired) be affected by the reduction of the other sample dimensions (l b and w) or by varying dimension ratios (l b {t and l b {w). Apparently, changing sample dimensions and their ratios add substantially to the complexity of the superimposed effect associated with the mechanisms discussed above. This might be one explanation for the deviations between the experimental data sets of [12,26]. A meaningful discussion regarding the impact of the initially available dislocation density is not possible due to the lack of data. Figure 11. Experimentally and numerically determined relation between flow stress σ ben and beam sample thickness t in microbending of Cu single crystals together with associated power fits. The deviations between the two experimental data sets might be explained by the different sample geometries: cantilever beams with rectangular cross section and high l b {t as well as l b {w ratios ( [12]); cantilever beams with trapezoidal cross section, high l b {t avg ratios, but low l b {w ratios ( [26]). The numerically determined flow stress values show a reasonable strengthening effect in the regime where GNDs dominate the micromechanical behavior of the crystal, i.e., in the range t Á 3 µm. For t À 3 µm, the impact of dislocation starvation and source limitation become crucial, leading to an even more pronounced increase in flow stress with decreasing beam thickness t. The data of Demir et al. refers to a flow stress measurement at 0.06 strain.
Both experimental data sets were fitted by a power function of the form at´b`c using the Levenberg-Marquardt algorithm, see Figure 11. Then, an error analysis is conducted for the data set of Motz et al. [12] in order to draw final conclusions regarding the statistical representativeness of the measurements. The standard deviation is calculated for each fitting parameter a, b, and c. The maximal and minimal standard deviation of each fitting parameter is considered independently to compute upper and lower confidence bounds while the other two parameters are kept constant in each case. The resulting confidence intervals are embedded in Figure 11 and are interpreted independently. The overall fit sensitivity associated with the prefactor a is small due to its narrow confidence interval over the entire beam thickness range. With respect to the exponent b, it is found that the power function is governed by the exponent b in the regime t ă 2 µm. As here, one would expect rather large experimental scatter, the power fit appears less sensitive compared to the regime t ą 2 µm. In fact, the confidence interval associated with the exponent b vanishes if t approaches «1 µm. Hence, we conclude that the determined exponent describes the experimentally measured relation between σ ben and t fairly well. The offset parameter c indicates a similar tendency as the exponent, i.e., for t ą 2 µm the power function becomes more sensitive with respect to c. In comparison to the simulation data which reflects a pure dependence on the beam thickness t, a good agreement is found for t Á 3 µm as the predicted data lies very close to the actual fit. Consequently, it is concluded that the mechanism associated with the geometrically necessary storage of dislocations governs the mechanical bending behavior in this range. For t À 3 µm, the impact of dislocation starvation and source limitation becomes obviously non-negligible. These findings fit also to DDD predictions of Hussein et al. [46] where single crystals with D ď 1.0 µm were found to be almost free of dislocations due to the limiting size whereas crystals with D ě 5.0 µm show pronounced dislocation activities.

Conclusions
Based on an extended gradient crystal plasticity model, the role of GNDs in the mechanical bending response of micron-sized, single crystalline Cu was investigated. The underlying model contains non-standard evolution relations for the edge and screw components of the slip system-based dislocation density vector, crystal-specific interaction relations, and higher-order gradient boundary conditions. The cantilever beam geometries considered in the numerical study allowed the examination of the strengthening effect associated with the geometrically necessary storage of dislocations solely as a function of the beam thickness which is a non-trivial task from an experimental point of view. Other size-dependent mechanisms such as dislocation starvation and source limitation were disregarded at the current stage. In particular, the influence of the beam thickness as well as the impact of a maximum permissible GND density was of primary importance. Besides, in relation to the coupled field problem between displacement and GND density degrees of freedom, a quantitative comparison between two different finite element formulations has been carried out. On the basis of our findings, we conclude: • The bending dominated deformation is captured more accurately by the mixed FE-formulation denoted as 20RI8FI. In contrast, the commonly applied linear FE-formulation (8FI8FI) overestimates the bending response for the size-independent as well as the size-dependent case.
The locking phenomenon only influences the predicted bending behavior (and not the predicted GND density) in the case of the 8FI8FI-element formulation. • The bending size effect is captured by the theory to the extent caused by geometrically necessary storage of dislocations. This size-dependent strengthening effect can be explained as follows: (i) Similar dislocation pile-ups have been found around the neutral plane where dislocations get stuck rapidly and lose the ability to accommodate the beam bending, independent of the beam size. The impact of the resulting back-stress effect on the bending response is nevertheless higher for the smallest beam as the bending stress is inversely proportional to the square of the beam thickness. The same holds for the flow stress computation in related cantilever beam bending experiments; (ii) In contrast to the distribution of the GND density, a much higher population of SSDs was found for the largest cantilever beam sample which indicates that the bending behavior is here mainly governed by random trapping processes. However, with decreasing beam thickness, these processes become less pronounced. This is supported by the fact that the magnitude of the SSD density becomes comparable to the one of the GND density in the case of the thinnest beam sample. Consequently, the impact of GNDs on the mechanical bending response is most pronounced in the thinnest beam sample. Accordingly, the location of maximum dislocation storage was found to shift from the sample surface towards the beam center when decreasing the beam thickness. • A physically motivated limitation of the GND density was incorporated into the model by modified evolution equations for the edge and screw GND density components. In the current crystal plasticity framework, this was done at the nodal level as GND densities were treated as additional degrees of freedom. This leads to a bending response with saturation-like hardening behavior-which is in accordance with experimental findings. At the same time, the smaller is stronger trend was conserved in accordance to the unrestricted case. In the end, a saturation limit of «1ˆ10 13 m´2 was found to match well the characteristics in the bending response of related experimental data where a flow stress saturation was obtained at about 3.5% normalized deflection. • Numerically determined flow stresses using a saturation limit of «1ˆ10 13 m´2 show a reasonable strengthening effect in the beam thickness range t Á 3 µm. The predicted flow stress of cantilever beam #3 is in great accordance with experimental data. The flow stress associated with cantilever beam #2 still shows an acceptable accuracy as it lies within the confidence interval of the related experimental data. For the thinnest beam sample, a considerable contribution from another size-dependent mechanism occurs. Based on this, it can be argued that GNDs dominate the micromechanical bending response in the thickness range t Á 3 while other mechanisms such as dislocation starvation and source limitation become crucial for t À 3 µm where an even more pronounced increase in flow stress is experimentally measured.