A continuum model for dislocation dynamics in three dimensions using the dislocation density potential functions and its application in understanding the micro-pillar size effect
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J u l A continuum model for dislocation dynamics in three dimensionsusing the dislocation density potential functions and its applicationto micro-pillars
Yichao Zhu a , Yang Xiang a, ∗ a Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay,Kowloon, Hong Kong
Abstract
In this paper, we present a dislocation-density-based three-dimensional continuum model,where the dislocation substructures are represented by pairs of dislocation density potentialfunctions (DDPFs), denoted by φ and ψ . The slip plane distribution is characterized bythe contour surfaces of ψ , while the distribution of dislocation curves on each slip planeis identified by the contour curves of φ which represents the plastic slip on the slip plane.By using DDPFs, we can explicitly write down an evolution equation system, which isshown consistent with the underlying discrete dislocation dynamics. The system includesi) A constitutive stress rule, which describes how the total stress field is determined in thepresence of dislocation networks and applied loads; ii) A plastic flow rule, which describeshow dislocation ensembles evolve. The proposed continuum model is validated throughcomparisons with discrete dislocation dynamics simulation results and experimental data.As an application of the proposed model, the “smaller-being-stronger” size effect observedin single-crystal micro-pillars is studied. A scaling law for the pillar flow stress σ flow againstits (non-dimensionalized) size D is derived to be σ flow ∼ log( D ) /D . Keywords:
Dislocation density, Crystal plasticity, Continuum model, Size effect,Micro-pillars, Finite elements ∗ Corresponding author
Email address: [email protected] (Yang Xiang)
Preprint submitted to Elsevier October 10, 2018 . Introduction
It is widely agreed that plasticity theories that properly integrate the accumulated knowl-edge in small-scale physics can facilitate the design of high-end materials. The continuumcrystal plasticity (CCP) theories (e.g. Rice, 1971; Peirce et al., 1983; Fleck and Hutchin-son, 1993; Nix and Gao, 1998; Gurtin, 2002) have shown their values in understanding theelasto-plastic behavior of crystals, but they are still phenomenological. On the other hand,the (three-dimensional) discrete dislocation dynamical (DDD) models take the dislocationmicrostructural evolution into account based on the fact that plastic deformation of crystalsis carried out by the motion of a large number of dislocations (e.g. Kubin et al., 1992; Zbibet al., 1998; Fivel et al., 1998; Ghoniem et al., 2000; von Blanckenhagen et al., 2001; Wey-gand et al., 2002; Xiang et al., 2003; Benzerga et al., 2004; Quek et al., 2006; Arsenlis et al.,2007; Rao et al., 2007; El-Awady et al., 2008; Tang et al., 2008; Senger et al., 2008; El-Awadyet al., 2009; Chen et al., 2010; Zhao et al., 2012; Zhou and LeSar, 2012; Ryu et al., 2013; Zhuet al., 2013; Zhu and Chapman, 2014b; Wang et al., 2014). In DDD models, dislocations aretreated as line singularities embedded into an elastic medium. The kinematics of individualdislocations is governed by a collection of laws for dislocation gliding, climb, multiplication,annihilation, reaction, etc., and the microstructural changes within crystals are then cap-tured by the evolution of dislocation curves. DDD models have been well applied to provideinsights in understanding many plastic deformation processes observed in micro- or nano-crystalline structures, such as in thin films and interfaces (e.g. von Blanckenhagen et al.,2001; Weygand et al., 2002; Quek et al., 2006; Zhou and LeSar, 2012; Wang et al., 2014) andin micro-pillars (e.g. Rao et al., 2007; El-Awady et al., 2008; Tang et al., 2008; Senger et al.,2008; El-Awady et al., 2009; Ryu et al., 2013). However, three-dimensional DDD modelsbecome too computationally intensive when the specimen size exceeds the order of severalmicrons.Therefore, a successful dislocation-density-based theory of plasticity (DDBTP) whoseassociated length scale lies between CCP’s and DDD’s is still highly expected. The devel-opment of DDBTP dates back to the works of Nye (1953), where a dislocation network is2epresented by a continuously distributed second-order tensor, known as the Nye disloca-tion tensor. Nowadays with more knowledge in physics taking place on smaller scales, asuccessful DDBTP should be constituted by laws that are consistent with the underlyingdiscrete dislocation dynamics from the following two aspects: i) A constitutive stress ruleto determine the stress field in the presence of a continuous dislocation density distributionand applied loads; ii) A plastic flow rule to capture the motion of dislocation ensembles (inresponse to the calculated stress field), which results in plastic flows in crystals.As the simplest dislocation configuration, systems of straight and mutually parallel dis-locations have been analyzed relatively well at the continuum level (e.g. Groma et al., 2003;Berdichevsky, 2006; Voskoboinikov et al., 2007; Kochmann and Le, 2008; Hall, 2011; Oz-top et al., 2013; Geers et al., 2013; Le and Guenther, 2014; Schulz et al., 2014; Zhu andChapman, 2014a; Le and Guenther, 2015). In this case, each dislocation can be treated asa point singularity in a plane that is perpendicular to all dislocations. As a result, the Nyedislocation density tensor is reduced to several scalar dislocation density functions. Thegeometric complexity of the dislocation networks is dramatically reduced in this case. How-ever, the development of three-dimensional DDBTP is still far from satisfactory despite anumber of valuable works (e.g. Nye, 1953; Kroener, 1963; Kosevich, 1979; Nelson and Toner,1981; Mura, 1987; Head et al., 1993; El-Azab, 2000; Acharya, 2001; Svendsen, 2002; Arsenlisand Parks, 2002; Sedl´aˇcek et al., 2003; Alankar et al., 2011; Sandfeld et al., 2011; Engelset al., 2012; Hochrainer et al., 2014; Li et al., 2014; Cheng et al., 2014). One of the mainbarriers in establishing a successful three-dimensional theory is due to the fact that the com-plex networks of curved dislocation substructures make the upscaling of discrete dislocationdynamics extremely difficult.To overcome such difficulties, Xiang (2009) introduced the idea of a coarse-grained dis-registry function (CGDF), which is defined to approximate the exact disregistry function(plastic slip) used in the Peierls-Nabarro models (Peierls, 1940; Nabarro, 1947; Xiang et al.,2008), by a smoothly varying profile without resolving details of dislocation cores. By thisway, the density distribution of a discrete curved dislocation network in a single slip plane(after local homogenization) can be simply represented by the scalar CGDF (more precisely,3he spatial derivatives of CGDF), and dislocation dynamics on the slip plane at the contin-uum level is explicitly formulated in terms of the evolution of the CGDF (Zhu and Xiang,2010). Using this representation of CGDF for dislocation density distribution, connectivitycondition of dislocations is automatically satisfied. The underlying topological changes ofdislocations are automatically handled by the evolution equation of the CGDF, and no lawfor dislocation annihilation needs to be further imposed. The dynamics of dislocations in thecontinuum model is derived from the DDD model, and the dislocation velocity in the contin-uum model depends on a continuum version of the Peach-Koehler force on dislocations. Ithas been rigorously shown by Xiang (2009) that in the continuum model, the Peach-Koehlerforce due to the resolved shear stress of a family of curved dislocations can be decomposedinto a long-range dislocation interaction force and a short-range self line tension force, andthey can both be expressed in terms of the spatial derivatives of CGDFs. The Frank-Readsource, which is one of the major mechanisms for dislocation multiplication, is also wellincorporated into this continuum framework (Zhu et al., 2014). As one application of thiscontinuum model using CGDFs, a two-dimensional Hall-Petch law, which relates the flowstress of a polycrystal not only to the physical dimension of its constituent grains, but alsoto the grain aspect ratio, is derived without any adjustable parameters (Zhu et al., 2014).In this paper, we generalize our previous single-slip-plane model to that for dislocationensembles in three-dimensions, where the density distribution of dislocations is locally co-determined by an in-plane dislocation density distribution and a slip plane distribution. Totake into account the spatial variation from these two aspects, we define a pair of dislo-cation density potential functions (DDPFs) for each active slip system. One DDPF ψ is employed to describe the slip plane distribution (after local homogenization) by itscontour surfaces, and the other DDPF φ is defined such that φ restricted on each slip planedescribes the plastic slip across the slip plane and identifies the density distribution of dislo-cation curves (after local homogenization) on that plane. Here we name φ and ψ by densitypotential functions, because the Nye dislocation density tensor is represented in terms ofthe spatial derivatives of these two functions. As our previous continuum model in a singleslip plane (Xiang, 2009), the major advantage of this three-dimensional continuum model4ies in its simple representation of distributions of curved dislocations (after local homoge-nization) using two scalar DDPFs, which automatically satisfies the connectivity conditionof dislocations.To derive the constitutive stress rule in the continuum framework with the DDPFs,we sequentially express the Nye dislocation density tensor, the plastic distortion and theelastic strain tensor in terms of the DDPFs. As in our previous continuum model in asingle slip plane (Xiang, 2009), the continuum Peach-Koehler force due to the resolvedshear stress consists of a long-range dislocation interaction force and a short-range self linetension force. The long-range stress field is determined by the derived constitutive stressrule and the equilibrium equations along with boundary conditions, and a finite element(FE) formulation is proposed to compute this long-range stress field. The local self linetension effect can be explicitly formulated in terms of the spatial derivatives of φ and ψ .The plastic flow rule is described by evolution equations of the DDPFs φ and ψ . For face-centered-cubic (fcc) crystals considered in this paper, the motion of dislocations is limitedto their respective slip planes, and the plastic flow rule is given by an evolution equationof φ . Frank-Read sources in our single-slip-plane continuum model (Zhu et al., 2014) isgeneralized to three-dimensional case with continuous distributions of a number of sources.The derived equations form a closed system evolving in time as summarized in Eqs. (54) to(57) in Sec. 2.9.With this continuum model, we investigate the size effect on crystalline strength widelyobserved in the uniaxial compression tests of monocrystalline micro-pillars (e.g. Uchic et al.,2004, 2009; Jang et al., 2012). Practically, an empirical power law is adopted to relatethe pillar flow stress σ flow to the pillar size D by σ flow ∼ D − m , where m is found to befrom 0 to 1, varying by study (Uchic et al., 2009). Typically there are two classes ofmodels proposed to rationalize this size effect. The first type falls into the family of the“dislocation starvation” models (Greer et al., 2005; Greer and Nix, 2006). They arguedthat a crystal small in size does not provide enough space for dislocation multiplicationand the flow strength gets increased as a result. The second category of models attributethe observed size effect to the stochastics of dislocation source lengths in the small-size5pecimens (Parthasarathy et al., 2007). There are also models using statistical approachesto reproduce the power law expression (e.g. Gu and Ngan, 2013). Analysis of heterogeneousdeformation in single crystal micropillars under compression has been performed by usinga hybrid elasto-viscoplastic simulation model which couples DDD model, and a scaling law σ flow ∼ D − n log D + αD − with two parameters α and 0 < n < σ flow ∼ bD log (cid:18) Db (cid:19) . (1)This relation is validated through comparison with experimental data conducted in severalfcc crystals for pillar size ranging from submicrons to tens of microns.The rest of this paper is organized as follows. In Sec. 2, the (three-dimensional) DDPFsare introduced, and this is followed by the derivation of the constitutive stress rule and theplastic flow rule needed at the continuum level. In Sec. 3, numerical schemes for the derivedequation system are presented. In Sec. 4, the derived continuum model is validated throughcomparisons with DDD simulation results. In Sec. 5, the continuum model is applied to studythe size effect arising in the uniaxial compression tests of single crystalline micro-pillars.To better illustrate the derivation of the continuum model, following notations are usedthroughout the article unless specified. The Cartesian coordinates are denoted by r =( x, y, z ). The i -th entry of a vector, for example r , is denoted by r i , and the ij -th entryof a second-order tensor, for example σ , is denoted by σ ij . Unless specified, the followingnotations are used: the vector gradient ( ∇ u ) ij = ∂u i /∂r j ; the cross product ( m × n ) i = P j,k =1 ǫ ijk m j n k with ǫ ijk the permutation tensor; the inner product of two vectors m · n = P i =1 m i n i ; the inner product of two second-order tensors α : β = P i,j =1 α ij β ij ; the innerproduct of a fourth-order tensor and a second-order tensor L : β = P k,l =1 L ijkl β kl ; themagnitude of a vector | u | = √ u · u ; the symmetric part of a second-order tensor sym( α ) =6 α + α T ) /
2; the outer product of two vectors ( a ⊗ b ) ij = a i b j ; and the row “curl” of asecond-order tensor ( ∇ × α ) ij = P k,l =1 ǫ jkl α il,k .
2. Continuum plasticity model based on dislocation density potential functions
In this section, we first review the continuum model for dislocation dynamics in one slipplane using a two-dimensional coarse-grained disregistry function (Xiang, 2009; Zhu andXiang, 2010; Zhu et al., 2014). Then by using the DDPFs, we build the three-dimensionalcontinuum model for dislocation dynamics and plasticity.
We consider describing a given discrete dislocation network in a single slip plane andwith the same Burgers vector b (e.g. the configuration on the top left of Fig. 1(a)) bya dislocation continuum. For this purpose, we use two field quantities at the continuum (a) (b)Figure 1: (a) Turning a discrete dislocation network in a single slip plane into a dislocation continuum.The discrete dislocation network with Burgers vector b is approximated by a family of smoothly varyingdislocation curves with local line direction l and dislocation spacing d in , averaged locally over a representativearea Ω ǫ of the discrete dislocation network, see Eq. (2). (b) Representation of the dislocation continuum in(a) (on the top right) using the CGDF φ . The i -th dislocation curve corresponds to the contour of φ with height φ = ib . l and the average dislocation spacing d in , whichequals the reciprocal of the net dislocation length per area ρ . The values of the two fieldquantities l and d in at each point in the continuum model come from the averaging over arepresentative area Ω ǫ of size ǫ in the discrete model, where d in << ǫ << D with D beingthe sample size in the continuum model. Under this assumption, all dislocations inside Ω ǫ can be treated as line segments, as schematically shown in Fig. 1(a). By superpositioningall dislocation segments inside Ω ǫ , we can obtain a super line segment denoted by L , i.e. L = P i s i for line segment s i ∈ Ω ǫ . The quantities of interest at the continuum level canthen be defined by l = lim ǫ/D → L | L | , d in = 1 ρ = lim ǫ/D → | Ω ǫ || L | , (2)where | Ω ǫ | is the area of Ω ǫ .By performing such averaging process everywhere, the original discrete dislocation net-work is turned into a dislocation continuum as schematically shown in Fig. 1(a) (on thetop right) and we can then use field variables to express it. Xiang (2009) introduced acoarse-grained disregistry function φ , such that the i -th dislocation curve in the averageddislocation continuum corresponds to the contour of φ with height φ = ib , see Fig. 1(b).With the function φ , the local dislocation line direction and inter-dislocation spacing canbe given by l = 1 |∇ φ | (cid:18) ∂φ ∂y , − ∂φ ∂x (cid:19) (3)and d in = 1 ρ = b |∇ φ | , (4)respectively, where ∇ = (cid:16) ∂∂x , ∂∂y (cid:17) . A fact used to derive these formulas is that the normalvector n of the i -th dislocation curve in the averaged dislocation continuum n = ∇ φ |∇ φ | , (5)which is the unit vector in the slip plane and normal to l .8ith the introduction of φ , we can capture the resolved shear stress including that dueto the long-range dislocation interaction by an integral τ = µ π Z R ( x − ˜ x ) ∂φ (˜ x, ˜ y ) ∂ ˜ x + ( y − ˜ y ) ∂φ (˜ x, ˜ y ) ∂ ˜ y (( x − ˜ x ) + ( y − ˜ y ) ) / d˜ x d˜ y + µν π (1 − ν ) b Z R (cid:18) b ∂φ (˜ x, ˜ y ) ∂ ˜ x + b ∂φ (˜ x, ˜ y ) ∂ ˜ y (cid:19) ( b ( x − ˜ x ) + b ( y − ˜ y ))(( x − ˜ x ) + ( y − ˜ y ) ) / d˜ x d˜ y, (6)where µ and ν are the shear modulus and the Poisson’s ratio, respectively, and a contributiondue to the local line tension effect τ = µbκ π ν − ν − ν − ν ( b ∂φ ∂x + b ∂φ ∂y ) /b r(cid:0) ∂φ ∂x (cid:1) + (cid:16) ∂φ ∂y (cid:17) log b/r c π r(cid:0) ∂φ ∂x (cid:1) + (cid:16) ∂φ ∂y (cid:17) + 1 , (7)where r c is a parameter depending on the dislocation core and κ = −∇· ( ∇ φ / |∇ φ | ) is thelocal (average) curvature of the dislocation. These stress formulas in the continuum modelwere derived rigorously from the discrete dislocation model by asymptotic analysis (Xiang,2009).In the single-slip-plane continuum model, φ measures the plastic slip across the slipplane in the direction of the Burgers vector, and the plastic flow is governed by an evolutionequation of φ (Zhu and Xiang, 2010; Zhu et al., 2014): ∂φ ∂t + v n s(cid:18) ∂φ ∂x (cid:19) + (cid:18) ∂φ ∂y (cid:19) = s . (8)Here v n is the dislocation moving speed along its normal direction (Zhu and Xiang, 2010) v n = m g ( τ + τ ) b, (9)where m g is the dislocation glide mobility and the applied stress can also be included inthe above equation. The right-hand term s in Eq. (8) formulates the effect due to thedislocation multiplication by Frank-Read sources (Zhu et al., 2014), which will be reviewedin detail in Sec. 2.8.1. 9 .2. Dislocation substructures in three dimensions represented by dislocation density poten-tial functions From now on, we present our continuum model for the dynamics of dislocation struc-tures in three dimensions. In this subsection, we introduce the representation of dislocationsubstructures in three dimensions using DDPFs. (a) (b)Figure 2: (a) A discrete dislocation network in three dimensions is approximated by a dislocation continuum,after local homogenization of dislocation ensembles within some representative volume Ω ǫ . These dislocationsin the network are associated with the same slip system with Burgers vector b and slip plane normal direction m . The dislocation continuum is characterized by the dislocation line direction l , the slip plane spacing d sl and the in-plane dislocation spacing d in averaged over the representative volume Ω ǫ , see Eqs. (11) and (12).(b) Representation of the dislocation continuum in (a) using a pair of DDPFs φ and ψ . The DDPF ψ isemployed such that the j -th slip plane is the contour plane of ψ of height ψ = jb . The DDPF φ is definedsuch that φ restricted on each slip plane describes the dislocation distribution on that plane, i.e., the i -thdislocation on the slip plane is the contour line φ = ib . We first focus on a given discrete dislocation network in a single slip system with Burgersvector b and slip plane normal direction m . We take a representative cuboid volume Ω ǫ with10ize ǫ , the Nye dislocation density tensor over Ω ǫ is (Nye, 1953; Kroener, 1963; Kosevich,1979) α = P i b ⊗ s i / | Ω ǫ | for dislocation line segment s i ∈ Ω ǫ , where | Ω ǫ | is the volume ofΩ ǫ , see Fig. 2(a). Denoting L = P i s i for all dislocation line segments s i ∈ Ω ǫ which is thesuperposition of all dislocation segments inside Ω ǫ , and assuming that ǫ << D where D isthe domain size of the continuum model, we write the Nye dislocation density tensor overΩ ǫ as α = ρ g b ⊗ l , (10)where l is the average dislocation line direction over Ω ǫ and ρ g is the dislocation densityover Ω ǫ (net length per unit volume): l = lim ǫ/D → L | L | , ρ g = lim ǫ/D → | L || Ω ǫ | . (11)See Fig. 2(a) for an illustration of this average. We characterize ρ g by two variables ρ g = 1 d sl d in , (12)where d sl is the slip plane spacing and d in is the in-plane dislocation spacing, see Fig. 2(a).By performing such averaging process everywhere, the original three-dimensional discretedislocation network is turned into a dislocation continuum in three dimensions as schemati-cally shown in Fig. 2(a). We introduce a pair of dislocation density potential functions (DDPFs) φ and ψ to represent this resulting dislocation continuum consisting of dislocationsin the same slip system. The DDPF ψ is employed to describe the distribution of the slipplanes in the averaged dislocation continuum, such that the j -th slip plane is the contourplane of ψ with ψ = jb . In this paper, we focus on fcc crystals, in which dislocations stay intheir respective slip planes. We further assume a uniform distribution of the slip planes withprescribed slip plane spacing d sl . Under these conditions, the DDPF ψ for the slip systemwith Burgers vector b and slip normal direction m can be written as ψ ( r ) = b m · ( r − r ) d sl , (13)where r is some point on the 0-th slip plane, see Fig. 2(b).11ith slip planes in the averaged dislocation continuum described by ψ , we introduceanother DDPF φ to describe the distribution of dislocations on each slip plane, such that φ restricted on each slip plane is the two-dimensional CGDF φ in the single-slip planemodel reviewed in the previous subsection. That is, the i -th dislocation curve on the j -thslip plane in the averaged dislocation continuum can be represented by { r | φ ( r ) = ib, ψ ( r ) = jb, for integers i, j } , (14)see Fig. 2(b). As φ in the single-slip plane model, here φ restricted on each slip planedescribes the plastic slip across the slip plane in the direction of the Burgers vector. Localgeometric quantities of dislocations and the Nye dislocation density tensor in the averageddislocation continuum are simply expressed in terms of these DDPFs φ and ψ , see thefollowing subsections.The major advantage of this three-dimensional continuum model characterized by DDPFslies in its simple representation of distributions of curved dislocations (after local homog-enization). This representation also automatically satisfies the connectivity condition ofdislocations, see Sec. 2.4. In addition, dislocation annihilation within the same slip plane isautomatically handled as in the previous single-slip-plane model (Xiang, 2009; Zhu and Xi-ang, 2010). When there are multiple slip systems activated, one can assign a pair of DDPFsto each active slip system.In body-centered-cubic (bcc) crystals or fcc crystals at high temperatures where dislo-cations are able to move out of their slip planes by climb, dislocations no longer stay onplanar slip planes. Dislocation networks in these crystals can still be represented under thecontinuum framework characterized by the DDPFs φ and ψ . In these cases, the contoursurface ψ ( r ) = jb may describe a curved surface rather than a plane, and the local slip planenormal direction in the averaged dislocation continuum is determined by m = ∇ ψ |∇ ψ | . (15) Now we present expressions of the local geometric quantities of dislocations in the aver-aged dislocation continuum in terms of the DDPFs φ and ψ . The slip plane normal direction12s given by Eq. (15). The local dislocation line direction l is l = ∇ φ × ∇ ψ |∇ φ × ∇ ψ | . (16)This is because the dislocation is contained in both the contour surfaces of φ = ib and ψ = jb (see Eq. (14)), thus it is perpendicular to both ∇ φ and ∇ ψ . The local dislocation normal n (which is in the slip plane) is also calculated in terms of DDPFs by n = m × l = ∇ ψ × ( ∇ φ × ∇ ψ ) |∇ φ × ∇ ψ ||∇ ψ | . (17)Moreover, by the Frenet-Serret formulas, we have κ n = ( l · ∇ l , l · ∇ l , l · ∇ l ), where κ isthe curvature of the local dislocation. Thus with the expressions for l and n in Eqs. (16)and (17), κ is also represented by the DDPFs as κ = n · ( κ n ) = X i =1 n i l · ∇ l i . (18)Finally, the in-plane dislocation spacing and the slip plane spacing are respectively d in = b |∇ φ × m | , d sl = b |∇ ψ | . (19)Here we have used the definition in Eq. (14) and the fact that the gradient of φ over the slipplane is ∇ in-plane φ = ∇ φ − ( m · ∇ φ ) m = ∇ φ × m . Note that the above formulas hold for ageneral DDPF ψ for the definition of dislocations in Eq. (14), including but not limited tothe particular expression in Eq. (13). After the local homogenization process described in Sec. 2.2, the Nye dislocation den-sity tensor α for a single slip system is calculated by Eq. (10) in terms of the averageddislocation line direction l and the net length per volume dislocation density ρ g over somerepresentative volume, and ρ g is given by Eq. (12) using the slip plane spacing d sl and thein-plane dislocation spacing d in . Further using the formulas in Eqs. (15), (16) and (19) withthe DDPFs φ and ψ , we have ρ g = |∇ φ × m | · |∇ ψ | b = |∇ φ × ∇ ψ | b . (20)13herefore, incorporating Eqs. (16) and (20) into Eq. (10) gives the following expressionfor the Nye dislocation density tensor based on the DDPFs α = b b ⊗ ( ∇ φ × ∇ ψ )= 1 b b (cid:16) ∂ψ∂z ∂φ∂y − ∂ψ∂y ∂φ∂z (cid:17) b (cid:0) ∂ψ∂x ∂φ α ∂z − ∂ψ∂z ∂φ∂x (cid:1) b (cid:16) ∂ψ∂y ∂φ∂x − ∂ψ∂x ∂φ∂y (cid:17) b (cid:16) ∂ψ∂z ∂φ∂y − ∂ψ∂y ∂φ∂z (cid:17) b (cid:0) ∂ψ∂x ∂φ∂z − ∂ψ∂z ∂φ∂x (cid:1) b (cid:16) ∂ψ∂y ∂φ∂x − ∂ψ∂x ∂φ∂y (cid:17) b (cid:16) ∂ψ∂z ∂φ∂y − ∂ψ∂y ∂φ∂z (cid:17) b (cid:0) ∂ψ∂x ∂φ∂z − ∂ψ∂z ∂φ∂x (cid:1) b (cid:16) ∂ψ∂y ∂φ∂x − ∂ψ∂x ∂φ∂y (cid:17) . (21)For the case with multiple slip systems, we use one pair of DDPFs φ λ and ψ λ for the λ -th slip system, and the Nye dislocation density tensor is expressed by α = X λ b λ ( b λ ) ⊗ ( ∇ φ λ × ∇ ψ λ ) . (22)It is easy to check that α given by Eq. (22) satisfies ∇ · α = , (23)which means P j =1 ∂α ij ∂r j = 0 for all i = 1, 2 and 3. This is the connectivity condition ofdislocations, which is due to the fact that dislocation lines can never begin or end inside thesample. In this subsection, we present the constitutive stress rule in our continuum plasticitymodel using DDPFs. The derivation is based on the constitutive relations commonly used inclassical dislocation-density-based models, including (when the specimen experiences smalldeformations): • The total distortion, which is the gradient of the displacement field u , can be decom-posed into an elastic distortion β e and a plastic distortion β p ∇ u = β e + β p . (24) • The Nye dislocation density tensor equals the curl gradient of the plastic distortion ∇ × β p = − α = − X λ b λ ( b λ ) ⊗ ( ∇ φ λ × ∇ ψ λ ) . (25)14 The stress field σ satisfies the Hooke’s law (e.g. the isotropic case): σ = 2 µ ǫ e + 2 µν − ν tr( ǫ e ) I , (26)where ǫ e is the elastic strain tensor related to the elastic distortion by ǫ e = sym ( β e ),tr( ǫ e ) is the trace of ǫ e , and I is the 3 × • The equilibrium condition in absence of body forces is ∇ · σ = . (27)By generalizing the expression for β p for the case of discrete dislocations (e.g., Sec. 1.6of Mura, 1987), and using the fact that φ describes the plastic slip across the slip plane inthe direction of the Burgers vector and ψ describes the distribution of the slip planes, in thecontinuum model we have β p = − X λ φ λ ( b λ ) ( b λ ⊗ ∇ ψ λ ) . (28)It can be checked that β p given by Eq. (28) satisfies Eq. (25). A comparison between Eq. (28)and its counterpart in continuum crystal plasticity theories suggests that the scalar γ λ = φ λ |∇ ψ λ | b λ (29)measures the magnitude of the plastic shearing in the direction of Burgers vector b λ in the λ -th slip system.By using Eqs. (24) and (28), we can express the elastic strain tensor by ǫ e = sym( ∇ u ) + X λ φ λ ( b λ ) sym( b λ ⊗ ∇ ψ λ ) . (30)Incorporating Eq. (30) with the Hooke’s law (26) and using the fact that b λ · ∇ ψ λ = 0 (theBurgers vector is always perpendicular to the normal direction of the slip plane), we obtain σ = L : ∇ u + 2 µ X λ φ λ ( b λ ) sym( b λ ⊗ ∇ ψ λ ) , (31)where the symmetric fourth-order tensor L is defined such that L : ∇ u = 2 µ (cid:18) sym ( ∇ u ) + ν − ν ( ∇ · u ) I (cid:19) . (32)15hen a solid body is purely elastic (dislocation free), φ λ = 0 and Eq. (31) becomes σ = L : ∇ u , which is exactly the constitutive relation in the linear elasticity theory.In summary, with given dislocation distributions described by DDPF pairs φ λ and ψ λ for λ varies over all the slip systems, the stress field σ is obtained by solving Eq. (31) andthe equilibrium condition in Eq. (27) over the sample Ω, subject to appropriate boundaryconditions on the boundary ∂ Ω. One boundary condition is the displacement boundarycondition imposed on part of the boundary denoted by ∂ Ω d : u | ∂ Ω d = u b , (33)another boundary condition is the traction boundary conditions imposed on the rest of theboundary denoted by ∂ Ω t : σ | ∂ Ω t · k = t b , (34)with k the outer unit normal to the surface ∂ Ω t . Here ∂ Ω = ∂ Ω d ∪ ∂ Ω t . It is well-known that there are various short-range dislocation interactions that playimportant roles in the plastic deformation processes, in addition to the long-range effectof dislocations. In this subsection, we first show that the continuum stress formulationpresented in the previous subsection is for the long-range dislocation interaction. We thenpresent the formulation for the short-range line tension effect in the DDPFs based threedimensional continuum model.When the medium is the whole three dimensional space R , using the Nye dislocationdensity tensor formula in Eq. (22) in our continuum model, the stress formulas given by theDDD model and the continuum model are as follows (Hirth and Lothe, 1982; Mura, 1987).In the DDD model, the stress field is σ dd ( r ) = µ π X λ N λ X i =1 Z γ λi sym (cid:18) b λ × ( r − ˜ r ) | r − ˜ r | ⊗ l d s (cid:19) + µ π (1 − ν ) X λ N λ X i =1 Z γ λi ( l · ( b λ × ∇ ))( ∇ ⊗ ∇ − I ∇ ) | r − ˜ r | d s, (35)16here γ λi denotes the i -th dislocation in the λ -th slip system, N λ is the number of dislocationsof the λ -th slip system, ˜ r goes over all points on γ λi , s is the arclength of γ λi , and ∇ denotesthe gradient with respect to r . In the continuum framework characterized by the DDPFs,the stress field is σ ( r ) = X λ µ π ( b λ ) Z R sym (cid:18) b λ × ( r − ˜ r ) | r − ˜ r | ⊗ ( ∇ φ λ (˜ r ) × ∇ ψ λ (˜ r )) (cid:19) d ˜ V + X λ µ π ( b λ ) (1 − ν ) Z R ( b λ · ∇ φ λ (˜ r ))( ∇ ψ λ (˜ r ) · ∇ )( ∇ ⊗ ∇ − I ∇ ) | r − ˜ r | d ˜ V , (36)where d ˜ V is an infinitesimal volume associated with ˜ r . Since the Nye dislocation density ofthe DDD model is averaged into the Nye dislocation density of the continuum model in thecoarse-graining process as described in Secs. 2.2 and 2.4, from Eqs. (35) and (36), we have σ dd −→ σ (37)in the coarse-graining process from the DDD model to the continuum model.When we consider a material with finite size, an image stress is added in the DDDmodel to accommodate the boundary conditions (Van der Giessen and Needleman, 1995).Following the same coarse-graining process, such a DDD stress solution gives a stress solutionthat satisfies all the equations and boundary conditions of the continuum model, which isthe same as that obtained by directly solving the continuum model due to the uniquenessof the solution. Therefore in this case, the stress field solved from the continuum model isalso the leading order approximation of the stress field of the DDD model.The stress field σ solved from the elasticity system in the continuum model describes thelong-range elastic interaction of dislocations. There are also various short-range dislocationinteractions that play important roles in the plastic processes. How to capture these short-range effects at the continuum level is an important issue in the development of dislocationbased plasticity theories.One way to incorporate the short-range dislocation interactions in the continuum modelis to add complementary resolved shear stresses that describes these effects, because thedynamics of dislocations depends on the resolved component of the stress field in the slip17lane (Hirth and Lothe, 1982). For the λ -th slip system, the resolved shear stress due tothe long-range dislocation interaction in the continuum model is τ λ long = b λ b λ · σ · ∇ ψ λ |∇ ψ λ | . (38)Recall that m λ = ∇ ψ λ / |∇ ψ λ | is the normal direction of the slip plane. The principle toinclude additional shear stresses in the continuum model is to obtain a better approximationto the shear stress in the DDD model τ λ dd −→ τ λ total = τ λ long + τ λ short (39)in the coarse-graining process from the DDD model to the continuum model, where τ λ short includes contributions to the resolved shear stress due to all the important short-rangeeffects, and τ λ dd = b λ b λ · σ dd · m λ is the resolved shear stress of the λ -th slip system using theDDD model.In this paper, we focus on the dislocation line tension effect which is one of the importantshort-range dislocation effects, and the total shear stress in the continuum model is τ λ total = τ λ long + τ λ self , (40)where τ λ self is the contribution to the resolved stress due to the line tension effect. The dis-location line tension effect plays crucial roles for example in dislocation multiplication byFrank-Read sources. For dislocation distributions in a single slip plane, we have used asymp-totic analysis to rigorously derive the expression of the line tension effect in the continuummodel from the DDD model (Xiang, 2009) as reviewed in Sec. 2.1. Here for dislocationdistributions in three dimensions represented by DDPFs, for the λ -th slip system, τ λ self isgiven by τ λ self = µb λ π (cid:18) ν − ν − ν − ν · |∇ ψ λ | ( b λ · ∇ φ λ ) ( b λ ) |∇ ψ λ × ∇ φ λ | (cid:19) κ · log (cid:18) b λ |∇ ψ λ | πr c |∇ ψ λ × ∇ φ λ | + 1 (cid:19) , (41)where the curvature κ of the local dislocation is calculated by Eq. (18). In order to derivethis formula, we have used the fact that the DDPF φ restricted onto each slip plane whichis a contour surface of the DDPF ψ reduces to our previous single slip plane model as wellas the expressions of geometric characterizations of the local dislocation given in Sec. 2.3.18here are other short-range effects to be included in τ short , such as that due to the short-range interaction between dislocations from different slip systems. These will be discussedelsewhere. In our continuum model characterized by DDPFs, the plastic flow rule is given by˙ φ λ + v λ n |∇ φ λ × ∇ ψ λ ||∇ ψ λ | = s λ (42)for the λ -th slip system, where ˙ φ λ = ∂φ λ /∂t , v λ n is the speed of the local dislocation in itsnormal direction in the slip plane, and s α formulates the dislocation generation by Frank-Read sources to be discussed in details in the next subsection. Being the three-dimensionalversion of Eq. (8), Eq. (42) is also established based on the conservation of plastic shearslips (Zhu and Xiang (2010), see also the level set DDD method in Xiang et al. (2003)) andthe fact that the gradient of φ λ over the slip plane is ∇ φ λ × m λ , where m λ = ∇ ψ λ / |∇ ψ λ | is the normal direction of the slip planes. Note that there is no need to assign extra rulefor dislocation annihilation on the same slip plane, which is automatically handled by thetopological changes in the contours of the DDPF φ λ during its evolution.The dislocation speed v λ n in Eq. (42) is determined by a mobility law following the DDDmodels: v λ n = m g b λ τ λ total = m g b λ (cid:0) τ λ long + τ λ self (cid:1) , (43)where m g is the dislocation glide mobility, τ λ long is the component of the long-range stressfield σ determined by the elasticity problem in Sec. 2.5 and resolved in the λ -th slip systemgiven in Eq. (38), and τ λ self is the effective stress due to the line tension effect given in Eq. (41)in the previous subsection.In general, boundary conditions are needed for Eq. (42). One extreme case is thatdislocations can exit the specimen freely, which can be described by the Neumann boundarycondition ∂φ λ ∂ k s = 0 where k s is the outer normal direction of the intersection of the slipplane and the specimen surface ( k s = m × ( m × k ) / | m × ( m × k ) | where k is the outer19ormal direction of the specimen surface). Another extreme case is that the dislocations areimpenetrable to a specimen surface, where φ λ is fixed on the specimen surface.The rate of Nye dislocation density tensor ˙ α , the rate of the plastic distortion ˙ β p andthe rate of the scalar plastic shear ˙ γ λ can all be easily calculated using ˙ φ λ in Eq. (42) andthe expressions of α , β p and γ λ in terms of φ λ and ψ λ in Eqs. (22), (28) and (29). Notethat in this paper, the DDPF ψ λ that describes the distribution of the slip planes is fixed.Many other quantities that are useful in understanding the plastic behavior of crystalscan also be expressed in terms of the DDPFs φ and ψ . For example, the total dislocationdensity (length per unit volume) within the specimen Ω can be formulated by ρ tot = 1 | Ω | X λ Z Ω |∇ φ λ × ∇ ψ λ | ( b λ ) d V + ρ , (44)where ρ is the dislocation density due to the initial distribution of dislocation segments ofthe Frank-Read sources (see the next subsection). Moreover, the total plastic strain rate,which is conventionally defined to be the rate of area swept by all dislocations multiplied bythe length of the respective Burgers vector per volume, can be expressed by˙ ǫ ptot = 1 | Ω | X λ Z Ω ˙ φ λ |∇ ψ λ | b λ d V. (45) In this subsection, we present the expression for the source term s λ in the plastic flowrule in Eq. (42) for a continuous distribution of Frank-Read sources in three dimensions,which is generalized from our previous continuum model for individual sources on a singleslip plane. We first review our continuum model for individual sources on a single slip plane (Zhuet al., 2014).A Frank-Read source is a dislocation segment pinned at its two ends. When the resolvedshear stress τ acting on it exceeds a critical value, known as the activation stress τ c , it will20eep injecting dislocation loops to the system (Hirth and Lothe, 1982). The time it takesfor a Frank-Read source to perform an operating cycle is known as the nucleation time t nuc .However, if observed at the continuum level, what we see is not the detailed loop-releasingprocess, but continuous dislocation flux originating from a small source region. In the single-slip plane model (Zhu et al., 2014) reviewed in Sec. 2.1, the operation of a Frank-Read sourceis controlled by three parameters all determined by the underlying DDD model: the sourceactivation stress τ c , the source operating rate which equals to 1 /t nuc , and the source regionΩ which is the two-dimensional region enclosed by a newly released dislocation loop.The nucleation stress τ c is evaluated by adopting the critical stress formula given by τ c = C s µb πl log (cid:18) lr c (cid:19) , (46)where C s depends on the source character and the Poisson’s ratio ν (with ν = 1 / C s = 1for an edge-oriented source and C s = 1 . l is the length of theFrank-Read source, and r c is a parameter depending on the dislocation core.The nucleation time t nuc is calculated to be t nuc = Q ch lm g b ( | τ | − τ c ) , (47)where τ is the resolved shear stress due to the long-range stress field, and Q ch dependsonly on the source orientation fitted from the DDD simulation ( Q ch = 6 . Q ch = 3 . l parallel to the x axis and centered at ( x s , y s ), the source region is approximately boundedby an ellipse Ω = (cid:26) ( x, y ) (cid:12)(cid:12)(cid:12)(cid:12) ( x − x s ) ( a l ) + ( y − y s ) ( a l ) ≤ (cid:27) , (48)where a and a are calculated to be 2 . . τ c , t nuc and Ω determined by Eqs. (46), (47) and (48), the Frank-Read sourceis incorporated into the single-slip plane model in the following sense. When the resolvedshear stress τ > τ c , a Frank-Read source keeps changing the value of φ at a speed of 1 /t nuc such that dislocation loops enclosing area Ω are continuously generated into the system.21athematically, the source term s in the single-slip plane continuum model in Eq. (8) isgiven by s = − m g b ( τ − sign( τ ) τ c ) Q ch l H ( | τ | − τ c ) · χ Ω , (49)where H ( z ) is the Heaviside function that equals 1 when z > χ Ω isthe characteristic function in Ω , i.e. χ Ω is 1 in Ω and vanishes elsewhere. Now we derive the expression of the Frank-Read sources in the three-dimensional contin-uum model. We first write down a formulation for individual Frank-Read sources in threedimensions by generalizing the single-slip plane model reviewed above, and then derive aformulation based on a source continuum for the three-dimensional continuum model. Theformulation is affiliated with a slip system, and we omit the slip system superscript λ forsimplicity of notations.We first write down a formulation for individual Frank-Read sources in three dimensions.Following the construction of three-dimensional dislocation continuum described in Sec. 2.2,we generalize the source region Ω in the single slip plane model given in Eq. (48) to threedimensions by a cylinder Ω = Ω × [ − d sl / , d sl / m , i.e. the intersection of the cylinderΩ with a slip plane is Ω , and the height of the cylinder is d sl in the slip plane normaldirection m , where d sl is recalled to be the averaged slip plane spacing.Since the physical dimension of the source region Ω is much smaller compared to thatof the domain size D , we further approximate the source region Ω by a regularized pointsource with same volume at the continuum level, i.e. χ Ω ( · ) ≈ | Ω | δ reg ( · ), where δ reg ( · )is a regularized Dirac function of the source region and | Ω | = πa a l d sl is the volumeof Ω . When there are S Frank-Read sources operating, a preliminary way to incorporatethem into the continuum model in Eq. (42) is to express the source term s as the sum ofcontributions from all the individual Frank-Read sources located at r k s , k = 1 , , · · · , S : s ind = − m g πa a b d sl S X k =1 l k ( τ − sign( τ ) τ kc ) Q k ch H ( | τ | − τ kc ) δ reg ( r − r k s ) . (50)22owever, in the continuum model, we do not want to resolve individual Frank-Readsources. For this purpose, we introduce a source continuum to approximate the collectiveeffect of all the individual Frank-Read sources. This is achieved by the following source termin the continuum model in Eq. (42): s = g ( r )( τ − sign( τ ) τ ( r )) H ( | τ | − τ ( r )) , (51)where at any point r , τ ( r ) is the source activation stress at r , and ( τ − τ ( r )) g ( r ) measuresthe rate of plastic shear slip initiated at r by the source continuum.Now we determine the functions τ ( r ) and g ( r ) from the discrete source model in Eq. (50).When the shear stress τ exceeds the critical stresses of all the sources that influence point r , all the Heaviside functions in Eqs. (50) and (51) can be dropped, a comparison betweenthem gives g ( r ) = − πm g a a b d sl X k l k Q k s δ reg ( r − r k s ) (52)and τ ( r ) = − m g µa a b d sl g ( r ) · X k (cid:18) C k s δ reg ( r − r k s ) Q k ch log (cid:18) l k r c (cid:19)(cid:19) , where g ( r ) = 0 . (53)Here we have used the expression of τ c in Eq. (46). The source continuum formula of s inEq. (51) with these expressions of τ ( r ) and g ( r ) still provides a good approximation to s ind in Eq. (50) when the shear stress τ does not exceed the critical stresses of all the sourcesthat influence point r , and it becomes exact again when the shear stress τ falls below all thecritical stresses of the sources that influence point r (In this case, s ( r ) = s ind ( r ) = 0).The functions τ ( r ) and g ( r ) for the source continuum can be calculated from the ar-rangement of the discrete Frank-Read sources. In this paper, the locations and parametersof the Frank-Read sources are given initially and do not change in the simulation, and wecalculate these two functions only once for the initial distribution of the sources.An example of the Frank-Read source continuum is given in Fig. 3. Fig. 3(a) shows a set ofindividual Frank-Read sources with their positions and lengths generated randomly followinga normal and a uniform distribution, respectively, and the profiles of its corresponding g ( r )on several selected slip planes are drawn in Fig. 3(b). It can be observed that g ( r ) attains23 yx z (a) −0.1 0 0.1 −0.10 0.1−0.500.5 y x z (b)Figure 3: Converting a number of Frank-Read sources to a source continuum in a cuboid sample. (a) Thediscrete sources. (b) The corresponding function g ( r ) in the source continuum obtained by Eq. (52) onseveral slip planes. The unit of the color bar is m g b . a relatively high value around the cuboid center, because the number density of individualsources are high in the same place in Fig. 3(a). To summarize, the derived continuum model based on the DDPFs is constituted mainlyby the following equations.1. A constitutive stress rule: Given a dislocation substructure described by φ λ and ψ λ ,the long-range stress field σ is determined by solving σ = 2 µ sym( ∇ u ) + ν tr( ∇ u )1 − ν I + X λ φ λ ( b λ ) sym( b λ ⊗ ∇ ψ λ ) ! (54)and the equilibrium condition ∇ · σ = , (55)with boundary conditions given in Eqs. (33) and (34).2. A plastic flow rule: The motion of dislocations belonging to the λ -th slip system isdescribed by an evolution equation for φ λ ∂φ λ ∂t + v λ n |∇ φ λ × ∇ ψ λ ||∇ ψ λ | = g λ ( r ) ( τ λ − sign( τ λ ) τ λ ( r )) H ( τ − τ ( r )) , (56)24here the dislocation speed v λ n is determined by the mobility law v λ n = m g b λ τ λ total = m g b λ (cid:0) τ λ long + τ λ self (cid:1) , (57)and the shear stress component resolved in the λ -th slip system τ λ long is calculatedfrom the long-range stress field σ by Eq. (38), τ λ self is the contribution due to the localdislocation line tension effect given in Eq. (41), and g λ ( r ) and τ λ ( r ) are two functionsfor the Frank-Read source continuum given by Eqs. (52) and (53), respectively.In addition, the Nye dislocation density tensor α , the total dislocation density in thesystem ρ tot , and the total plastic strain rate ˙ ǫ ptot are expressed in terms of the DDPFs φ and ψ by Eqs. (22), (44), and (45), respectively. In our DDPF-based continuum model, we can write a total free energy of the system as E = E elastic + E self , (58)where E elastic is the elastic energy E elastic = 12 Z Ω ∇ u + X λ φ λ b λ ( b λ ) ⊗ ∇ ψ λ ! : L : ∇ u + X λ µφ λ ( b λ ) sym( b λ ⊗ ∇ ψ λ ) ! d V − Z Ω t t · u d S (59)and E self is the energy due to the dislocation line tension effect E self = X λ Z Ω µ |∇ ψ λ × ∇ φ λ | π (cid:18) ν − ν b λ · ∇ φ λ ( b λ | m λ × ∇ φ λ | ) log b λ πr c | m × ∇ φ λ | (cid:19) d V. (60)This total free energy in Eq. (58) in its form agrees with existing theories using otherrepresentations of dislocation densities (e.g. Nelson and Toner, 1981; Berdichevsky, 2006;Le and Guenther, 2014) and our continuum dislocation dynamics model in a single slipplane (Xiang, 2009).Since the DDPFs ψ λ for all λ stay constant during the deformation process, we can write E = E ( u , φ , φ , · · · , φ λ , · · · ). It can be calculated that inside Ω, we have δ E δ u = −∇ · σ , (61)25nd δ E δφ λ = |∇ ψ λ | b λ (cid:0) τ λ long + τ λ self (cid:1) , (62)Recall that the resolved shear stress τ λ long is from the long-range stress field σ given byEq. (38), and τ λ self is the contribution due to the local dislocation line tension effect given inEq. (41). Note that the variation of E self with respect to φ λ gives the correct leading ordercontribution in τ λ self .Using these variations of the free energy, the governing equations in our continuum modelcan be written as = δ E δ u , (63) ∂φ λ ∂t = − L λ δ E δφ λ + s λ . (64)The first equation gives the equilibrium condition in Eq. (27), which together with theconstitutive relation in Eq. (31) and the boundary conditions in Eqs. (33) and (34) form theelasticity system. (These boundary conditions can also be obtained from the variation of theelastic energy in Eq. (59).) The second equation describes the plastic flow rule (dynamicsof dislocations), where f φ λ = − δ E δφ λ can be understood as the configurational force and L λ = m g ( b λ ) |∇ φ λ ×∇ ψ λ ||∇ ψ λ | is the kinetic coefficient. These equations are consistent with thegoverning equations in the DDD models (Hirth and Lothe (1982) and those references inthe introduction).When the system evolves to its equilibrium state, δ E /δφ λ = 0, which gives rise to themicro force balance state τ λ long + τ λ self = 0. This together with Eq. (63) agree with the modelsfor equilibrium states of distributions of straight dislocations (e.g. Berdichevsky, 2006; Leand Guenther, 2014).
3. Numerical implementation of the continuum model
In this section, we briefly discuss the numerical implementation of our continuum modelas summarized in Sec. 2.9. We will focus on solving for the long-range stress field σ from26qs. (54) and (55) with boundary conditions and the plastic flow rule described by theevolution of DDPF φ λ in Eq. (56).In the simulations in this paper, the computational domain is chosen to be a cuboidΩ = [ − D/ , D/ × [ − D/ , D/ × [ − L/ , L/ z direction, whichis the loading direction, and on the top surface u | z = L/ = u b0 ( t ) is imposed as a result ofcompression. The shear force is set to be free on these two surfaces. On the other four sidesurfaces, traction free boundary conditions are imposed. The DDPF ψ λ in this paper takesthe form in Eq. (13) with inter-slip plane distance d sl = 100 b , and does not change in theplastic deformation process. The long-range stress field satisfying Eqs. (54) and (55) subject to boundary conditionsin Eqs. (33) and (34) yields a weak form that R Ω σ long : ( ∇ v )d V = R ∂ Ω t t b · v d S for any testvector function v ∈ { v | v = , on ∂ Ω d } . Replacing the stress field by the constitutive stressrule given by Eq. (54), we obtain the weak form for the displacement field u Z Ω ∇ v : L : ∇ u d V = Z ∂ Ω t t · v d S − µ X λ Z Ω φ λ sym( b λ ⊗ ∇ ψ λ ) : ( ∇ v )d V, (65)where the forth-order tensor L is defined in Eq. (32).In our simulations, Ω is meshed by C3D8 bricks. We then discretize the weak form inEq. (65) to get a linear algebraic equation system as K FE u FE = f FE , (66)where u FE is a vector of 3 N dimensions containing all nodal values of u with N beingthe total number of nodes, K FE is known as the stiffness matrix, and f FE is assembled bydiscretizing the right hand side of Eq. (65).If the last term is omitted, Eq. (65) is the weak formulation that is widely used in the FEformulation in linear elasticity. Hence many tools well-developed for solving purely elasticproblems, such as meshing, assembling and inversion of K FE , can be inherited by the FE27ormulation proposed here. The contribution from φ λ and ψ λ in Eq. (66) can be treated thesame way as a body force to the linear system. φ λ We discretize the evolution equation of DDPF φ λ in Eq. (56) using a finite differencescheme. Here the grid points of φ λ are chosen coinciding with the vertices of the C3D8bricks. Combined with the mobility law in Eq. (57), Eq. (56) can be written as ∂φ λ ∂t + m g b λ (cid:0) τ λ long + τ λ self (cid:1) | m λ × ∇ φ λ | = s λ . (67)The temporal derivative of φ λ is approximated by the forward Euler scheme. For the spatialderivatives of φ λ , we follow the methods in Xiang et al. (2003) that the first-order upwindscheme is used for those in the term associated with the long-range stress field τ λ long | m λ ×∇ φ λ | , and the central difference scheme is used to calculate those spatial derivatives in theterm τ λ self | m λ × ∇ φ λ | where τ λ self is given by Eq. (41).Moreover, the regularized δ -function in the source term s λ (Eqs. (51)-(53)) is given by δ reg ( r ) = 1∆ s · ππ − (cid:18) cos π | m λ × r | ∆ s + 1 (cid:19) · s (cid:18) cos π ( m λ · r )∆ s + 1 (cid:19) (68)for all r ∈ { r || m λ × r | < ∆ s , | m λ · r | < ∆ s } , where ∆ s and ∆ s are two smoothingparameters.
4. Numerical examples
In this section, the derived continuum model is validated through comparisons with DDDsimulations. All results presented are obtained by using 10 × ×
20 C3D8 bricks in thefinite element discretization of the cuboid simulation domain described at the beginning ofSec. 3.
This example is aimed to provide an illustration of the continuum model. For the cuboidsimulation cell Ω, L = 24000 b and D = 9600 b . A Frank-Read source of length l = 400 b . × − µ is placed at the center of Ω, and a 0 .
3% constant strainis applied by compression on the top surface. All simulations start with a dislocation-freestate. −0.2 0 0.2 −0.200.2−0.500.5 yx z (a) t = Lm g µb −0.2 0 0.2 −0.200.2−0.500.5 yx z (b) t = Lm g µb −0.2 0 0.2 −0.200.2−0.500.5 yx z (c) t = Lm g µb −0.2 0 0.2 −0.200.2−0.500.5 yx z (d) t = Lm g µb Figure 4: Snap shots of the distribution of dislocation curves generated by an operating Frank-Read sourcethat locates at the center of the simulation domain Ω. The source releases dislocation loops in response toa constant strain by compression on the top surface of Ω. The length unit of Ω is its height L . In Fig. 4, we plot the contour lines of the DDPF φ on one of the slip planes, which giverough locations of the dislocation curves. It is observed that in response to the applied strain,the source keeps releasing dislocation loops, which exit Ω from its side surfaces. As a result,the resolved shear stress drops during this loop-releasing process, and so does the pressureon the top surface as shown in Fig. 5. These findings agree with the common understandingabout the role played by a Frank-Read source: It releases dislocation loops so as to softenthe materials. When the time t is about 14000 L/ ( m g µb ), the resolved shear stress finallyfalls below the source activation stress, the source is then deactivated, see Fig. 5.The displacement vector u on the domain boundary ∂ Ω is shown in Fig. 6. It can beseen that the orientations of the displacement gradients are in the normal direction of theactivated slip plane that contains the Frank-Read source, and they are localized near thisslip plane. This is because the dislocation loops that have left the specimen form surface29
500 1000 1500 200002468 x 10 −3 Time (Ld/( µ m g b )) S t r e ss ( µ ) Resolved shear stressat the source centerPressure on top surface
Figure 5: As the Frank-Read source keeps releasing dislocation loops, both the resolved shear stress atthe source and the pressure on the top surface of the simulation domain drop. When the time t is about14000 L/ ( m g µb ), the resolved shear stress falls below the source activation stress (indicated by the dashed-dotted line), and the source is thus deactivated. −0.2 0 0.2 −0.20 0.2−0.500.5 yx z −1−0.500.51x 10 −3 (a) u −0.2 0 0.2 −0.20 0.2−0.500.5 yx z −505x 10 −4 (b) u −0.2 0 0.2 −0.20 0.2−0.500.5 yx z −3−2.5−2−1.5−1−0.50x 10 −3 (c) u Figure 6: The three components of the displacement vector u on the simulation cell surface ∂ Ω. The lengthunit is L , the height of Ω. steps on the activated slip plane (in a smooth sense). To further validate the continuum model, we compare its numerical results with the DDDsimulation results obtained by El-Awady et al. (2008). Following their DDD simulations,30e choose a micro-pillar of nickel. The loading axis is < > , and a single slip systemis activated with slip direction [0¯11] and slip normal (111). The Schmid’s factor m s is thus0.4050. The shear modulus is 76GPa, the Poisson’s ratio is 0 .
31, and the length of theBurgers vector is b = 0 . − . The dislocation glide mobility m g is unspecified by El-Awady et al. (2008), we here follow Senger et al. (2008) to let m g = 10 / (Pa · s). In our simulations, the micro-pillars are chosen to be cuboids as describedat the beginning of Sec. 3 for the ease of the finite element implementation. The samplesizes here are defined to be the size of the cuboid base D . The (height to base size) aspectratio of the micro-pillars is L/D = 3. We choose D = 0 . µ m or 1 µ m.Sample Mean source Standard Max source τ Flow stresslength ( µ m) deviation ( µ m) length ( µ m) (MPa) (MPa)1 0.1835 0.1104 0.4832 170.7 420.02 0.2403 0.0935 0.4423 132.8 331.53 0.2039 0.0992 0.3340 145.5 357.94 0.1928 0.1083 0.3711 158.9 388.8 Table 1: Statistics of the individual sources whose lengths are obtained randomly following a uniformdistribution within [20nm, D ] for samples of size D = 0 . µ m. The locations of the sources are also determinedrandomly following the uniform distribution over the micro-pillar. Sample Mean source Standard Max source τ Flow stresslength ( µ m) deviation ( µ m) length ( µ m) (MPa) (MPa)1 0.4066 0.2252 0.9196 83.8 196.42 0.4047 0.2412 0.9995 91.0 230.23 0.4596 0.2333 0.9605 80.7 204.04 0.4709 0.2275 0.9349 80.3 201.2 Table 2: Same as Table 1 for samples of size D = 1 µ m. In our simulations, the lengths of the Frank-Read sources are generated randomly fol-lowing a uniform distribution within [20nm, D ]. The initial dislocation density ρ is also31andomly generated within the range 1 . ∼ × m − , and these pre-existing dislocationsegments of the sources are uniformly assigned to the twelve slip systems of fcc nickel. Thestatistics of the initial source distributions for various samples are listed in Tables 1 and 2.With these isolated Frank-Read sources generated, the corresponding source continuum isdeveloped following Eqs. (51)-(53). For the source character parameter C s in the criticalstress formula in Eq. (46), we choose it to be 1 .
35 and 2 .
02 for a single edge source and ascrew source, respectively, to agree with the simulation set-up in El-Awady et al. (2008).In Tables 1 and 2, we also show the values of τ , which is defined as the minimumvalue of the activation stress τ ( r ) of the Frank-Read source continuum given by Eq. (53),over the points where g ( r ) = 0 inside the pillar. The minimum activation stress τ will beused in the derivation of the scaling law in the next section. Engineering strain (%) E ng i nee r i ng s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (a) D = 0 . µ m Engineering strain (%) E ng i nee r i ng s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (b) D = 1 µ mFigure 7: Stress-strain curves obtained by the simulations of our continuum model. The vertical bars denotethe ranges of the flow stress predicted by the DDD simulations of El-Awady et al. (2008), whose averagedvalues are shown by the black dots. The stress-strain curves obtained by our continuum model are shown in Fig. 7 for samplesof sizes D = 0 . µ m and D = 1 µ m. Good agreement with the DDD results by El-Awadyet al. (2008) are observed. Firstly, both methods predict an initially elastic regime and analmost perfectly plastic regime, where work-hardening effect is barely observed. Secondly,both methods suggest that the engineering stress stays roughly unchanged or oscillate around32ome constant value in the regime of perfect plastic deformation, and this stress is measuredas the flow stress of the micro-pillars. The “smaller-being-stronger” size effect on crystallinestrength, which is indicated by the flow stress, is observed in both our simulations andin the DDD simulations of El-Awady et al. (2008). Moreover, statistical effects in theflow stress are seen in the simulations of both methods. Such statistical effects have alsobeen examined in other literature (e.g. El-Awady et al., 2009; Li et al., 2014). To makequantitative comparisons between the results of our model and those of the DDD simulationsof El-Awady et al. (2008), the respective ranges and averaged values for the flow stressrecorded by El-Awady et al. (2008) are drawn in Fig. 7. It can be seen that the flow stressescalculated by our continuum model all fall into the respective ranges predicted by their DDDsimulations.These comparisons show that our continuum model is able to provide an excellent sum-mary of the corresponding underlying dynamics of discrete dislocations. In the next section,we will use it to quantitatively study the size effect on strength observed in the uniaxialcompression tests of micro-pillars.
5. Size effect on strength of single-crystalline micro-pillars
Now we investigate the “smaller-being-stronger” size effect in micro-pillars using ourcontinuum model. For the initialization of the Frank-Read sources, we follow Shishvan andVan der Giessen (2010) that, in analogy to the distribution of grain sizes in polycrystalswhich has been experimentally measured, the Frank-Read source size l follows a log-normaldistribution with the probability density function p ( l ) = 1 √ πσ sd l e − (log l − log l m)2 √ σ (69)with two parameters l m and σ sd to be determined. The parameter l m which can be consideredas the effective mean source length should decrease with the pillar size D . The locationsof these source segments follow the uniform random distribution over the micro-pillar. The33ingle-ended sources (Parthasarathy et al., 2007) are included, and this happens when partof the source segments are outside the pillar.To validate the initial source distribution in Eq. (69) and determine the parameters, thenumerical results by using our continuum model are compared with the experimental dataof nickel micro-pillars by Dimiduk et al. (2005). In our simulations, most parameters arechosen the same as those used in Sec. 4.2 with the following exceptions in accordance withDimiduk et al. (2005): the loading axis is set along [269] and the singly active slip system isof the slip direction [101] and slip normal (¯111), the Schmid factor is m s = 0 .
48, the shearmodulus is 78GPa, and the aspect ratio of the pillar (
L/D ) is chosen randomly between 2and 3. The total density of all source segments is 3 × − m − following Dimiduk et al.(2005).To fit their experimental results, we choose the effective mean source length l m = α m D ,where α m = 1 /
15. The standard deviation σ sd is determined under the assumption that theprobability of a source segment that is longer than D is no more than 10 − , which gives σ sd = 0 .
4. The computed stress-strain curves with these values of parameters are shown inFig. 8 for different samples varying in size. It can be seen that the obtained flow stressesare in good agreement with the experimental data of Dimiduk et al. (2005). The size effecton strength is also clearly observed from our simulation results.
To describe how the pillar strength depends on its size D , we propose a formula σ flow = µb πm s α eff D log (cid:18) α eff Dr c (cid:19) + σ . (70)where σ is some constant stress, m s is recalled to be the Schmid factor, and α eff is adimensionless constant such that α eff D measures the length of the weakest source in pillarof size D . Eq. (70) suggests a scaling law σ flow ∼ bD log (cid:18) Db (cid:19) . (71)A comparison of Eq. (70) with experimental data (to be shown at the end of this subsection)suggests that Eq. (70) pocesses a wider effective zone, that is, it is valid for many (at least34 Shear strain (%) S hea r s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (a) D = 1 µ m Shear strain (%) S hea r s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (b) D = 2 . µ m Shear strain (%) S hea r s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (c) D = 5 µ m Shear strain (%) S hea r s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (d) D = 10 µ m Shear strain (%) S hea r s t r e ss ( M P a ) Sample 1Sample 2Sample 3Sample 4 (e) D = 20 µ mFigure 8: Stress-strain curves obtained by simulations using our continuum model with initial arrangementof Frank-Read sources whose length distribution following the log-normal one given by Eq. (69). Foursimulations are performed for each size of the micro-pillar D that ranges from 1 µ m to 20 µ m. The verticalbars identify the ranges of the experimentally measured flow stress by Dimiduk et al. (2005). τ (MPa) σ f l o w ( M P a ) µ m2.4 µ m5 µ m10 µ m20 µ m (a) log( α eff D/r c )/D τ m i n0 ( M P a ) µ m2.4 µ m5 µ m10 µ m20 µ m (b)Figure 9: (a) Linear dependence between the flow stress σ flow and the minimum activation stress of thesource continuum τ given in Eq. (72) (the dashed line), with fitted parameter σ . (b) The dependenceof the minimum activation stress τ on the sample size D given in Eq. (74) (the dashed line), where α eff = 1 / r c = 0 . b . These relations are validated by numerical results obtained using the continuummodel for pillars with different sizes (specified in the legend) shown by the dots in (a) and (b). We first show that the flow stress is related to the minimum activation stress τ by σ flow = τ m s + σ , (72)where σ is some constant stress, and recall that m s is the Schmid factor. This meansthat the flow stress is determined by the weakest Frank-Read source inside the pillar. Com-parisons of the results of this linear relation and the computed flow stresses by using thecontinuum model for pillars with different sizes are shown in Fig. 9(a), in which the pa-rameter σ is fitted from the numerical results. The excellent agreement seen in thesecomparisons demonstrates that Eq. (72) provides a good quantitative description for theflow stress σ flow in terms of τ .The next step is to relate τ to the sample size D . We follow the formula of the36ctivation stress of a single Frank-Read source to assume τ = C s µb πl eff log (cid:18) l eff r c (cid:19) , (73)where C s is dependent on the source character chosen to be 1, l eff can be considered as aneffective source length and it is assumed to be a fraction of D . Using simulation results ofthe continuum model for pillars with different sizes, fitting τ against l eff with reference toEq. (73) gives that l eff = α eff D , where the parameter α eff = 1 /
6. Physically, this means thatgiven a micro-pillar of size D , the longest source length is most likely α eff D . Hence τ isrelated to D by τ = µb πα eff D log (cid:18) α eff Dr c (cid:19) . (74)This relation is validated by numerical results obtained using the continuum model for pillarswith different sizes as shown in Fig. 9(b).Combining Eqs. (72) and (74), we obtain the scaling law in Eq. (70). Predictions ofEq. (70) are in good agreement with the experimental data as shown in Fig. 10 for singlecrystal nickel, aluminium and copper pillars. In each comparison in the figure, the resolvedflow stress τ flow = m s σ flow is plotted against the pillar size D , with α eff = 1 / r c = 0 . b and σ being a fitted parameter in the scaling law.It is noted that Eq. (74) has also been employed to predict the flow stress in polycrys-talline thin-film structures (von Blanckenhagen et al., 2001; Gruber et al., 2008). In theirworks, α eff is suggested to be 1 / / In this subsection, we discuss other plastic behaviors of micro-pillars that can be capturedby the continuum model.Firstly, recall that in our continuum model, the dislocation networks within the specimen(after local homogenization) are described by the DDPFs as defined in Sec. 2.2 and illustratedby Fig. 2, that is, the dislocation distribution on each slip plane is described by the contourcurves of the DDPF φ , and the distribution of the slip plane is represented by another DDPF ψ . Fig. 11 shows some snapshots of the dislocation substructures in some simulations using37 −5 −4 −3 −2 Sample size D ( µ m) S hea r s t r e ss τ f l o w ( µ ) Experiment Ni<269>log(D)/D scaling law (a) Nickel −1 −5 −4 −3 −2 Sample size D ( µ m) S hea r s t r e ss τ f l o w ( µ ) Experiment Al<123>Experiment Al<135>log(D)/D scaling law (b) Aluminium −1 −4 −3 −2 Sample size D ( µ m) S hea r s t r e ss τ f l o w ( µ ) Experiment Cu<111>log(D)/D scaling law (c) Copper < > −1 −4 −3 −2 Sample size D ( µ m) S hea r s t r e ss τ f l o w ( µ ) Experiment Cu<123>log(D)/D scaling law (d) Copper < > Figure 10: The log( D ) /D scaling law for the size-dependent effect of the flow stress of a micro-pillar given inEq. (70) resolved in the slip planes is examined by the experimental data for single crystal nickel (Dimiduket al., 2005), aluminium and copper (Uchic et al., 2009) micro-pillars. the continuum model for pillars of size D = 2 . µ m, 5 µ m, and 10 µ m. The distributions ofdislocation curves in the figure look smoothly varying. This is because the continuum modelonly resolves the dislocation microstructures in an average sense (see Sec. 2.2).We can also keep track of the two quantities that are commonly used to describe thestate of the specimen in a plastic deformation process: the total dislocation density ρ tot given by Eq. (44) and the plastic strain rate ˙ ǫ ptot given by Eq. (45). Numerical results of ρ tot and ˙ ǫ ptot obtained by simulations using our continuum model are plotted in Fig. 12.These results along with the stress-strain curves presented in Figs. 7 and 8 suggest that thefollowing evolution process may take place inside the micro-pillars when being compressed38 yx z (a) D = 2 . µ m −0.1 0 0.1 −0.100.1−0.500.5 yx z (b) D = 5 µ m −0.2 0 0.2 −0.200.2−0.500.5 yx z (c) D = 10 µ mFigure 11: Snapshots of the dislocation substructures (in an average sense as discussed in Sec. 2.2) in somesimulations using our continuum model for pillars of size (a) D = 2 . µ m, (b) 5 µ m, and (c) 10 µ m. Thesnapshots are taken from the perfect plastic deformation regime as shown by the stress-strain curves inFigs. 7 and 8. Engineering strain (%) ρ t o t ( µ m − ) µ m2.4 µ m5 µ m10 µ m20 µ m (a) Engineering strain (%) P l a s t i c s t r a i n r a t e ( s − ) µ m2.4 µ m5 µ m10 µ m20 µ m (b)Figure 12: Evolution of (a) the total dislocation density determined by Eq. (44) and (b) the plastic strainrate determined by Eq. (45) in micro-pillars with various sizes, obtained by using the continuum model. ǫ ptot converge to a same value for samples of various size.This is because this converged values are determined by the applied strain rates which arethe same for all samples in the compression tests. −0.20 0.2−0.200.2−0.500.5 xy z (a) 2.5% strain −0.20 0.2−0.200.2−0.500.5 xy z (b) 9% strain −0.20 0.2−0.200.2−0.500.5 xy z (c) 12% strainFigure 13: Shapes of a deformed micro-pillar of size D = 5 µ m under various applied strain: the cuboidsformed by the dashed-lines describe the original shapes of the pillar. By using the continuum model, we are also able to track the shape changes of themicro-pillars. Given u the displacement field, r + u is the position of a point, whose initialposition is at r . Examples of the profile of a deformed pillar during compression are shown inFig. 13. Note that in our simulations, spatial variation in the plastic shear slips induced bya distribution of Frank-Read sources has been locally averaged over many discrete sources inthe source continuum formulation. The pillar profiles in Fig. 13 agree in the averaged sensewith the experimental observations (e.g. Fig. 4 in Dimiduk et al. (2005)). With smaller40umber of sources, relatively strongly non-uniform displacement fields along z direction canbe observed in our simulations, see Fig. 14 for an example, which more accurately agreewith the slip-band structures observed in the experiments. For an extreme case, if we checkthe surface of a micro-pillar containing only one operating Frank-Read source as shown inFig. 6, localized displacement gradient is clearly observed around the activated slip plane.These agree with the conclusion of Akarapu et al. (2010) drawn from simulations using ahybrid elasto-viscoplastic model which couples DDD that such localized deformation, whichis strong for small-size pillars, is due to the heterogeneous dislocation distributions that leadto clusters of sources. More efficient numerical implementation method of our continuummodel with fine meshes will help resolve more accurately the pillar profile observed in theexperiments, which will be developed in the future work. −0.5 0 0.52.53.54.55.5 x 10 −4 z u (a) u −0.5 0 0.5−404x 10 −4 z u (b) u −0.5 0 0.5−3−2−10 x 10 −3 z u (c) u Figure 14: Displacement field u = ( u , u , u ) along a vertical edge of a micro-pillar with size D = 1 µ m and8 Frank-Read sources under 0 .
3% applied strain. The length unit is L .
6. Conclusion
In this paper, we have presented a dislocation-density-based continuum model to studythe plastic behavior of crystals, in which the dislocation substructures are represented bypairs of DDPFs. A discrete dislocation network in three dimensions is approximated bya dislocation continuum after local homogenization of dislocation ensembles within somerepresentative volume. For the λ -th slip plane in the dislocation continuum, a DDPF φ λ isdefined such that φ λ restricted on each slip plane describes the plastic slip across the slip41lane and identifies the distribution of dislocation curves by its contour lines, and anotherDDPF ψ λ is employed to describe the slip plane distribution by its contour surfaces. Themajor advantage of this three-dimensional continuum model lies in its simple representa-tion of distributions of curved dislocations. The connectivity condition of dislocations isautomatically satisfied with this representation of dislocations.Based on the DDPFs, the plastic deformation process of crystals can be formulatedby a system of equations as given by Eqs. (54)–(57). We have shown that the equationsystem provides an effective summary over the underlying discrete dislocation dynamics,including the long-range elastic interaction of dislocations, dislocation line tension effect, anddislocation multiplication by Frank-Read sources. Numerically, a finite element formulationis proposed to compute the long-range stress field.The continuum model is validated by comparisons with DDD simulations and experi-mental data. As one application of the continuum model characterized by DDPFs, the sizeeffect on the strength of micro-pillars is studied, and the pillar flow stress is found scalingwith its (non-dimensionalized) pillar size D by log( D ) /D .Future work may include generalizations of the continuum model characterized by DDPFsto incorporate the out-of-slip-plane dislocation motion (cross-slip or climb) in which the con-tour surfaces of the DDPF ψ may be curved and evolved in the plastic deformation. Wewill also take into account in the continuum model other short-range dislocation interactionssuch as dislocation reactions and junction formation, as well as dislocation interactions withother defects (e.g. Xiang and Srolovitz, 2006; Chen et al., 2010; Zhu and Xiang, 2014).Efficient numerical implementation method will also be considered in the future work. Acknowledgement
This work was partially supported by the Hong Kong Research Grants Council GeneralResearch Fund 606313. 42 eferences
A. Acharya. A model of crystal plasticity based on the theory of continuously distributed dislocations.
J.Mech. Phys. Solids , 49:761–784, 2001.S. Akarapu, H.M. Zbib, and D.F. Bahr. Analysis of heterogeneous deformation and dislocation dynamics insingle crystal micropillars under compression.
Int. J. Plast. , 26:239–257, 2010.A. Alankar, P. Eisenlohr, and D. Raabe. A dislocation density-based crystal plasticity constitutive modelfor prismatic slip in -titanium.
Acta Mater. , 59:7003–7009, 2011.A. Arsenlis and D. M. Parks. Modeling the evolution of crystallographic dislocation density in crystalplasticity.
J. Mech. Phys. Solids , 50:1979–2009, 2002.A. Arsenlis, W. Cai, M. Tang, M. Rhee, T. Oppelstrup, T. G. Hommes, T. G. Pierce, and V. V. Bulatov.Enabling strain hardening simulations with dislocation dynamics.
Modelling Simul. Mater. Sci. Eng. , 15:553–595, 2007.A. A. Benzerga, Y. Br´echet, A. Neeldeman, and E. Van der Giessen. Incorporating three-dimensionalmechanisms into two-dimensional dislocation dynamics.
Modelling Simul. Mater. Sci. Eng. , 12:159–196,2004.V. L. Berdichevsky. On thermodynamics of crystal plasticity.
Scripta Mater. , 54:711–716, 2006.Z. Chen, K. T. Chu, D. J. Srolovitz, J. M. Rickman, and M. P. Haataja. Dislocation climb strengthening insystems with immobile obstacles: Three-dimensional level-set simulation study.
Phys. Rev. B , 81:054104,2010.B. Cheng, H. S. Leung, and A. H. W. Ngan. Strength of metals under vibrations - dislocation-density-function dynamics simulations.
Philos. Mag. , pages 1–21, 2014.D. M. Dimiduk, M. D. Uchic, and T. A. Parthasarathy. Size-affected single-slip behavior of pure nickelmicrocrystals.
Acta Mater. , 53:4065–4077, 2005.J. A. El-Awady, S. B. Biner, and N. M. Ghoniem. A self-consistent boundary element, parametric dislocationdynamics formulation of plastic flow in finite volumes.
J. Mech. Phys. Solids , 56:2019–2035, 2008.J. A. El-Awady, M. Wen, and N. M. Ghoniem. The role of the weakest link mechanism in controlling theplasticity of micropillars.
J. Mech. Phys. Solids , 57:32–50, 2009.A. El-Azab. Statistical mechanics treatment of the evolution of dislocation distributions in single crystals.
Phys. Rev. B , 61:11956–11966, 2000.P. Engels, A. X. Ma, and A. Hartmaier. Continuum simulation of the evolution of dislocation densitiesduring nanoindentation.
Int. J. Plast. , 38:159–169, 2012.M. Fivel, L. Tabourot, E. Rauch, and G. Canova. Identification through mesoscopic simulations of macro-scopic parameters of physically based constitutive equations for the plastic behaviour of fcc single crystals.
J. Phys. IV France , 8:151–158, 1998. . A. Fleck and J. W. Hutchinson. A phenomenological theory for strain gradient effects in plasticity. J.Mech. Phys. Solids , 41:1825–1857, 1993.M. G. D. Geers, R. H. H. Peerlings, M. A. Peletier, and L. Scardia. Asymptotic behaviour of a pile-up ofinfinite walls of edge dislocations.
Arch. Ration. Mech. Anal. , 209:495–539, 2013.N. M. Ghoniem, S. H. Tong, and L. Z. Sun. Parametric dislocation dynamics: a thermodynamics-basedapproach to investgations of mesoscopic plastic deformation.
Phys. Rev. B , 61:913–927, 2000.J. R. Greer and W. D. Nix. Nanoscale gold pillars strengthened through dislocation starvation.
Phys. Rev.B , 73, 2006.J. R. Greer, W. C. Oliver, and W. D. Nix. Size dependence of mechanical properties of gold at the micronscale in the absence of strain gradients.
Acta Mater. , 53:1821–1830, 2005.I. Groma, F. F. Csikor, and M. Zaiser. Spatial correlations and higher-order gradient terms in a continuumdescription of dislocation dynamics.
Acta Mater. , 51:1271–1281, 2003.P. A. Gruber, J. Bohm, F. Onuseit, A. Wanner, R. Spolenak, and E. Arzt. Size effects on yield strengthand strain hardening for ultra-thin cu films with and without passivation: A study by synchrotron andbulge test techniques.
Acta Mater. , 56(10):2318–2335, 2008.R. Gu and A. H. W. Ngan. Dislocation arrangement in small crystal volumes determines power-law sizedependence of yield strength.
J. Mech. Phys. Solids , 61:1531–1542, 2013.M. E. Gurtin. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessarydislocations.
J. Mech. Phys. Solids , 50:5–32, 2002.C. L. Hall. Asymptotic analysis of a pile-up of regular edge dislocation walls.
Mater. Sci. Eng. A , 530:144–148, 2011.A. K. Head, S. D. Howison, J. R. Ockendon, and S. P. Tighe. An equilibrium-theory of dislocation continua.
SIAM Rev. , 35:580–609, 1993.J. P. Hirth and J. Lothe.
Theory of dislocations . Wiley, New York, 2nd edition, 1982.T. Hochrainer, S. Sandfeld, M. Zaiser, and P. Gumbsch. Continuum dislocation dynamics: Towards aphysical theory of crystal plasticity.
J. Mech. Phys. Solids , 63:167–178, 2014.D. C. Jang, X. Y. Li, H. J. Gao, and J. R. Greer. Deformation mechanisms in nanotwinned metal nanopillars.
Nature Nanotechnology , 7:594–601, 2012.D. M. Kochmann and K. C. Le. Dislocation pile-ups in bicrystals within continuum dislocation theory.
Int.J. Plast. , 24:2125–2147, 2008.A. Kosevich. Crystal dislocations and the theory of elasticity. In
Dislocations in Solids, Vol. 1 , pages 33–141.North-Holland, Amsterdam, 1979.E. Kroener. Dislocation: a new concept in the continuum theory of plasticity.
J. Math. Phys. , 42:27–37,1963. . P. Kubin, G. Canova, M. Condat, B. Devincre, V. Pontikis, and Y. Br´echet. Dislocation microstructuresand plastic flow: a 3d simulation. Solid State Phenom. , 23/24:455–472, 1992.K. C. Le and C. Guenther. Nonlinear continuum dislocation theory revisited.
Int. J. Plast. , 53:164–178,2014.K. C. Le and C. Guenther. Martensitic phase transition involving dislocations.
J. Mech. Phys. Solids , 79:67–79, 2015.D. S. Li, H. M. Zbib, X. Sun, and M. Khaleel. Predicting plastic flow and irradiation hardening of ironsingle crystal with mechanism-based continuum dislocation dyanmics.
Int. J. Plast. , 52:3–17, 2014.T. Mura.
Micromechanics of Defects in Solids . Dordrecht: Martinus Nijhoff, 1987.F. R. N. Nabarro. Dislocations in a simple cubic lattice.
Proc. Phys. Soc. , 59:256–272, 1947.D. Nelson and J. Toner. Bond-orientational order, dislocation loops, and melting of solids and smectic-aliquid crystals.
Phys. Rev. B , 24:363–387, 1981.W. D. Nix and H. Gao. Indentation size effects in crystalline materials: a law for strain gradient plasticity.
J. Mech. Phys. Solids , 46:411–425, 1998.J. F. Nye. Some geometrical relations in dislocated crystals.
Acta Metall. , 1:153–162, 1953.M. S. Oztop, C. F. Niordson, and J. W. Kysar. Length-scale effect due to periodic variation of geometricallynecessary dislocation densities.
Int. J. Plast. , 41:189–201, 2013.T. A. Parthasarathy, S. I. Rao, D. M. Dimiduk, M. D. Uchic, and D. R. Trinkle. Contribution to size effectof yield strength from the stochastics of dislocation source lengths in finite samples.
Scripta Mater. , 56:313–316, 2007.R. Peierls. The size of a dislocation.
Proc. Phys. Soc. , 52:34–37, 1940.D. Peirce, R. J. Asaro, and A. Needleman. Material rate dependence and localized deformation in crystallinesolids.
Acta metall. , 31:1951–1976, 1983.S. S. Quek, Y. Xiang, Y. W. Zhang, D. J. Srolovitz, and C. Lu. Level set simulation of dislocation dynamicsin thin films.
Acta Mater. , 54:2371–2381, 2006.S. I. Rao, D. M. Dimiduk, M. Tang, T. A. Parthasarathy, M. D. Uchic, and C. Woodward. Estimating thestrength of single-ended dislocation sources in micron-sized single crystals.
Philos. Mag. , 87:4777–4794,2007.J. R. Rice. Inelastic constitutive relations for solids: An internal-variable theory and its application to metalplasticity.
J. Mech. Phys. Solids , 19:433–455, 1971.I. Ryu, W. D. Nix, and W. Cai. Plasticity of bcc micropillars controlled by competition between dislocationmultiplication and depletion.
Acta Mater. , 61:3233–3241, 2013.S. Sandfeld, T. Hochrainer, M. Zaiser, and P. Gumbsch. Continuum modeling of dislocation plasticity:Theory, numerical implementation, and validation by discrete dislocation simulations.
J. Mater. Res. , 26: Modelling Simul. Mater. Sci. Eng. , 22, 2014.R. Sedl´aˇcek, J. Kratochv´ıl, and E. Werner. The importance of being curved: bowing dislocations in acontinuum description.
Philos. Mag. , 83:3735–3752, 2003.J. Senger, D. Weygand, P. Gumbsch, and O. Kraft. Discrete dislocation simulations of the plasticity ofmicro-pillars under uniaxial loading.
Scripta Mater. , 58:587–590, 2008.S. Shao, N. Abdolrahim, D. F. Bahr, G. Lin, and H. M. Zbib. Stochastic effects in plasticity in smallvolumes.
Int. J. Plast. , 52:117–132, 2014.S. S. Shishvan and E. Van der Giessen. Distribution of dislocation source length and the size dependentyield strength in freestanding thin films.
J. Mech. Phys. Solids , 58:678–695, 2010.B. Svendsen. Continuum thermodynamic models for crystal plasticity including the effects of geometrically-necessary dislocations.
J. Mech. Phys. Solids , 50:1297–1329, 2002.H. Tang, K. W. Schwarz, and H. D. Espinosa. Dislocation-source shutdown and the plastic behavior ofsingle-crystal micropillars.
Phys. Rev. Lett. , 100, 2008.M. D. Uchic, D. M. Dimiduk, J. N. Florando, and W. D. Nix. Sample dimensions influence strength andcrystal plasticity.
Science , 305:986–989, 2004.M. D. Uchic, P. A. Shade, and D. M. Dimiduk. Plasticity of micrometer-scale single crystals in compression.
Annu. Rev. Mater. Res. , 39:361 – 386, 2009.E. Van der Giessen and A. Needleman. Discrete dislocation plasticity - a simple planar model.
ModellingSimul. Mater. Sci. Eng. , 3:689–735, 1995.B. von Blanckenhagen, P. Gumbsch, and E. Arzt. Dislocation sources in discrete dislocation simulations ofthin-film plasticity and the hall-petch relation.
Modelling Simul. Mater. Sci. Eng. , 9:157–169, 2001.R. E. Voskoboinikov, S. J. Chapman, J. R. Ockendon, and D. J. Allwright. Continuum and discrete modelsof dislocation pile-ups. i. pile-up at a lock.
J. Mech. Phys. Solids , 55:2007–2025, 2007.J. Wang, R.F. Zhang, C.Z. Zhou, I.J. Beyerlein, and A. Misra. Interface dislocation patterns and dislocationnucleation in face-centered-cubic and body-centered-cubic bicrystal interfaces.
Int. J. Plast. , 53:40–55,2014.D. Weygand, L. H. Friedman, E. Van der Giessen, and A. Needleman. Aspects of boundary-value problemsolutions with three-dimensional dislocation dynamics.
Modelling Simul. Mater. Sci. Eng. , 10:437–468,2002.Y. Xiang. Continuum approximation of the peach-koehler force on dislocations in a slip plane.
J. Mech.Phys. Solids , 57:728–743, 2009.Y. Xiang and D. J. Srolovitz. Dislocation climb effects on particle bypass mechanisms.
Philos. Mag. , 86: ActaMater. , 51:5499–5518, 2003.Y. Xiang, H. Wei, P. B. Ming, and W. E. A generalized peierls-nabarro model for curved dislocations andcore structures of dislocation loops in al and cu.
Acta Mater. , 56:1447–1460, 2008.H. M. Zbib, M. Rhee, and J.P. Hirth. On plastic deformation and the dynamics of 3d dislocations.
Int. J.Mech. Sci. , 40:113–127, 1998.D. G. Zhao, H. Q. Wang, and Y. Xiang. Asymptotic behaviors of the stress fields in the vicinity of dislocationsand dislocation segments.
Philos. Mag. , 92:2351–2374, 2012.C. Z. Zhou and R. LeSar. Dislocation dynamics simulations of plasticity in polycrystalline thin films.
Int.J. Plast. , 30-31:185–201, 2012.X. H. Zhu and Y. Xiang. Continuum model for dislocation dynamics in a slip plane.
Philos. Mag. , 90:4409–4428, 2010.X. H. Zhu and Y. Xiang. Continuum framework for dislocation structure, energy and dynamics of dislocationarrays and low angle grain boundaries.
J. Mech. Phys. Solids , 69:175–194, 2014.Y. C. Zhu and S. J. Chapman. A natural transition between equilibrium patterns of dislocation dipoles.
J.Elast. , 117:51–61, 2014a.Y. C. Zhu and S. J. Chapman. Motion of screw segments in the early stage of fatigue testing.
Mater. Sci.Eng. A , 589:132–139, 2014b.Y. C. Zhu, S. J. Chapman, and A. Acharya. Dislocation motion and instability.
J. Mech. Phys. Solids , 61:1835–1853, 2013.Y. C. Zhu, H. Q. Wang, X. H. Zhu, and Y. Xiang. A continuum model for dislocation dynamics incorporatingfrank-read sources and hall-petch relation in two dimensions.
Int. J. Plast. , 60:19–39, 2014., 60:19–39, 2014.