A discrete slip plane model for simulating heterogeneous plastic deformation in single crystals

In small-scale mechanical tests, such as micropillar compression tests, plastic deformation is often localized in narrow slip traces. These slip traces result from a few dislocation sources with relatively low nucleation stresses that are present in the material. In order to accurately simulate such small-scale experiments, the stochastics of the underlying dislocation network must be taken into account, which is usually done by performing discrete dislocation dynamics simulations. However, their high computational cost generally restricts these simulations to small and simple geometries and small applied displacements. Furthermore, effects of geometrical changes are usually neglected in the small strain formulation adopted. In this study, a discrete slip plane model for simulating small-scale experiments on single crystals is proposed, which takes the most important characteristics of dislocation plasticity for geometries in the micrometer range into account, i.e. the stochastics and physics of dislocation sources. In the model, the properties of all lattice planes are sampled from a probability density function. This results in a heterogeneous ﬂow stress within a single crystal, unlike the uniform properties assumed in conventional crystal plasticity formulations. Moreover, the slip planes can be grouped together in bands via a weakest-link principle. The resulting equations are implemented in a standard crystal plasticity ﬁnite element model, using a ﬁnite deformation formulation. Within this setting, only the collective dislocation motion on glide planes is modeled, resulting in a signiﬁcantly lower computational cost compared to frameworks in which the dynamics of individual dislocations are considered. This allows for simulating multiple realizations in 3D, up to large deformations. A small case study on micropillar compression tests is presented to illustrate the capabilities of the model. (cid:1) 2021 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
Small-scale experiments on metal alloys have great scientific value because they can be used to test the mechanical response of isolated features in a microstructure, such as individual phases in multi-phase materials (Ghassemi-Armaki et al., 2014;Tian et al., 2020), or the effect of isolated grain boundaries by testing bi-crystals (Kirchlechner et al., 2017;Malyar et al., 2017).Such experiments have provided invaluable insights in the often complex and unexpected microdeformation mechanisms that control the macroscopic material behavior of the corresponding bulk polycrystalline materials from which the microspecimens were extracted.Moreover, the extracted local information, e.g. the critical resolved shear stress of individual slip systems, can be used to identify parameters in polycrystalline plasticity models used to improve advanced steel grades (Ghassemi-Armaki et al., 2013;Chen et al., 2014).One of the most prominent experimental techniques for testing crystals at the microscale is the micropillar compression test developed by Uchic, Dimiduk and coworkers (Uchic et al., 2004;Uchic and Dimiduk, 2005;Dimiduk et al., 2005).Since then, a large range of micropillar compression tests have been done on a variety of metals and alloys (Frick et al., 2008;Greer et al., 2005;Brinckmann et al., 2008;Zaiser et al., 2008).Besides compression, also tensile tests on single crystals have been performed (Kiener et al., 2008;Du et al., 2018).For more extensive overviews on small-scale testing the reader is referred to (Greer and de Hosson, 2011;Kraft et al., 2010;Uchic et al., 2009).
A prime reason for small-scale testing is the fact that the behavior of single crystals with dimensions in the (sub-) micron range is clearly different compared to that of the bulk material.Only a few localized slip traces are commonly observed.Furthermore, yield stresses of comparable samples can vary significantly.In larger samples, plasticity is dominated by the collective behavior of dislocations.At the scale of individual crystals, individual dislocation sources and/or obstacles play a more dominant role.As a result, plastic slip is concentrated on the weakest glide planes.Consequently, the stochastics of these weakest links may dominate the behavior of microcrystals (Norfleet et al., 2008;El-Awady et al., 2009;Parthasarathy et al., 2007;Pan et al., 2015).Other remarkable observations in microcrystals are size effects, despite the absence of strain gradients, and scale-free intermittent strain bursts (Miguel et al., 2001;Csikor et al., 2007).
In order to exploit the full potential of microscale experiments that allow for a detailed analysis of the complex dislocation mechanisms, each experiment should be complemented by large-strain simulations which account for the microstructure heterogeneities, stochastic behavior and geometrical effects due to shape changes.
Discrete dislocation dynamics (DDD) simulations have significantly contributed to the understanding of the physical processes of plasticity in microcrystals (Tang et al., 2008;Rao et al., 2008;Senger et al., 2008;El-Awady et al., 2009;Cui et al., 2014;El-Awady, 2015).In DDD, the dislocation network is modeled by discrete line segments and their evolution over time is simulated.However, because all dislocation segments need to be tracked individually, DDD simulations entail a significant computational cost.Therefore, simulations are mostly limited to small strains as well as small and simple geometries or are simplified in two dimensions.Otherwise, they typically require massive parallelization to obtain results in a reasonable time.To the best of the authors knowledge, the 2D framework developed by Irani and co-workers (Deshpande et al., 2003;Irani et al., 2015) is the only DDD framework formulated in a finite strain setting.This restricts the possibilities to use DDD for studying the influence of (i) geometric effects including large changes in specimen shapes, (ii) lattice rotations, (iii) influence of statistical variations and (iv) changes in the evolution of the local loading conditions arising in microscale experiments.
Besides DDD, continuum plasticity frameworks have been used for analyzing small-scale experiments on single crystals.For example, Zhang et al. (2006) used finite element simulations with isotropic plasticity to develop guidelines for the design of microcompression experiments.Similarly, Kiener et al. (2009) used isotropic plasticity simulations to study the influence of deformation constraints imposed at the top and bottom of micropillars.While isotropic plasticity simulations have clear merit, the discrepancy with small-scale experiments revealing slip on specific glide planes is obvious.Other groups have used crystal plasticity finite element (CPFE) simulations, in which the crystalline nature of the material is taken into account through the tensorial deformation resulting from plastic slip on pre-defined slip systems, to partially close this gap.For example, Raabe et al. (2007) used CPFE simulations to study the effects of the initial orientation, sample geometry and friction in micropillar compression tests.Similarly, Shade et al. (2009) studied lateral constraint effects, focussing on crystals orientated in single-slip and loaded in compression.Du et al. (2018) used CPFE to explain slip activity in microtensile tests of interstitial free ferrite.While CPFE does take the crystal orientation and associated dislocation glide in specific slip directions into account, it is still a continuum framework in which one generally assumes homogeneous properties throughout individual grains.Consequently, it does not include the discrete and stochastic nature of the underlying slip planes and dislocation structure, which is a key characteristic of small-scale tests on single crystals and which is adequately captured in computationally much more expensive DDD simulations.
Various modeling frameworks have been proposed that partially fill the gap between DDD and crystal plasticity theories.Zaiser and Moretti (2005) introduced a 2D model considering a single slip system which combined random stress fluctuations as a function of space and strain with long-range interaction stresses mediated by Green's function, resulting in plastic strain hetero-geneities.In the 1D gradient plasticity framework of Zhang and Aifantis (2011) a pillar was divided into elastic and plastic layers, where each plastic layer has its own yield stress.In this higherorder continuum theory a third order stress arises which accounts for significant changes in internal stresses as the plastic layers yield.By requiring this third order stress to be continuous across boundaries when two adjacent layers deform plastically, strain burst were predicted.Implementing their framework in a cellular automaton made it possible to also capture the stochastic variation between stress-strain curves (Konstantinidis et al., 2014).In the theoretical crystal plasticity model of Lin et al. (2015) the plastic strain is composed of a series of strain bursts that follow a power-law distribution, while the yield stress is based on the single-arm dislocation source model of Parthasarathy et al. (2007).However, none of these frameworks is able to simulate complex geometries together with plastic deformation resulting from slip on multiple slip systems.The dislocation-based CPFE framework introduced by Lin et al. (2016) was used to study the size-dependent deformation morphology of pillars.However, in that study localization into a single slip band had to be triggered by incorporating strain-softening into the slip resistance evolution law and by introducing a material imperfection in the center of the pillar.
The present study presents a model in which the collective behavior of dislocations gliding over discrete slip planes is modeled, without considering individual dislocations.Instead, it is assumed that the discrete slip planes exhibit a relative displacement (slip) which is smooth in time and space along the plane.Large differences in the amount of slip may however exist between neighboring planes.The properties of the slip planes are connected to the stochastics of the underlying dislocation sources.This results in a heterogeneous flow stress and a stochastic stress-strain response of the crystal.It is shown that the computational effort can be reduced significantly via a weakest-link principle that only considers the weakest glide plane in a slip band of finite width.With small adaptations, the resulting equations of this discrete slip plane model can be implemented in a standard CPFE model using a finite deformation formulation, enabling large strain simulations.The framework allows the simulation of small-scale experiments in the micrometer range, where the behavior is dominated by the weakest links and in which the plastic deformation is localized in a finite number of slip traces.
The structure of the paper is as follows.The model and corresponding assumptions will be introduced in Section 2. In Section 3 the assumptions are assessed, whereby the mesh independency of the model is demonstrated.Finally, in Section 4 a small case study is presented to illustrate its added value and its capabilities in terms of modeling small-scale experiments on single crystals.

Model assumptions
Dislocation glide is considered to be the main mechanism underlying plastic deformation.This mechanism results in crystallographic slip, i.e. the sliding of two crystal regions with respect to each other along a slip plane.Instead of modeling all dislocations individually, as is done in DDD simulations, slip is modeled as a continuous process in time and in space along a slip plane, similar to the perspective taken in crystal plasticity.However, unlike crystal plasticity, the plastic slip is concentrated in discrete, parallel slip planes and the amount of slip on these planes may vary significantly, even between neighboring planes.Initially, all discrete atomic slip planes of the possibly active slip systems are taken into account.The spacing of these slip planes is equal to where a is the lattice constant of the material and h; k and l are the Miller indices of the plane.The thin bands in between the discrete planes are considered to behave as a purely elastic continuum.An example is shown in Fig. 1.Here, the slip planes of two possible slip systems are depicted.Slip steps on both slip systems are indicated with lengths v 1 and v 2 , for which the resulting geometry is shown.

Model equations
The constitutive equation used to describe the plastic slip of adjacent atomic planes along the slip plane between them, as described in the previous section, is based on the widely adopted power-law relation for the plastic shearing rate introduced by Hutchinson and Hill (1976).Here, it is assumed that the amount of slip, v, on a certain slip system can be described by a viscoplastic (rate-dependent) relation, i.e. the velocity, _ v, between the crystallographic blocks separated by the slip plane of a slip system is dependent on the resolved shear stress, s, which is acting on the same slip system according to where s is the slip resistance, _ v 0 is a reference velocity and r is a rate sensitivity parameter.Note that when the rate sensitivity parameter approaches zero, a rate-independent response is recovered.
The initial value of the slip resistance, s 0 , contains three contributions, following the models by Parthasarathy et al. (2007) and Lin et al. (2015).These are a nucleation stress, s nuc , required to nucleate new dislocations, a lattice friction term, s fric , and a forest strength due to the initial dislocation density, q dis : where G and b are respectively the shear modulus of the material and the length of the Burgers vector.The three contributions in Eq.
(2) will be elaborated on in the next section (Section 2.3).Due to the evolution of the dislocation forest and dislocations encountering obstacles, the slip resistance usually increases during plastic deformation.This effect is accounted for by a phenomenological hardening law, as introduced by Bronkhorst et al. (1992) for crystal plasticity, given by where s 1 is the saturation value of the slip resistance, k 0 is the initial hardening rate and a is the hardening exponent.Note that this evolution is local to the slip plane, i.e. interaction effects such as latent hardening are neglected.

Dislocation source stochastics, introduced through the slip resistance
For the lattice friction and the initial dislocation density in Eq. ( 2), it is assumed that the spatial variation and variation between slip systems is negligible.Therefore, the same value is adopted for all slip planes in the crystal lattice.However, the nucleation stress can vary strongly between different slip planes and slip systems.It is relatively easy to nucleate dislocations on slip systems that contain dislocation sources, while nucleating dislocations on slip systems without a source requires a high stress since they have to be nucleated at the free surface (Xu et al., 2013;Hu et al., 2019).
In the proposed model, the nucleation stress is assumed to be constant over the length of a slip plane, but varies per atomistic slip plane.The variation in nucleation stresses is described by a probability density function, pðs nuc Þ.This distribution is composed of two underlying distributions: one distribution, g þ ðs nuc Þ, for slip planes that contain a dislocation source and another distribution, g À ðs nuc Þ, for slip planes without a dislocation source.The full distribution is obtained by employing the probability of finding a dislocation source on a slip plane, f src , resulting in where f src can be estimated by the dislocation source density q src , via with A the area of a lattice plane in the sample.Note that in these definitions of pðs nuc Þ and f src it is assumed that dislocation sources are distributed homogeneously over all lattice planes and that a single lattice plane will never contain more than one dislocation source.The latter assumption is reasonable when f src ( 1, which is the case for typical values of the dislocation source density q src . Different distributions can be used to describe the variation in dislocation source strengths.For example, (Shishvan and van der Giessen, 2010) used a log-normal distribution in their DDD simulations of copper thin films.El-Awady et al. (2009) used different Weibull distributions for the source length of Frank-Read (FR) sources in DDD simulations of nickel micropillars, while a more physics-based route was followed in the model of Parthasarathy et al. (2007), who proposed a distribution of the source length of single-arm (SA) sources.The latter two distributions are expressed in terms of source length.Therefore, they are rewritten here in Fig. 1.Schematic drawing of consecutive slip on two discrete slip planes.First, slip takes place along a slip plane of slip system 1, resulting in a slip step of length v1.Next, a slip step of length v2 takes place on a plane of slip system 2.This results in a deformed geometry with visible slip traces.
terms of source strength.The strength of a dislocation source can be estimated based on the line tension and is related to its characteristic length, k, by This relation can be used for both FR sources and SA sources.
Only the geometrical factor, a, varies for different types of sources (Pichaud et al., 1978;Rao et al., 2007).Two examples of nucleation stress distributions for lattice planes, indicated by A and B, are shown in Fig. 2. Here, the shapes of the left peaks are governed by underlying distributions g þ ðs nuc Þ, which thus reflect the stochastic nature of the nucleation stress of dislocation sources.For the distribution marked A, a log-normal distribution is used for g þ ðs nuc Þ, while a SA source distribution is adopted for the distribution denoted B. The right peaks in the figure are governed by the underlying distributions g À ðs nuc Þ, which are taken as narrow normal distributions at high stress around the theoretical material strengths for both distributions A and B. Since the behavior will be dominated by the weakest planes, the choice for g À ðs nuc Þ generally has a negligible effect.It only becomes important in the so-called dislocation starvation regime, which is only observed for pillars with diameters into the nanometer range (e.g.< 160 nm Shan et al., 2008), although the exact range of diameters in which dislocation starvation can be observed depends on the initial dislocation density (El-Awady, 2015) and the applied strain rate (Hu et al., 2019).The distribution denoted A in Fig. 2 will be explored in Section 3, while the distribution denoted by B will be explored in Section 4. Both source distributions will be discussed in more detail in these sections.

Model reduction
Modeling each individual discrete plane on which slip may occur in numerical simulations would require an excessively fine discretization, which would be computationally prohibitive.Therefore, a simple method to reduce the number of modeled planes is proposed.This method is based on the weakest-link principle.
First, instead of considering all slip planes of a slip system separately, planes are grouped together into bands of thickness l.The geometry is thus divided into multiple potential slip bands.This is done for every possible slip system independently, which means that slip bands of different slip systems may intersect.Consequently, every material point belongs to as many slip bands as there are slip systems.This is schematically shown in Fig. 3, where the slip planes of Fig. 1 are grouped together in a number of slip bands.Note that for a slip band width of l ¼ d the original atomistic model is recovered, i.e. each slip band then contains only a single slip plane.
Following the weakest-link principle, dislocation nucleation happens only in one or a few of the weakest dislocation sources.As a result, these weakest dislocation sources dominate the behavior of a single crystal.Therefore, it is assumed that within the group of slip planes within a slip band, slip only occurs on the weakest plane, i.e. the slip plane with the lowest nucleation stress.This single slip step is then distributed over the thickness of the slip band, which hence becomes elasto-plastic instead of elastic.
Following this concept, an amount of slip v on the weakest slip plane in a slip band is smeared out over a thickness l.This results in a shear strain in the slip direction with magnitude c ¼ v l .Eq. ( 1) may hence be rewritten as Similarly, Eq. ( 3) now reads As a further refinement of the modeling, cross hardening may be accounted for by introducing the hardening matrix q.This matrix has a value of 1 on the diagonal and a value of q n offdiagonally.The final result is given by where N is the total number of slip systems taken into account and a and b are used to denote a particular slip system.The form of Eqs.
(7) and ( 9), by construction, is identical to hardening laws used in many standard crystal plasticity models (Cuitiño and Ortiz, 1992;Kalidindi et al., 1992) based on the multiplicative split of the deformation gradient tensor into an elastic part, F e , and a plastic part, F p : and where the plastic velocity gradient tensor is given by wheres a 0 and ña 0 are respectively the slip direction and the slip plane normal of slip system a in the lattice configuration.The implementation of the model can thus take benefit of existing CPFE codes.Note that only the nucleation stress of the weakest slip plane in a slip band has to be known, since the whole slip band will experience the same nucleation stress.Therefore, it is no longer necessary to sample the nucleation stress of all slip planes within the band from the distribution defined in Eq. (4).Instead, order statistics is used to find the probability density function of the weakest plane in a slip band (David and Nagaraja, 2003): where Pðs nuc Þ is the cumulative distribution function of pðs nuc Þ and n is the number of slip planes grouped together in a single slip band.Sampling the nucleation stress for a slip band can thus be efficiently done from this weakest-link distribution, p min ðs nuc Þ.

Finite element implementation
As stated earlier, Eqs. ( 7) and ( 9) are common expressions used in rate-dependent CPFE models based on the multiplicative split of Fig. 2. Two different probability density functions for the nucleation stress, each having two peaks, one below 100 MPa, the other above 10 GPa.The shape of the left peaks, i.e. g þ ðsnucÞ, reflects the stochastic nature of the nucleation stress of dislocation sources, while the shape of the right peaks, i.e. g À ðsnucÞ, corresponds to the distribution of the nucleation stress required to nucleate dislocations without a dislocation source.Inset: zoom of the left peaks.the deformation gradient.Therefore, the implementation of the "lumped" discrete slip plane model is straightforward.Only the division of the geometry in different slip bands, as depicted in Fig. 3, requires special treatment.Since there are as many slip band orientations as there are slip systems, it is not possible to create a mesh that conforms to all these slip bands.Therefore, a regular mesh is used instead and for each integration point (not each element) it is determined, for each slip system, to which slip band it belongs based on its global position in the reference configuration.This means that in a single material integration point, the different slip systems will generally have different properties.Similarly, the properties at integration points within a single element will almost always be different.
A Saint Venant-Kirchhoff type of model is adopted for the elastic response of the material.The multiplicative split of Eq. ( 10) first maps the material from its initial configuration to a stress-free intermediate configuration through F p .The second Piola-Kirchhoff stress tensor defined in this intermediate configuration, S e , is computed by where 4 C is the fourth-order elasticity tensor and : denotes a double contraction.The resolved shear stress on a slip system in the intermediate configuration is now calculated as The finite element model is implemented in the form of a fortran user subroutine in the commercial finite element package Msc.Marc (Marc user's manual, 2014).The update of state variables F p and s is done implicitly via a trapezoidal integration scheme.

Model assessment
In this section the impact of grouping slip planes together in bands on the computed response is assessed.In particular, we study the influence of the thickness of the bands.This is done based on an idealized one-dimensional version of the model, with only one slip system, so that the comparison can be made with the reference case in which all planes are modeled individually, i.e. l ¼ d.At the end of the section, the finite element implementation is verified by a mesh convergence study.

Geometry and material properties
Throughout this section, we consider a problem which is inspired by the microtensile tests on a single crystal of ferrite as performed by Du et al. (2018).The dimensions of the gauge section of the specimen are 9 Â 3 Â 2 lm.The crystal considered here, denoted as Grain 2 in Du et al. (2018), is orientated with respect to the loading axis such that the highest Schmid factor has a value of 0.495, while the second highest Schmid factor takes a value of 0.445.
In this example, the adopted nucleation stress distribution corresponds to the distribution used by (Shishvan and van der Giessen, 2010); it is shown as distribution A in Fig. 2. The source distribution is a log-normal distribution in which 99.7% of the probability density function lies between the physically minimum and maximum values of the nucleation stress.The minimum value of the nucleation stress corresponds to the largest possible Frank-Read source with length k max that fits in the sample (equal to the smallest dimension of the geometry), i.e. s min nuc ¼ aGb kmax (see Eq. ( 6)).The maximum value of the nucleation stress corresponds to the theoretical strength of the crystal, s max nuc ¼ Gb 2pd (Hull and Bacon, 2011).With this choice for s min nuc and s max nuc the mean, l, and standard deviation, r, of the distribution are given by and The non-source distribution, g À ðs nuc Þ, is chosen as a normal distribution with mean s max nuc and a standard deviation of 0:01ðs max nuc À s min nuc Þ, such that it does not interfere with g þ ðs nuc Þ.However, as stated before, the precise choice of g À ðs nuc Þ has a negligible effect on the results.The resulting distribution, pðs nuc Þ, was already shown in Fig. 2 as Distribution A. All other relevant material properties can be found in Table 1.Note that the initial hardening rate, k 0 , as well as the saturation value of the slip resistance, s 1 , scale with the initial slip resistance, s 0 , and hence also vary per slip plane.

One-dimensional implementation
In order to explore some basic characteristics of the model, a one-dimensional (1D) implementation of the model with a simplified elastic response is used.This simplification allows us to model all atomistic slip planes of a slip system, i.e. the case l ¼ d, which would not be possible with a finite element implementation.Furthermore, given the low computational cost of the 1D model, many realizations can be simulated, so that mean values and standard deviations can be statistically assessed.
Only the slip system with the highest Schmid factor is considered, since it is expected that most slip activity will occur on this slip system.Furthermore, by neglecting geometrically nonlinear effects, it can be assumed that the elastic response of the whole system is given by where u p ðtÞ is the total prescribed displacement in the direction of the loading axis, v i is the slip displacement on the ith plane, h is the angle between the slip direction of the considered slip system and the loading axis, E is the Young's modulus and L is the length of the 1D domain.The resolved shear stress is related to the uniaxial stress via the Schmid factor, m, by Eq. ( 17), together with Eqs.(1), ( 3) and ( 18), upon sampling the distribution pðs nuc Þ as discussed above, can be solved with a standard solver for a set of nonlinear ordinary differential equations, e.g. the forward Euler method.

Slip band width
First, the effect of grouping atomistic slip planes together in a slip band of a certain width with properties governed by the weakest link will be examined.For computational efficiency of the 3D model it would be beneficial if the response is not too sensitive to the ratio of the band width to the atomistic slip plane distance, l=d, so that a large band width can be adopted.Therefore, the error made by grouping planes together in bands of width l is assessed here, using the 1D implementation of the model.
There are 23304 atomistic slip planes of the considered slip system that fit in the 9 lm long gauge section of the specimen.First, the nucleation strengths of all these slip planes are sampled from the distribution given by Eq. ( 4).For the realization obtained, the weakest planes are subsequently determined for a range of band widths, starting with a width of l ¼ 23304d (i.e. a single band) and then successively dividing it by two.This whole procedure is repeated for 100 different realizations.
The average responses, together with the standard deviations, obtained for the different band widths l are plotted in Fig. 4a; the average and standard deviation of the reference solutions in which every slip plane is tracked (i.e.l ¼ d) is also shown.This figure reveals that both the average response and the standard deviation quickly converge upon refining the band width.With a band width of l ¼ 2913d, i.e. the gauge section is divided into 8 slip bands, the maximum deviation in stress, relative to the response where all individual slip planes are taken into account (l ¼ d) is only 1.5 percent.Furthermore, the error introduced by weakest-link principle is much smaller than the standard deviation, even if all slip planes in the gauge section of the specimen are grouped together in just a single band.
Next, one of the 100 realizations with 23304 slip planes is examined in more detail.In Fig. 4b, the shear strain of this realization is plotted as a function of the longitudinal position for the different values of l at the same prescribed displacement u p ðtÞ.Consequently, the area below all curves is the same.The solid blue line (l ¼ d) represents the case where all individual atomistic slip planes are taken into account, where the most significant shear strain peaks are shown.Four extremely narrow peaks are visible, which indicates that almost all slip is localized on four atomistic slip planes.This implies that there are four relatively weak dislocation sources present on the primary slip system.The exact heights of these four peaks are not depicted in the figure, but are corresponding, from left to right respectively, to slip steps of v ¼ 0:43 lm, v ¼ 0:11 lm, v ¼ 0:11 lm and v ¼ 0:02 lm, i.e. most plastic deformation is concentrated in the most left peak.For l ¼ 23304d (brown dashed line), the strain is smeared out over the whole length of the geometry.In this case the response is fully governed by the weakest plane of the entire geometry, i.e. the model assumes that all plastic deformation is concentrated on the most left weak slip plane, but smears this slip across the entire length of the specimen.Refining once results in two strain plateaus, each smeared out over half the length of the geometry (purple dashed line).In this case the response is governed by the weakest plane in the left half and the weakest plane in the right half.More plastic deformation has occurred in the left half of the geometry, since this half contains a weaker plane than the right half of the geometry.For a slip band width of l ¼ 2913d (green dashed line) all four weak slip planes appear in different bands.This means that all four weak slip planes are taken into account by the model.For this specific case, this slip band width thus adequately captures the response of all atomistic slip planes.Refining the band width further (orange dashed line) results in more localized strain peaks.However, the number of strain peaks and the response of the system is no longer affected by this refinement.
In general, an error is made if more than one relatively weaker atomistic planes on which significant plastic deformation occurs fall within the same slip band.In such cases only the slip occurring on the weakest of these planes is taken into account and that on the other ones is neglected.The probability of such an error decreases with decreasing band width, which explains the convergence in Fig. 4a.The minimum band width that should be adopted is dependent on the adopted distributions and parameters.For example, the probability of finding a dislocation source, f src , has a significant effect on the band width which should be adopted because it determines the number of weak lattice planes present in the sample.

Mesh convergence
Based on the results presented in the previous section, a band width of l ¼ 2913d is adopted, i.e. 8 slip bands over the specimen length.Instead of using the 1D idealization, the model is now solved with the finite element method for a 3D specimen geometry of 9 Â 3 Â 2 lm.The geometry is discretized with hexahedral elements with 20 nodes and integrated with 8 Gaussian quadrature points (reduced integration), because this element is known to be less sensitive to volumetric locking compared to other standard continuum elements (Hughes, 1987).Initially, an element size of 0.5 lm is used, i.e. 8 Â 6 Â 4 cubic elements, which equals 0:85l.
Subsequently, the discretization is successively refined by a factor of 2 in all directions for the same realization of slip planes.We consider only one slip system to be able to compare with the 1D model.The considered realization is identical to the one that is used for Fig. 4b.The resulting stress-strain curves are presented in Fig. 5a.With the initial element size of 0:85l, the response is too stiff.The error is relatively large because the discretization is not conforming with the slip bands.Upon refinement, the response converges.Based on this analysis, an element size of h ¼ 0:212l, i.e. approximately 1=5th of the band width, seems appropriate.The result obtained with the idealized 1D model is also shown, by the solid purple line in the diagram.Both results adequately match.Mesh convergence has also been confirmed for simulations with all slip families included (not shown).
In Fig. 6 the contours of the accumulated slip are shown as computed with the coarsest mesh (h ¼ 0:85l) and the finest mesh (h ¼ 0:115l).Most slip is localized into a single band.Furthermore, the accumulated slip in that band has a much lower value for the coarse mesh than for the fine mesh.This is also visible in Fig. 5b, where the accumulated slip profile over the length of the geometry is plotted for different element sizes and as well as for the 1D idealization.This confirms that the element size needs to be sufficiently small compared to the band width and that an element size of h ¼ 0:25l is an appropriate choice.

Application to micropillar compression of nickel
As an application of the model we consider the microcompression tests on nickel micropillars, as done by Dimiduk et al. (2005).These micropillars are loaded in the ½269 direction of the FCC crystal.3D simulations up to large strains are performed.Because of the relatively low computational cost, multiple realizations, simulated with multiple sets of boundary conditions, can be considered.

Strength distribution and size effects
The measured yield stresses of pillars with different diameters show an apparent size effect, in the sense that the pillars with a smaller diameter result in a higher strength.The five slip systems with the highest Schmid factors are shown in Table 2.The adopted source nucleation stress distribution is based on the statistical model of single-arm (SA) sources by Parthasarathy et al. (2007).It has already been shown by these authors that such a distribution can adequately predict the size effect observed in micropillars.The distribution is based on the assumption that dislocation sources in such small circular pillars are dislocation lines which are pinned at one point, while the other end of the dislocation is at the free surface.SA sources thus only need one pinning point.These pinning points are assumed to be randomly distributed in the lattice planes.For cylindrically shaped specimens the lattice planes are ellipse-shaped due to their orientation with respect to the axis of the cylinder.The shortest distance from a pinning point to the edge of the ellipse is taken as the characteristic length of the SA source, k, which is related to the nucleation stress of the source through Eq. ( 6).This is schematically shown in Fig. 7a, where 3 arbitrary pinning points are depicted.
The probability density function of the nucleation stress of a SA source is now given by where R is the radius of the pillar and / is the angle between the slip plane normal and the pillar axis/loading direction.The nonsource nucleation stress is again taken as a normal distribution around the theoretical strength.The full distribution used is denoted by Distribution B in Fig. 2. The parameters of the model have been chosen such that the experimentally observed size effect, by Dimiduk et al. (2005), is also predicted with the model presented here.The full parameter set used may be found in Table 3. Parameters k 0 ; s 0 ; s 1 and a were used to fit the trend of the stress-strain response.f src mainly influences the scatter in the yield stress.All other parameters are physical quantities which are known or can be measured.
The 1D idealization of Section 3.2 is adopted, for which only the primary slip system is considered.Simulations have been done for 4 different pillar diameters.For every diameter 100 simulations have been performed and the stress at 0.2% strain is taken as engineering yield stress.The results are shown in Fig. 7b.Here, also a power-law is fitted on the simulation data.The obtained exponent of À0.68 is in good agreement with literature on micropillar compression tests, where an exponent between À0.6 and À0.9 is reported (Greer and de Hosson, 2011;Kraft et al., 2010;Uchic et al., 2009).Furthermore, the scatter in the simulation results is in reasonable agreement with the scatter observed in experiments.Standard deterministic plasticity models are only able to predict the average behavior of these experiments.The significant scatter suggest that individual samples can behave significantly different compared to this average behavior, which highlights the need to include stochastic effects in the model.

Slip activity
For the 3D finite element simulations, all 12 f111g < 110 > slip systems are taken into account.A rectangular pillar geometry with dimensions 5 Â 2 Â 2 lm is adopted.Because small deviations in boundary conditions may have a significant effect on such small pillars, two different sets of boundary conditions are applied.We consider the first case where the specimen is compressed but all other displacements and rotations remain free.This is done by prescribing the displacement in the loading direction in an average way.The full set of boundary conditions in these simulations is indicated as Case A in Fig. 8. Here, < u x > denotes the area average of the displacement in x-direction, i.e.

< u x
In Figs.9a and 10a the deformation at an applied strain of À5 % is shown for two different realizations.The colormaps show the maximum value of the shear strain out of all 12 slip systems.For the first realization (Fig. 9a), most strain is localized in one slip band of the ð111Þ½10 1 slip system.This is also the primary slip system, with a Schmid factor of 0.48.Some slip activity on the ð1 1 1Þ½101 slip system is also observed.This is the secondary slip system, with a Schmid factor of 0.40.This slip system is activated because its weakest plane has an initial slip resistance of 34 MPa, which is significantly lower than the second weakest plane of the primary slip system, which has a initial slip resistance of 62 MPa.This shows that Schmid's law in micropillars can already break down because of the microstructural stochastics of dislocation source strengths.The same point was demonstrated by Ng and Ngan (2008), who used an extension of the statistical model of Parthasarathy et al. (2007).One should also note that slip on this slip system would not be predicted with conventional CPFE, which again demonstrates the importance of including stochastic behavior at these small scales.For the second realization (Fig. 10a), all strain is localized in a single band of the primary slip system.
The same samples are also simulated with more constraining boundary conditions.The displacements at the top and bottom are now fully prescribed, as indicated by Case B in Fig. 8.This may be a more realistic set of boundary conditions as there is typically a significant amount of friction between the loading plates and the sample, since this improves the geometrical stability of micropillar compression tests (Zhang et al., 2006;Raabe et al., 2007).
The results of the simulations with the more constraining boundary conditions for the same two realizations considered before are shown in Figs.9b and 10b.The same two slip bands of the ð111Þ½10 1 slip system are again observed in both realizations,      together with the ð1 1 1Þ½101 slip system for the first realization.However, prominent slip activity in a band of the ð1 1 1Þ½110 slip system can also be seen for both realizations.This slip system has a Schmid factor of 0.23, which is only the fifth highest Schmid factor (Table 2).The fact that it is nevertheless activated can be explained by considering the projections of the slip directions on the y-z plane, which are shown in the polar plot of Fig. 11.Note that the slip directions are defined such that the x-components are negative, since these are the directions in which slip will occur in a compression test.Both the bottom and top of the micropillar are constrained in such a way that no displacement in this plane is allowed.Therefore, any displacement in y-or z-direction that would occur by slip on the primary slip system has to be compensated for -preferably by slip activity on some other system, or otherwise elastically.The slip direction of the primary slip system is ½10 1 (blue arrow).A slip system with ½110 as slip direction (purple arrow) is best suited to compensate the displacement of the primary slip system in the y-z plane because its projected vector on this plane points in almost the opposite direction.As a result, a band of the ð1 1 1Þ½110 slip system containing a relatively weak dislocation source will slip as soon as slip on a plane of the primary slip system occurs, despite its low Schmid factor of 0.23.This shows that it is important to analyze boundary conditions and to take geometrical effects into account when analyzing micropillar compression tests.Studying these effects would not be possible with 3D DDD simulations since the computational cost restricts these simulations to be done up to small strains only.Furthermore, geometrical changes and rotations are not taken into account due to the small strain formulation underlying most DDD frameworks.

Summary and conclusions
A model has been presented for the simulation of small-scale experiments on single crystals.The model is based on the sliding of discrete slip planes, instead of considering all dislocations individually.The properties of the slip planes are sampled from a probability distribution, which is based on the underlying stochastics of the dislocation sources.By making use of the weakest-link principle, lattice planes can be clustered together in slip bands.This allows the model to be implemented in a standard 3D CPFE framework accounting for large deformations.In this way, both crystal heterogeneity and geometrical effects can be taken into account.Moreover, it was shown that the weakest-link principle can be used with relatively large band widths.This considerably reduces the computational cost of the model.
Finally, the application of model was demonstrated on a compression tests of nickel micropillars.Here, it was shown that an appropriate distribution for the dislocation source strength allows one to capture the experimentally observed size effects.Additionally, slip is localized in only a few slip bands, which is characteristic for small-scale experiments on single crystals with dimensions up to a few micrometers.The added value of the model was demonstrated by analyzing the difference in slip system activity between two different sets of boundary conditions.Such simulations are relevant for the micropillar compression community, but cannot be done with standard CPFE which does not include stochastic influences of the underlying dislocation structure which have a dominant role at these small scales.
The (orders of magnitude) lower computational cost compared to DDD and the ability to take rotations and other finite deformation effects into account make the model attractive for analyzing small-scale experiments with dimensions in the order of micrometers and above.
A limitation of the framework in its current form is that it does not account for intermittent strain bursts, which are often attributed to the spontaneous shutdown of active dislocation sources and the activation of new dislocation sources (Rao et al., 2008;Tang et al., 2008;Cui et al., 2014).Instead, it is assumed that all dislocation sources are initially present in the material and the hardening behavior is smoothly captured.This aspect, and the extension to multiple grains are topics of further research.

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

Fig. 3 .
Fig.3.Division of the geometry into slip bands of width l for 2 different slip systems.The discrete slip steps, as depicted in Fig.1, are now distributed over these slip bands.

Fig. 4 .
Fig. 4. (a) The average stress-strain responses and standard deviations for different band widths obtained with the 1D idealization of the model.(b) The longitudinal strain profile over the length of the geometry for a single simulation obtained with the 1D idealization of the model.The strain profile corresponding to l ¼ d (blue solid line) consists of four very narrow peaks with significant amounts of slip (v ¼ 0:43 lm, v ¼ 0:11 lm, v ¼ 0:11 lm and v ¼ 0:02 lm) marking the locations where plastic slip has actually taken place.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.The slip contours obtained for the microtensile problem of Section 3 with element sizes h ¼ 0:85l (top) and h ¼ 0:115l (bottom).Note that the deformation is shown on true scale.

Fig. 7 .
Fig. 7. (A) A slip plane in a cylindrical micropillar of diameter D under an angle / with respect to the loading direction.Three pinning points and their shortest distance to the edge of the ellipse are depicted.(b) Yield stress as a function of pillar diameter.The simulation results are obtained with a 1D implementation of the model.Note that both axes are logarithmic.

Fig. 8 .
Fig. 8. Two different sets of boundary conditions, termed A and B, implemented in the finite element simulations of the compression tests of the nickel micropillars of 5 Â 2 Â 2 lm.

Fig. 9 .
Fig. 9. Contour bands of the maximum accumulated slip of all 24 slip systems for the two sets of boundary conditions defined in Fig. 8 (a): Case A and (b): Case B for realization 1.Note that the deformation is shown on true scale.

Fig. 10 .
Fig. 10.Contour bands of the maximum accumulated slip of all 24 slip systems for the two sets of boundary conditions (a): Case A and (b): Case B for realization 2. Note that the deformation is shown on true scale.

Table 1
Model parameters used for simulating the tensile tests of ferrite.

Table 2
The five slip systems with the highest Schmid factors of a FCC crystal loaded in the ½269 direction.

Table 3
Model parameters used for simulating compression tests of nickel.