1 Introduction

Over the past years, computational solid mechanics has become an integral part of theoretical materials science. Significant attention has been focused on mesoscale continuum mechanics where material properties and structural dimensions are important. Such formulations are intermediate between a direct atomistic simulation and an unstructured continuum description of the deformation processes. A variety of theoretical frameworks are emerging to describe inelastic deformation at the mesoscale: in this study we shall focus on one of these methods viz. discrete dislocation plasticity (DDP). In discrete dislocation plasticity, the dislocations are treated as line singularities in an elastic solid. A many-body interaction problem involving the discrete dislocations needs to be solved together with a complimentary more conventional elasticity boundary value problem.

Effect of structural size on the nominal material strength has been observed in many experimental investigations and backed by theoretical predictions. Sources of plastic size effects in crystalline materials include:

  1. (1)

    Geometrically necessary dislocations (GNDs). This has been the focus of phenomenological higher order plasticity theories [1, 2] and is used to explain size effects in bending, indentation, etc.

  2. (2)

    Constrained dislocation glide by internal interfaces such as grain boundaries. This is typically the main reason for size effects in phenomena such as the Hall-Petch effect [3, 4].

  3. (3)

    Source limited plasticity. Continuum plasticity assumes that at any point where a flow condition is met, plastic deformation will take place. Nevertheless, this does not need to be the case at very small sizes such as frictional asperity contact [5].

  4. (4)

    Dislocation starvation. The size effects in micro-pillar compression have been primarily associated with this effect [6,7,8] where dislocations exit the specimen at a rate faster than they are nucleated, and thus there is no build-up of dislocations within the specimen.

Following the pioneering work of Van der Giessen and Needleman [9], the DDP method has been shown to successfully predict numerous observations of plasticity size effects due to the above mentioned phenomena at the micron and sub-micron length scale. These include size effects in composites [10], bending [11], indentation [12], uniaxial compression [13], and under constrained shear [14]. The framework has also been used to investigate crack growth under monotonic [15] and fatigue loading [16]. Numerical methods to extend the framework to three-dimensional (3D) problems [17] and quasi-3D or the so-called 2.5D [18] have also been developed. In their study, Yasin et al. [19] and Vattré et al. [20] have presented a framework coupling continuum elasticity with 3D discrete dislocation dynamics. This 3D framework has been used to investigate a range of problems including the analysis of micro-cracks behaviour in high-cycle fatigue conditions [21] and the behaviour of micro-pillars under uniaxial compression [22]. The 3D framework has also been applied in more complex settings such as the modelling of polycrystalline materials [23], indentation [24], radiation hardening [25], and fatigue crack propagation [26]. All these studies have been conducted in a small strain context wherein finite deformation effects such as lattice rotations and changes in the geometry of the body are not accounted for.

Bending of specimens has been extensively used to investigate size effects in crystalline plasticity since the early work of Stölken and Evans [27]. These experiments attempted to impose a pure bending stress field on the specimens, but did not directly measure the moment-curvature relation. More recently, small scale bending experiments were conducted by first using focussed ion beam (FIB) to cut cantilever beams from single crystals and then applying tip loads using nano-indentation instruments [28, 29]. These experiments enable the direct measurement of load versus displacement relations and are used to infer plasticity size effects. All bending experiments involve significant lattice rotations and changes in the specimen geometry. Discrete dislocation plasticity investigations of such experiments [11, 30, 31] have neglected these finite strain effects and usually attributed size effects in bending to GNDs.

The primary question addressed in this study is whether lattice rotations and the associated curving of the slip planes affect the influence of GNDs in governing plasticity size effects in bending. To address this question we shall employ the recently developed finite strain DDP formulation of Irani et al. [32]. The outline of the paper is as follows: First, the small strain and finite strain DDP formulations are briefly described. Next, the finite strain and small strain DDP predictions of the cantilever bending problem are presented for geometrically self-similar beams. These predictions are then used to rationalise size effects in cantilever bending experiments.

2 Formulation of the problem

The bending/tensile response of single crystals is investigated here using both small and finite strain discrete dislocation plasticity in a two-dimensional (2D) setting. The crystals are taken to be elastically isotropic with Young’s modulus E and Poisson’s ratio \(\nu \). Plane strain conditions are assumed with deformations restricted to the \(X_1-X_2\) plane. Here, we briefly summarise both the small and finite strain formulations. For more details see the references cited.

2.1 Small strain discrete dislocation plasticity

The small strain DDP framework is well-established (see for example Ref. [11]) and here, we only summarise some of the salient points. In the DDP formulation plastic deformation and flow, is described by the nucleation and glide of discrete edge dislocations. Consider a 2D body under plane strain conditions such that at time t, the body contains N edge dislocations, with the dislocations being represented as line singularities in an elastic medium, with Burgers vector b. In this method, in order to obtain the solution, the problem is decomposed in two additive parts and the field quantities are computed using superposition. First, the singular \(({\tilde{\;}})\) field associated with the N dislocations are calculated analytically. Typically, solutions of dislocations fields in an infinite medium are used to represent the \(({\tilde{\;}})\) field, but it is equally possible to use other solutions such as those for dislocations in a half-space. Next, the complete solution is obtained by adding an image \(({\hat{\;}})\) field that ensures the boundary conditions are satisfied. Thus, the displacements, strains, and stresses are expressed as

$$\begin{aligned} u_i={\hat{u}}_i+{\tilde{u}}_i,\qquad \varepsilon _{ij}=\hat{\varepsilon }_{ij}+{\tilde{\varepsilon }}_{ij},\qquad \sigma _{ij}={\hat{\sigma }}_{ij}+{\tilde{\sigma }}_{ij}, \end{aligned}$$
(1)

respectively, where the \(({\tilde{\;}})\) field is the sum of the fields of the individual dislocations (labelled by \(1 \le I \le N\)) in their current positions, i.e.

$$\begin{aligned} {\tilde{u}}_i\equiv \sum _{I=1}^N {\tilde{u}}_i^{(I)},\qquad {\tilde{\varepsilon }}_{ij}\equiv \sum _{I=1}^N\tilde{\varepsilon }_{ij}^{(I)},\qquad {\tilde{\sigma }}_{ij}\equiv \sum _{I=1}^N{\tilde{\sigma }}_{ij}^{(I)}. \end{aligned}$$
(2)

The image \(({\hat{\;}})\) field in the other hand is obtained by solving a small strain elasticity boundary value problem [9].

In this method, at the beginning of a calculation the crystal is stress- and dislocation-free. The long-range interactions of the dislocations are accounted for through their elastic fields and for the short-range interactions constitutive rules are prescribed. New dislocation pairs are generated by simulating Frank-Read sources. In two dimensions, this is mimicked by discrete point sources randomly distributed on discrete slip planes. These sources generate a dislocation dipole with their Burgers vectors aligned with the slip plane direction. This occurs when the magnitude of the Peach-Koehler force \(f^{(I)}\) on source I exceeds a critical value \(\tau _\mathrm{nuc}b\) during a time period \(t_\mathrm{nuc}\). The sign of the dipole is determined by the sign of the resolved shear stress along the slip plane. Furthermore, the distance between the two nucleated dislocations, \(L_\mathrm{nuc}\), is taken such that the attractive stress field, which the dislocations exert on each other is equilibrated by a shear stress of magnitude \(\tau _\mathrm{nuc}\). Annihilation of two opposite signed dislocations on a slip plane occurs when they are within a material-dependent critical annihilation distance, \(L_e\). The magnitude of the glide velocity \(v^{(I)}\) along the slip direction of dislocation I is taken to be linearly related to the Peach-Koehler force \(f^{(I)}\) through the drag relation

$$\begin{aligned} v^{(I)}=\frac{1}{B}\,f^{(I)}, \end{aligned}$$
(3)

where B is the drag coefficient. Obstacles to dislocation motion are modelled as points associated with a slip plane. The dislocations gliding on the slip plane of an obstacle, get pinned as they try to pass through this obstacle. Moreover, obstacles release the pinned dislocations when the Peach-Koehler force on the obstacle exceeds \(\tau _\mathrm{obs}b\).

2.2 Finite strain discrete dislocation plasticity

The small strain DDP framework neglects lattice rotations, as well as the changes in the geometry of specimens due to deformation. These effects can be significant in the bending of relatively small specimens. Thus, we also investigate the bending of single crystals using the finite strain discrete dislocation plasticity framework [32, 33]. Analogous to the small strain DDP method, the finite strain DDP framework assumes that (1) lattice strains remain small away from the dislocation cores, and (2) the elastic properties are unaffected by slip. However, in contrast to the small strain calculations, the finite strain framework accounts for: (1) finite deformation-induced lattice rotations, and (2) the effect of shape changes due to slip on the momentum balance. Compared to the small strain formulation, the finite strain DDP formulation is less established, and thus we present this method more elaborately; nevertheless, readers are referred to Ref. [32] for further details.

Fig. 1
figure 1

Deformation of the crystal envisioned as slip in which the crystal lattice is unaffected, followed by the rotation and stretching of the material. For simplicity of representation, only one slip system is illustrated

Given the nonlinearity of the finite strain formulation, the total displacement rate is written as the superposition of the analytically known \(({\tilde{\;}})\) field of dislocations and the \(({\hat{\;}})\) field that enforces the boundary conditions, i.e.

$$\begin{aligned} {\dot{u}}_i(X_j,t)={\dot{{\hat{u}}}}_i(X_j,t)+{\dot{{\tilde{u}}}}_i(X_j,t), \end{aligned}$$
(4)

where \(X_j\) denotes the position of a material point in the initial configuration. The material \({\dot{u}}_{i,j}^\mathrm{m}\) and lattice \({\dot{u}}^\mathrm{e}_{i,j}\) velocity gradients are then calculated as

$$\begin{aligned} {\dot{u}}_{i,j}^\mathrm{m} ={\dot{{\hat{u}}}}_{i,j}+{\dot{{\tilde{u}}}}^\mathrm{d}_{i,j}, \end{aligned}$$
(5)

and

$$\begin{aligned} {\dot{u}}^\mathrm{e}_{i,j} ={\dot{{\hat{u}}}}_{i,j}+{\dot{{\tilde{u}}}}_{i,j}, \end{aligned}$$
(6)

respectively, where \((\ )_{,i}\) denotes the spatial gradients with respect to the initial configuration, \(\partial (\ )/\partial X_i\). The gradient \({\dot{{\tilde{u}}}}^\mathrm{d}_{i,j}\) is calculated by numerically differentiating the velocity field \({\dot{{\tilde{u}}}}_i\) with respect to \(X_j\) using the finite element shape functions in the initial configuration and thereby including the contributions from the deformations due to slip. On the other hand, \({{\dot{{\tilde{u}}}}}_{i,j}\) is calculated by analytically differentiating \({\dot{{\tilde{u}}}}_i\) with respect to \(X_j\) and thus does not include slip contributions. Hence, \({{\dot{{\tilde{u}}}}}_{i,j}\) are the lattice velocity gradients. Analogously, the Cauchy stress field \(\sigma _{ij}\) is written as

$$\begin{aligned} \sigma _{ij}={\hat{\sigma }}_{ij}+{\tilde{\sigma }}_{ij}, \end{aligned}$$
(7)

where \({\tilde{\sigma }}_{ij}\) is the sum of the stresses due to all dislocations in the body analogous to Eq. (2) while the stress field \({\hat{\sigma }}_{ij}\) corrects for the boundary conditions. It then follows that the weak form of the equilibrium equation over the current domain \(\varOmega ^*\) is

$$\begin{aligned} \int _{\varOmega ^*}{\hat{\sigma }}_{ij}\,{\hat{\eta }}_{i;j}\,\mathrm{d}\varOmega =\int _{S^{*}_T}(T_i^*-{\tilde{T}}_i^*)\,{\hat{\eta }}_i\,\mathrm{d}S, \end{aligned}$$
(8)

where \(T_i^*=\sigma _{ij}n^*_j\) are the tractions specified on the boundary \(S^*_T\) with outward normal \(n^*_j\) in the current configuration and \({\tilde{T}}_i={\tilde{\sigma }}_{ij}n^*_j\). In addition, the displacements are specified on the boundary \(S^*_U\) and enter the solution through boundary conditions. Here, \((\ )_{;i}\equiv \partial (\ )/\partial x_i\), i.e. spatial gradients with respect to the current configuration and \({\hat{\eta }}_i\) are a set of continuous and arbitrary trial displacement fields. With the deformation gradient given by \(F_{ij}=\delta _{ij}+u_{i,j}^\mathrm{m}\), the 2nd Piola-Kirchhoff stress \(S_{ij}\) is related to \(\sigma _{ij}\) via

$$\begin{aligned} S_{ij}=F^{-1}_{im}\, \sigma _{mn}\,F^{-1}_{jn}\,J, \end{aligned}$$
(9)

where \(\delta _{ij}\) is the Kronecker delta and \(J\equiv \mathrm{det}(F_{ij})\). The superposition (7) implies that

$$\begin{aligned} {\hat{S}}_{ij}=F^{-1}_{im}\,{\hat{\sigma }}_{mn}\,F^{-1}_{jn}\,J, \end{aligned}$$
(10a)

and

$$\begin{aligned} {\tilde{S}}_{ij}=F^{-1}_{im}\,{\tilde{\sigma }}_{mn}\,F^{-1}_{jn}\,J, \end{aligned}$$
(10b)

such that \(S_{ij}={\hat{S}}_{ij}+{\tilde{S}}_{ij}\). Using Eqs. (5) and (10a), the weak form of the equilibrium statement (8) is written in the undeformed configuration as

$$\begin{aligned}&\int _{\varOmega }\left( {\hat{S}}_{ij}\,{\hat{\eta }}_{i,j}+{\hat{S}}_{ij}\,{\hat{u}}_{k,i}\,{\hat{\eta }}_{k,j}\right) \mathrm{d}\varOmega \nonumber \\&\quad = \int _{S^*_T}(T_i^*-{\tilde{T}}_i^*)\,{\hat{\eta }}_i\,\mathrm{d}S - \int _{\varOmega }\left( {\hat{S}}_{ij}\,{\tilde{u}}^\mathrm{d}_{k,i}\, {\hat{\eta }}_{k,j}\right) \mathrm{d}\varOmega , \end{aligned}$$
(11)

where \(\varOmega \) denotes the domain in the undeformed configuration. This weak form of the equilibrium statement with traction boundary conditions on \(S^*_T\) and displacement boundary conditions on \(S_U^*\) is amenable to solution for the \(({\hat{\;}})\) field by a conventional finite element technique. Unlike the small strain formulation, this complimentary problem for the \(({\hat{\;}})\) field is nonlinear and needs to be solved iteratively as explained in Ref. [32].

With lattice strain assumed to be small, the stresses are assumed to be related to strains via the linear elastic Hooke’s law. However, in the finite strain setting this requires special care which is described here in some detail. Here, analogous to conventional crystal plasticity, the total deformation gradient \(F_{ij}=F_{im}^eF_{mj}^p\) is decomposed into two sequential operations as shown schematically in Fig. 1. First, the material slips on the currently active slip planes. These slips do not affect the orientation or the structure of the crystal lattice. The deformation at this stage is described by a deformation gradient \(F^p_{ij}\). Second, the lattice and material deform together and the elastic response of the lattice and any rigid body rotations are realised via the deformation gradient \(F^e_{ij}\equiv \delta _{ij}+u^{e}_{i,j} \). The constitutive relation for lattice elasticity is, therefore, written in terms of a material stress \(\Sigma _{ij}\) in the intermediate configuration of Fig. 1 as

$$\begin{aligned} \Sigma _{ij}=F^{e\,-1}_{im}\,\sigma _{mn}\,F^{e\,-1}_{jn}\,J^e=\hat{\Sigma }_{ij}+\tilde{\Sigma }_{ij}, \end{aligned}$$
(12a)

where \(J^e\equiv \mathrm{det}(F^e_{ij})\) and Eq. (7) then gives

$$\begin{aligned} \hat{\Sigma }_{ij}=F^{e\,-1}_{im}\,{\hat{\sigma }}_{mn}\,F^{e\,-1}_{jn}\,J^e, \;\;\;\mathrm{and}\;\;\; \tilde{\Sigma }_{ij}=F^{e\,-1}_{im}\,{\tilde{\sigma }}_{mn}\,F^{e\,-1}_{jn}\,J^e. \end{aligned}$$
(12b)

With lattice strain assumed to be small, the stress \(\Sigma _{ij}\) in the intermediate configuration is assumed to be related to the lattice Green-Lagrange strain \(E_{ij}\) via the linear elastic Hooke’s law

$$\begin{aligned} \Sigma _{ij}=L_{ijkl}E_{kl}, \end{aligned}$$
(13a)

with

$$\begin{aligned} E_{ij}=\frac{1}{2}\left( F^e_{mi}F^e_{mj}-\delta _{ij}\right) . \end{aligned}$$
(13b)

Combining Eqs. (12) and (13) it follows that

$$\begin{aligned} {\hat{\Sigma }}_{ij}=L_{ijkl}E_{kl}-F^{e\,-1}_{im}\,{\tilde{\sigma }}_{mn}\,F^{e\,-1}_{jn}\,J^e. \end{aligned}$$
(14)

Subsequently, by using Eqs. (10a), (12b), (14) as well as the fact that \(F_{ij}=F_{im}^eF_{mj}^p\), the 2nd Piola-Kirchhoff stress \({\hat{S}}_{ij}\) follows as

$$\begin{aligned} {\hat{S}}_{ij}=F^{p\,-1}_{im}\left( L_{mnqp}E_{qp}-F^{e\,-1}_{mp}\, {\tilde{\sigma }}_{pq}\,F^{e\,-1}_{nq}\,J^e\right) F_{jn}^{p\,-1}\,J^p, \end{aligned}$$
(15)

where \(J^p\equiv \mathrm{det}(F^p_{ij})\).

2.2.1 Finite strain discrete dislocation plasticity constitutive rules

We now summarise the plane strain short-range constitutive rules, highlighting the differences between the small strain and finite strain formulations. In the finite strain model the dislocations are no longer confined to a fixed slip plane due to both finite lattice rotations and slip on intersecting slip systems. Thus, the basic entity is a slip system (i.e. the orientation of the slip plane normal and the slip direction) rather than a slip plane. Moreover, because of lattice rotations, the orientation of a nucleated dislocation dipole (the two dimensional analogue of a nucleated loop) varies with position. The local lattice rotation is given by

$$\begin{aligned} \varphi = \sin ^{-1}(R_{21}), \end{aligned}$$
(16)

where the rotation tensor \(R_{ij}\) is written in terms of the right stretch tensor \(U_{ij} \) as \(R_{ij} = F_{ik}^{e}\,U_{kj}^{-1}\). Thus, the orientation of a dislocation on a slip system oriented at \(\phi ^{(\alpha )}\) in the undeformed configuration is \(\phi ^{(\alpha )}+\varphi \) in the current configuration.

At time t there are N dislocations in the body and the Peach-Koehler force \(f^{(I)}\) on dislocation I is given by

$$\begin{aligned} f^{(I)}=\left( {\hat{\sigma }}_{ij} + \sum _{J \ne I} {\tilde{\sigma }} ^{(J)}_{ij} \right) \, b^{*(I)}_j m^{*(\alpha )}_i, \end{aligned}$$
(17)

where \({\tilde{\sigma }}^{(J)}_{ij}\) is the stress field of dislocation J at the position of dislocation I while \(m^{*(\alpha )}\) is the unit vector normal to slip plane \(\alpha \) and \(b^{*(I)}_j\) the Burgers vector of dislocation I in the current configuration. The Peach-Koehler force includes the long-range interactions with all other dislocations in the material. This force determines the evolution of the dislocation structure, accounting for glide, generation, annihilation, pinning at and releasing from the obstacles. Similar to the small strain case, the magnitude of the glide velocity along the current slip direction of dislocation I is taken to be linearly related to the Peach-Koehler force such that the velocity \(v_i^{(I)}\) of dislocation I is given as

$$\begin{aligned} v_i^{(I)}= \frac{1}{B}\,f^{(I)}s^{*(\alpha )}_i, \end{aligned}$$
(18)

where B is the drag coefficient and \(s_i^{*(\alpha )}\) is the unit vector along the slip direction of slip system \(\alpha \) in the current configuration. Here, it is assumed that the drag coefficient B is constant throughout the body. We also do not account for any changes in the resistance to dislocation motion near a free surface associated with the energy required to create a new free surface when the dislocation exits.

New dislocation pairs are generated by simulating Frank-Read sources. In two dimensions, this is mimicked by discrete point sources on a slip system which generate a dislocation dipole with their Burgers vectors aligned with \(s_i^{*(\alpha )}\). This occurs when the magnitude of the Peach-Koehler force at that source exceeds a critical value \(\tau _\text {nuc}b\) during a time period \(t_\text {nuc}\). The sign of the dipole is determined by the sign of the resolved shear stress along the current slip direction. The distance between the two dislocations at nucleation, \(L_\text {nuc}\), is taken such that the attractive stress field that the dislocations exert on each other is equilibrated by a shear stress of magnitude \(\tau _\text {nuc}\). Annihilation of two opposite signed dislocations on slip system \(\alpha \) occurs when they are sufficiently close together. This is modelled by eliminating the two dislocations when they are within a material-dependent critical annihilation distance \(L_{e}\). Unlike in the small strain formulation where only opposite signed dislocations on a given slip plane can annihilate each other, in the finite strain context opposite signed dislocations on a given slip system can annihilate each other. Thus, annihilation of two opposite signed dislocations on a particular slip system occurs when they are within a distance equal to \(L_{e}\), irrespective of their current slip planes. Obstacles to dislocation motion are modelled as points associated with a slip system. Dislocations on the slip system of an obstacle, get pinned as they try to pass through that point. Again, unlike the small strain case, dislocations and obstacles are associated with a slip system rather than a slip plane. Thus, dislocations that glide on the slip system of an obstacle and are within a distance equal to \(L_{e}\) from it, get pinned to that obstacle. Pinned dislocations can only pass through an obstacle when their Peach-Koehler force exceeds an obstacle dependent value \(\tau _\text {obs}b\).

2.3 Specification of the cantilever beam bending boundary value problem

The cantilever beam analysed is sketched in Fig. 2 and comprises a “half dog-bone” type geometry similar to the specimens in experiments in which the cantilever is cut out from a large single crystal [28]. The length of the entire specimen is \(L_R+a\) with a gauge section of length a, width W, and a root of width \(W_R\). We fix the origin of the coordinate system at the cantilever root such that \(X_1\)-axis coincides with the axis of symmetry of the specimen and the loading is applied at the cantilever tip at \((X_1,X_2)=(L_R+a,W/2)\). The boundary conditions imposed to simulate cantilever bending are

$$\begin{aligned} {\dot{u}}_2 = -{\dot{\delta }},\;{\dot{T}}_1=0, \quad&\text {on} \quad \; (X_1,X_2) = (L_R+a,W/2), \end{aligned}$$
(19a)

along with traction-free boundary conditions on all surfaces other than the cantilever root where a built-in boundary condition is specified as

$$\begin{aligned} {\dot{u}}_1 = {\dot{u}}_2 = 0,\; \quad&\text {on}\quad \; X_1 = 0. \end{aligned}$$
(19b)

The traction-free boundary conditions also imply that dislocations can exit the specimen from all edges other than the edge at \(X_1=0\) where the zero-displacement constraints are specified. The work conjugate force P to the applied displacement imposes a bending moment \(M\equiv Pa\) at the root of the gauge section of the cantilever beam and a loading rate \(|{\dot{\delta }}|/a = 500 /\mathrm{s}\) is employed.

Fig. 2
figure 2

Sketch of the boundary value problem of the tip loading of the cantilever beam. The sketch shows the crystal in the double slip configuration and includes the coordinate system employed and labels of the critical dimensions

Fig. 3
figure 3

Sketch of uniaxial tension/compression of the dog-bone specimen. The coordinate system employed and critical dimensions are included in this figure where the crystal is shown in the double slip configuration

Two specimen geometries are analysed in this study.

  1. (1)

    Short beams. These beams have dimensions with ratios \(a/W=3\), \(W_R/W=3\), and \(L_R/W=5.6\). The taper angle as defined in Fig. 2 is then \(\alpha =20^\mathrm{\circ }\). The three geometrically self-similar specimens that are analysed are specified by the widths \(W=0.5\), 1.0, \(1.5\;\upmu \hbox {m}\).

  2. (2)

    Long beams. The only difference in this case is that the gauge section aspect ratio is increased to \(a/W=6\) with the beam root dimension ratios kept fixed at the same values, i.e. \(W_R/W=3\) and \(L_R/W=5.6\), such that again \(\alpha =20^\mathrm{\circ }\). For this case we also analyse three geometrically self-similar specimens with \(W=0.5\), 1.0, \(1.5\;\upmu \hbox {m}\).

We use a conventional finite element method to solve the complementary problem for the \(({\hat{\;}})\) fields. Thus, interior material points do not become boundary points and traction free boundary conditions and the displacements are prescribed on the same material points throughout the deformation history. In all calculations a uniform finite element grid with constant strain triangles with minimum side length \(e=0.05\;\upmu \hbox {m}\) is used and time step of \({\varDelta {t}} = 0.5\) ns is employed to resolve the dislocation dynamics. It is worth emphasising that the finite strain calculations are considerably limited by the distortion of the finite element mesh as the deformations are typically highly localised to a limited number of slip planes. Thus, large global deformations are not realisable in these calculations. In addition, calculations are conducted for three realisations of obstacle and source distributions and the reported results are averages over these three realisations.

2.4 Reference properties

The following set of material properties is used in both the small strain and finite strain DDP simulations. The crystal is assumed to be elastically isotropic with Young’s modulus \(E = 70\) GPa and Poisson’s ratio \(\nu = 0.33\), which are representative values for aluminium. In the undeformed configuration, the slip planes are spaced 100b apart, where \(b = 0.25\) nm is the magnitude of the Burgers vector. Simulations are reported for two crystallographic configurations of the crystal: (1) in the crystal oriented for single slip there is a single slip system oriented at \(\phi _1=30^{\circ }\) with respect to the \(X_1\)-axis in the undeformed configuration and (2) in the crystal orientated for symmetric double slip, there are two slip systems at \(\phi _1=30^{\circ }\) and \(\phi _2=-30^{\circ }\) with respect to the \(X_1\)-axis in the undeformed configuration.

In all calculations, the specimen is initially dislocation-free, but dislocations are generated from discrete sources placed on the slip planes within the crystal. The sources are taken to have a Gaussian strength distribution with a mean source strength \({\bar{\tau }}_\mathrm{nuc}=50\) MPa and a standard deviation of 10 MPa with a nucleation time \(t_\mathrm{nuc}=10\) ns. These Frank-Read sources are randomly distributed on the slip planes with a density \(\rho _\text {src} \approx 30 \; \upmu \text {m}^{-2}\). Similarly, obstacles are also distributed over the slip planes with a density \(\rho _\text {obs} \approx 60 \; \upmu \text {m}^{-2}\). These obstacles pin dislocations as long as the shear stress on the obstacle is below the obstacle strength \(\tau _\text {obs}=150\) MPa. The drag coefficient for glide is \(B=10^{4}~\text {Pa}\cdot \text {s}\), a representative value for several FCC crystals [34] and the critical distance for annihilation is \(L_e = 6b\).

Fig. 4
figure 4

Finite strain DDP predictions of the variation of the uniaxial tension and compression response of the \(W=0.5\) and \(1.0\;\upmu \hbox {m}\) crystals in the a single slip and b double slip configuration. The results are shown in terms of the magnitude of applied nominal stress \(\sigma _\mathrm{nom}\) versus the nominal strain \(\varepsilon _\mathrm{nom}\). Here, solid and dashed lines represent the results obtained for uniaxial tension and compression, respectively

Fig. 5
figure 5

Small strain DDP predictions of the variation of the uniaxial tension and compression response of the \(W=0.5\) and \(1.0\;\upmu \hbox {m}\) crystals in the a single slip and b double slip configuration. The results are shown in terms of the magnitude of applied nominal stress \(\sigma _\mathrm{nom}\) versus the nominal strain \(\varepsilon _\mathrm{nom}\)

3 Comparison of finite strain and small strain DDP predictions

Pure bending of crystals using discrete dislocation plasticity has been investigated by a number of authors initiated by the pioneering work of Cleveringa et al. [11]. However, most small-scale experiments are unable to generate pure bending loading with experiments typically carried out in a cantilever beam setting [28]. In cantilever beam bending, stress gradients exist along both the beam axis and width. Moreover, in this setting the imposed loads also generate shear stresses. In their study, Tarleton et al. [31] performed small strain DDP calculations of cantilever beam bending and argued that the GND structure formed in cantilever beam bending differed from that observed in pure bending. Finite strain effects such as lattice rotations can have a significant influence on the GND structures in bending and in this section, we aim to quantify these effects.

3.1 Uniaxial tensile/compressive response of the crystal

Before proceeding to discuss the bending response of the specimens, it is instructive to set a baseline by describing the response of the crystals under uniaxial tension and compression. To realise this aim, dog-bone specimens with a geometry very similar to those used in the cantilever beam bending simulations are employed. This specimen geometry is sketched in Fig. 3 and is identical to the “short beam” specimens expect for the fact that a complete dog-bone is used. In these samples the \(W_R\) wide grip sections exist at both ends of the gauge section with dimensions \(a\times W\), such that the total length of the specimen is \(L_t=2L_R+a\). Using the coordinate system sketched in Fig. 3, tension/compression loading is imposed by applying boundary conditions

$$\begin{aligned} {\dot{u}}_1&= {\dot{U}},\;{\dot{T}}_2=0, \quad \quad \text {on} \quad \; X_1 = L_t, \end{aligned}$$
(20a)
$$\begin{aligned} {\dot{u}}_1&= -{\dot{U}},\;{\dot{T}}_2=0, \quad \text {on} \quad \; X_1 = 0, \end{aligned}$$
(20b)

in addition to traction-free boundary conditions on all other surfaces of the specimen. Rigid body motion is constrained by specifying displacement \(u_2=0\) on one additional material point in the specimen. A loading rate \(|{\dot{U}}|/a=500/\mathrm{s}\) is employed and we report results for the tension and compression of two geometrically self-similar specimens with \(W=0.5\) and \(1.0\;\upmu \hbox {m}\). The results are reported in terms of the applied nominal stress \(\sigma _\mathrm{nom}\) calculated as

$$\begin{aligned} \sigma _\text {nom}= \frac{1}{W}\left| \int _{S_R} T_1 \mathrm{d}s\,\right| . \end{aligned}$$
(21)

Here, the integration is performed along the boundary \(S_R\) that lies along \(X_1 =L_t\) as a function of a measure of nominal strain defined as \(\varepsilon _\mathrm{nom}\equiv 2|U|/a\). Calculations are conducted for three realisations of obstacle and source distributions and the tension/compression curves presented are averages over these three realisations.

The finite strain DDP tension and compression responses of the \(W=0.5\) and \(1.0\;\upmu \hbox {m}\) crystals oriented for single and double slip are plotted in Fig. 4a, b, respectively. The corresponding small strain DDP predictions are included in Fig. 5. The following four key features emerge from these results:

  1. (1)

    The small strain DDP results are nearly identical to the finite strain predictions, except for the single slip \(W=0.5\;\upmu \hbox {m}\) specimen where the small strain results show a larger hardening regime.

  2. (2)

    There is negligible tension/compression asymmetry due to any finite strain effect for both the single and double slip cases.

  3. (3)

    The crystals oriented for single slip display an initial elastic response followed by a mildly hardening plastic behaviour. By contrast, in the double slip configuration the crystals display nearly no plastic hardening.

  4. (4)

    The crystals oriented for single slip have a small plastic size effect with the strength of the \(W=1.0\;\upmu \hbox {m}\) specimens slightly lower compared to the \(W=0.5\;\upmu \hbox {m}\) specimens. No discernible size effect is observed for the crystals oriented in double slip.

The absence of a clear specimen size effect on the tension/compression response is in contrast to the results previously reported by Deshpande et al. [13] and Balint et al. [35]. This difference is rationalised as follows. The tension/compression size effect is primarily due to the so-called dislocation starvation effect that requires dislocations to exit the specimens faster than the nucleation rate of dislocations. Thus, dislocation starvation is only present at a low obstacle density so that newly nucleated dislocations are not held-up at obstacles and can exit the specimens from the free surfaces rapidly. Here, a rather high obstacle density equal to twice the source density is employed, which implies that even the smallest specimens are not in the dislocation starvation regime.

Fig. 6
figure 6

Sketch of the \(W=1.0\;\upmu \hbox {m}\) tensile specimen with the window over which the dislocation distributions are shown in b and c. Finite strain DDP predictions of the dislocation and \(\sigma _{11}\) distributions in the \(W=1.0\;\upmu \hbox {m}\) specimen at \(\varepsilon _\mathrm{nom}=0.008\) for the b single slip and c double slip configuration. These distributions are shown over the central section illustrated in a

The high obstacle density also results in the mildly hardening response seen for the crystals oriented for single slip. Recall that the same source and obstacle densities are used in both the single and double slip configurations. Thus, the number of obstacles per slip plane is larger in the single slip configuration case. The consequence of this is seen in Fig. 6. In Fig. 6b, c, we include the dislocation distributions in the single and double slip cases, respectively, at an applied tensile strain of \(\varepsilon _\mathrm{nom}=0.008\) for the \(W=1.0\;\upmu \hbox {m}\) specimens. The corresponding distributions of the stress \(\sigma _{11}\) is also presented in these figures. These distributions, as shown in Fig. 6a, are plotted over a central section of length \(L_0 = 4a/3 = 4\;\upmu \hbox {m}\). In the single slip configuration, the high obstacle density employed blocks the nucleated dislocations. This results in dislocation pile-ups as seen in Fig. 6b while no such pile-ups occur in the crystals oriented for double slip as observed in Fig. 6c. Eventually, these pile-ups lead to the mildly hardening responses observed in the stress-strain curves of the single slip case.

Fig. 7
figure 7

Comparison between finite and small strain DDP predictions of the evolution of the bending stress \(\sigma _b\) with the measure of rotation \(\theta \) for crystals in a single slip and b double slip configurations. Results are shown for three selected specimen sizes W of the short beams. Here, solid and dashed lines represent the results obtained by small strain and finite strain DDP calculations, respectively

3.2 Cantilever beam bending

The uniaxial tension/compression responses of the specimens presented above suggest that the crystals in both their configurations are approximately elastic-perfectly plastic and display negligible size dependence under uniaxial loading over the range of specimen sizes considered here. In conventional continuum plasticity the uniaxial response is sufficient to characterise the bending response of the crystal. Employing this assumption, we expect the plastic collapse moment \(M_{f}\) for the cantilever beams of Fig. 2 to be given by

$$\begin{aligned} M_{{f}}=\frac{\sigma _{f} W^2}{4}, \end{aligned}$$
(22)

where \(\sigma _{f}\) is the plastic flow strength of the crystal in uniaxial tension/compression as deduced from Fig. 4. In DDP simulations, the plastic flow strengths under bending and uniaxial loading are not necessarily the same. Thus, in order to obtain a measure of the differences between the conventional continuum plasticity predictions and those obtained from DDP, we define a bending stress measure

$$\begin{aligned} \sigma _b\equiv \frac{4M}{W^2}, \end{aligned}$$
(23)

where \(M\equiv Pa\) is the imposed bending moment at the root of the gauge section of the specimen due to the applied tip load P. In the following, we present the results by describing the evolution of the bending stress \(\sigma _b\) as a function of a measure of beam rotation \(\theta \equiv |\delta |/a\). In this section, we will focus on understanding the impact of finite deformations on the bending response of the material, and this aim will be achieved by comparing the predictions of finite strain and small strain DDP formulations for \(\sigma _b\).

Predictions of \(\sigma _b\) versus \(\theta \) are presented in Fig. 7a, b for the single and double slip configurations, respectively, of the short beams. These figures show that after an initially elastic response, there is a knee in the \(\sigma _b\) versus \(\theta \) curve corresponding to the onset of dislocation activity and the initiation of plastic deformation. However, unlike the uniaxial loading case, the bending stress displays a hardening plastic response. Moreover, the bending stress in the plastic regime is size dependent with the smaller specimens exhibiting a larger bending stress. This is consistent with a wide body of both experimental [27, 28] and computational [11, 30, 31] investigations that suggest that the bending stress is size dependent in the micron size regime. We shall explore the origins of this size effect in more details in Sect. 4. Very briefly, the size effect is primarily a manifestation of the so-called GNDs. As will be seen subsequently, in these simulations the GNDs are primarily dislocation pile-ups that cause back-stresses and thereby plastic hardening. These dislocation pile-ups and their associated back-stresses lead to the size effects seen here. Moreover, the plastic hardening is more pronounced in the crystals in the single slip configuration (Fig. 7a) compared to the double slip case (Fig. 7b). In order to rationalise this observation we plot the dislocation structures in the \(W=0.5\) and \(1.5\;\upmu \hbox {m}\) specimens at \(\theta =0.045\) in Fig. 8a, b, respectively, for the crystals oriented in the single slip configuration. The corresponding plots for crystals oriented in the double slip configuration are included in Fig. 8c, d, respectively. In Fig. 8 contours of the lattice rotation \(\varphi \) are also included. Longer dislocation pile-ups are observed in the specimens with only one active slip system and we attribute the stronger plastic hardening to these longer pile-ups. This higher hardening rate is also consistent with the fact that plasticity is more constrained when only one active slip system is present compared to the double slip case and also that under uniaxial loading the single slip configuration exhibited a higher plastic hardening rate compared to the double slip case (Fig. 4).

In Fig. 7, predictions of the small strain DDP formulation are included along with the finite strain predictions. While the small strain DDP predictions are very similar to their finite strain counterparts, the small strain predictions are slightly stronger in all cases. These differences can be understood by noting the lattice rotation distributions plotted in Fig. 8. These lattice rotations tend to curve the dislocation slip planes as shown schematically in Fig. 9 where some selected slip planes are shown in both the undeformed and deformed configurations for the double slip case. The curving of the slip planes tends to relax the back-stresses created by dislocation pile-ups and hence reduce the plastic hardening. The effect of lattice rotations and curving of the dislocation slip planes is only included in the finite strain DDP formulation and hence the small strain predictions are slightly stronger compared to the corresponding finite strain predictions. However, as mentioned in Sect. 2.3, the finite strain computations presented here are limited by mesh distortion effects such that the lattice rotations obtained for the maximum applied cantilever tip displacements were less than \(3^\mathrm{\circ }\) (i.e. \(|\varphi |\le 0.05\)). Over this limited range of lattice rotations, the predicted finite strain effects are small, but we anticipate these effects to become significant if computations could be conducted for larger cantilever tip displacements. We emphasise here that the observed differences between the small strain and finite strain DDP formulations are primarily due to finite strain discrete dislocation plasticity effects rather than finite strain elastic effects: a comparison between finite strain and small strain elastic bending of the specimen is shown in “Appendix” section to emphasise this point.

Fig. 8
figure 8

Finite strain DDP predictions of the dislocation structure and lattice rotation \(\varphi \) in the short beam specimens at an applied \(\theta =0.045\) in the current configuration. Results are shown for the single slip case specimens of size a \(W=0.5\;\upmu \hbox {m}\) and b \(W=1.5\;\upmu \hbox {m}\) while the double slip case with specimen sizes \(W=0.5\) and \(1.5\;\upmu \hbox {m}\) are shown in c and d, respectively

4 Size effects in cantilever beam bending

Over the range of the cantilever tip displacements investigated here, finite strain effects were shown not to be significant. We thus proceed to investigate the underlying plasticity mechanisms for cantilever beam bending using the small strain DDP formulation. We first restrict our attention to the short beam specimens and then continue with exploring the effect of the beam aspect ratio a / W.

Predictions of \(\sigma _b\) versus \(\theta \) are included in Fig. 7a, b, respectively, for specimen sizes \(W=0.5\), 1.0, \(1.5\;\upmu \hbox {m}\). As was mentioned in Sect. 3, these results indicate that the bending stress is strongly hardening in the plastic regime, and the single slip configuration is stronger compared to the double slip case. Moreover, a plastic size effect exists with the bending stress in the plastic regime being larger for the smaller specimens.

In order to clarify these trends we define a flow bending stress \(\sigma _{f}\) as the value of \(\sigma _b\) at \(\theta =0.045\). Predictions of the variation of \(\sigma _{f}\) with W for the geometrically self-similar short beams is plotted in Fig. 10. The error bars in Fig. 10 indicate the scatter in the results based on the different realisations of the source and obstacle distributions computed here. The variation of the bending flow stress with W is reasonably well described by a relation of the form

$$\begin{aligned} \sigma _{f}=\sigma _0\left( \frac{W}{W_0}\right) ^{-n}, \end{aligned}$$
(24)

where \(\sigma _0\) is the flow stress of a specimen of width \(W=W_0\), i.e. \(\sigma _{f}\) increases with decreasing W. This curve is included in Fig. 10 and with \(W_0=1.0\;\upmu \hbox {m}\) and \(n\approx 0.3\), it is seen to fit the data with reasonable accuracy for both the slip configurations of the short beam specimens.

Fig. 9
figure 9

Sketch of the curving of the slip planes in the end loaded cantilever beam due to lattice rotations. The sketch shows the specimen in the double slip configuration in both the undeformed (dashed lines) and deformed (solid lines) states

Fig. 10
figure 10

Summary of the small strain DDP predictions of the variation of the bending flow strength \(\sigma _{f}\) with specimen size. Predictions are included for both the single and double slip configurations and the short and long beam specimens. The error bars indicate the variation in the results between the different realisations of the source and obstacle distributions computed here. The lines are fit to the data using Eq. (24). Here, solid and dashed lines represent the results obtained for short and long beams, respectively

Fig. 11
figure 11

Small strain DDP predictions of the dislocation structure and stress \(\sigma _{11}\) at an applied \(\theta =0.045\) in the single slip case. Results are shown for the short beam specimens of size a \(W=0.5\;\upmu \hbox {m}\) and b \(W=1.5\;\upmu \hbox {m}\) while the corresponding distributions in the long beam specimens of sizes \(W=0.5\) and \(1.5\;\upmu \hbox {m}\) are shown in c and d, respectively

The distributions of dislocations and contours of the stress component \(\sigma _{11}\) are plotted in Fig. 11 at an applied rotation \(\theta =0.045\). The results shown in Fig. 11a, b are for the \(W=0.5\) and \(1.5\;\upmu \hbox {m}\) short beam specimens, respectively, with the crystals oriented for single slip. These figures demonstrate that in the cantilever beams the dislocation activity is concentrated around the root of the gauge section of the specimen. The vast majority of these dislocations have the same sign as their oppositely signed counterparts have exited the specimen through the free surfaces. This implies that in these specimens there is a significant net Burgers vector due to the dislocations and hence a large GND density in Nye’s terminology [36]. These GNDs are dislocation pile-up structures and result in back-stresses on the dislocation sources counteracting their further operation. The size effect is largely associated with these GNDs as argued extensively in Ref. [1, 27].

4.1 Effect of cantilever aspect ratio

Consider the cantilever beam sketched in Fig. 2. The beam can either collapse by plastic bending with the formation of a plastic hinge at the root of the gauge section at \(X_1=L_R\) or by shear collapse of the gauge section. To illustrate this competition, consider the continuum isotropic plastic limit where the beam material has a uniaxial strength \(\sigma _Y\) and a shear strength \(\tau _Y\) for shearing in the \(X_1-X_2\) plane. Then the ratio of the shear collapse load \(P_\mathrm{s}\) to the bending collapse load \(P_\mathrm{b}\) is given as

$$\begin{aligned} \frac{P_\mathrm{s}}{P_\mathrm{b}}=4\frac{\tau _Y}{\sigma _Y}\left( \frac{a}{W}\right) , \end{aligned}$$
(25)

i.e. as the aspect ratio a / W increases, the shear collapse load increases with respect to the bending collapse load and hence, the deformation of the beam will become increasingly bending dominated. This basic notion holds irrespective of the details of the assumed plastic constitutive properties of the beam material. Thus, in order to investigate the contribution of the shear force on the plastic properties reported above, we proceed to investigate the bending response of the long beams with a gauge section aspect ratio \(a/W=6\) (as detailed in Sect. 2.3).

Predictions of \(\sigma _{f}\) for the long beam specimens are included in Fig. 10 along with the short beam results. It can be observed that the overall trends are very similar with \(\sigma _{f}\) decreasing with increasing W and the double slip case being slightly weaker compared to the single slip configuration. However, as it is shown in Fig. 10 the bending flow strength of the long beam specimens is appreciably smaller compared to the equivalent short beam samples. This is rationalised by noting that the long beams also have a lower dislocation density compared to the short beam case. We attribute this to the fact that for the same applied bending stress, the longer beams have lower shear forces (and stresses) acting upon them and this results in a lower density of statistically stored dislocations (SSDs). To clarify this, we plot in Fig. 11c, d the dislocation structure and the stress distributions in the \(W=0.5\) and \(1.5\;\upmu \hbox {m}\) long beam specimens, respectively, at \(\theta =0.045\). The dislocation structures in the long beam specimens are almost exclusively comprised of same-signed dislocations while there are relatively more opposite signed dislocations in the short beam specimens; see Fig. 11a, b. This confirms that the long beam specimens have a lower density of SSDs compared to the short beam case. Consequently, the lower dislocation density accounts for the lower bending flow strength in the long beam specimens, as dislocation motion is less inhibited by a dense dislocation structure.

The results for the long beam specimens demonstrate that the reduction in the SSD density does not change the trend of the variation of \(\sigma _{f}\) with W. Moreover, the bending flow stress is reasonably well described by Eq. (24) with \(W_0=1.0\;\upmu \hbox {m}\) and \(n\approx 0.3\). Thus, these results suggest that although the value of \(\sigma _0\) in Eq. (24) depends on the specimen aspect ratio, the value of \(\sigma _{f}/\sigma _0\) only depends on specimen width W.

5 Discussion

Initiated by the pioneering work of Stölken and Evans [27], numerous studies have confirmed that the bending strength of crystalline materials increases with decreasing specimen size; see for example Ref. [28]. This so-called size effect has traditionally been attributed to the fact that the GND density, which is associated with the imposed curvature in bending, increases with decreasing specimen size. The vast majority of models employed to capture this size effect are the so-called strain gradient plasticity models (see for example Ref. [37]), which explicitly include the contribution of GNDs to the strength and thereby capture the observed size effect. These continuum models thus have an in-built length scale as an additional parameter, which is calibrated either via experiments or lower length scale models. Another class of closely related models uses the concept of the Nye [36] dislocation density tensor to estimate the GND density, which enhances the physical basis of these higher order plasticity theories [38].

The drive to improve the mechanistic basis of the higher order plasticity theories has led to the development of what is now known as continuum dislocation plasticity (CDP) with the aim of describing the behaviour of the ensembles of dislocations via continuum field equations. For example, Sandfeld et al. [39] presented a numerical implementation of such a theory and applied it to the plane-strain micro-bending problem. Later, Le and Nguyen [40] (also see Refs. [41, 42]) proposed analytical descriptions for bending of single crystal beams within the CDP setting. In their work, in addition to the usual GND size effects, they also included the effect of specimen size dependent dislocation nucleation.

In the current study we employ discrete dislocation plasticity wherein the dislocation structures and the associated stresses are generated as a natural outcome of the boundary value problem. The main difference between the discrete dislocation plasticity and continuum models discussed above, is that in the discrete setting, GNDs and SSDs cannot be differentiated without including an artificial averaging length scale. Moreover, in the cantilever bending problem analysed here the lattice curvature varies over the length of the beam. Thus, it is not possible to quantify the GND density in a sensible manner from these discrete dislocation plasticity calculations. Therefore, making direct comparisons between the discrete dislocation plasticity and continuum plasticity descriptions is not instructive. Nevertheless, there is an excellent agreement between the discrete dislocation plasticity predictions and measurements of the size effects besides a good qualitative agreement with the above mentioned continuum plasticity predictions.

6 Concluding remarks

Cantilever beam configurations are commonly used to analyse the bending response of micro-sized specimens. Here, we have analysed the bending of tip loaded single crystal cantilever beam specimens with one and two active slip systems using DDP. Both finite strain and small strain DDP formulations are employed to investigate the influence of lattice rotations and changes in specimen geometry on the bending response of the specimens.

Over the range of specimen sizes analysed here, the bending stress versus applied tip displacement response has a strong plastic hardening component. This hardening rate increases with decreasing specimen size, and thus the bending flow strength is size dependent. The hardening rates are slightly lower when the finite strain DDP formulation is employed as curving of the slip planes is accounted for in the finite strain formulation. This curving relaxes the back-stresses in the dislocation pile-ups that result due to the bending stress field and thus reduces the hardening rate. However, over the range of applied tip displacements analysed here the small strain DDP formulation is shown to capture with excellent accuracy all the key features of the bending response of end loaded cantilever beams.

The small strain DDP formulation is then used to analyse size effects in these cantilever beam specimens. In line with the well-known pure bending case, the bending stress is shown to display a plastic size dependence. Moreover, the bending flow strength of beams with a larger aspect ratio is shown to be appreciably smaller compared to the equivalent short beam samples. We attribute this to the fact that for the same applied bending stress, the longer beams have lower shear forces (and stresses) acting upon them, and this results in a lower density of statistically stored dislocations: this reduction in the dislocation density lowers the bending strength.