Computational Modeling of Dislocation Slip Mechanisms in Crystal Plasticity: A Short Review

The bridge between classical continuum plasticity and crystal plasticity is becoming narrower with continuously improved computational power and with engineers’ desire to obtain more information and better accuracy from their simulations, incorporating at the same time more effects about the microstructure of the material. This paper presents a short overview of the main current techniques employed in crystal plasticity formulations for finite element analysis, as to serve as a point of departure for researchers willing to incorporate microstructure effects in elastoplastic simulations. We include both classical and novel crystal plasticity formulations, as well as the different approaches to model dislocations in crystals.


Introduction
Analysis of the plastic deformations of metals is very important in many aspects of the design of metallic goods, for example, in the design of the alloys which will be used in the components, in their manufacturing process (e.g., shape forming), and in the behaviour of the material in service, both to sustain the ordinary actions and to behave in a controlled, predictable way (fracture, fatigue, crash-worthiness, etc.) [1,2].
Plastic slip is the most common plastic deformation mechanism in crystalline solids, e.g., most of the metals. This slip takes place as a plastic shearing in specific planes and directions governed by the arrangement of the atoms in a regular crystal structure [3]. This arrangement produces anisotropic elastic and plastic properties depending on the packing of atoms into planes. Compact planes slide in particular directions, a movement which is facilitated by the presence of dislocations (defects in the crystal arrangement, e.g., by the absence of atoms in the structure, line defects). Plastic slips in crystals is a plastic process that does not change the volume. However, for materials whose dominant deformation mechanism is slip, it is usually the origin of yield stress and hardening in macroscopic plasticity. Then, the study of crystal plasticity is fundamental to understand in a rational way the macroscopic behavior of the polycrystal, and how different aspects of the microstructure affect this behavior and the observed critical quantities like yield stress, elastoplastic anisotropy, anisotropic local and global deformation, hardening, thermal properties, etc. [4]. In the late 20th and beginning of the the 21st Century, after the work of Pierce et al. [5], thanks to improved computational power, crystal plasticity finite element simulations have received an important boost because these multiscale procedures, in which the studied phenomena has a length scale of µm, can now be addressed computationally. Then, computational homogenization schemes [6] are a typical procedure to bring the microscale complexity to the macroscale in a simple, standarized way. Thus, current computational power increasingly facilitates this type of analysis, bringing classical materials science to the engineering office [7].
Whereas the scale at which dislocation movements takes place is much smaller than the continuum scale (or precisely because of that, which facilitates the requirement of separation of scales), many effects may be studied at a mesoscale if we consider that a given stress integration point represents a homogenization of the grain, or part of the grain, so crystal plasticity computations give an implicit Representative Volume Element (RVE) of the dislocation theory, being the direction of the Burgers vector the sliding direction, and its modulus the computed sliding times the length of the RVE [8,9].
The choice of crystal plasticity theory as the simulation tool for different inelastic phenomena depends itself on the objectives of such simulation. Several aspects of the behavior of materials are easily extracted from this type of simulations. Examples are texture evolution (e.g., anisotropy symmetries evolution as a result of the evolution of crystallographic planes and grain structure), effects of grain size, effects of twinning or impurities, speed of deformation, influence of grain shape (e.g., after cold forming), etc. Furthermore, damage in these materials due to void nucleation in grain boundaries has also been studied through crystal plasticity simulations (e.g., [10,11], see also [12]). The coupling between the crystal plasticity and phase-field damage is also proposed to study the fracture [13,14]. There are other options, mainly computational, to model plasticity-based phenomena at different scales, specially through the combination of one scale with a larger one in which the desired effect needs to be accounted for. Examples are atomistic models [15,16], dislocation dynamics models [17], and continuum models [18,19]. However, crystal plasticity is probably the most used framework for incorporating the grain and crystal structure as well as the inherent effects of dislocations in the continuum modelling. Then, in this paper after a brief review of the physical aspects, we perform a review of the main ingredients and formulations employed today in crystal plasticity simulations. We finish with some demonstrative simulations as examples of the possibilities of crystal plasticity. The purpose of the review is to serve as an introductory compact presentation of current mostly used formulations for researchers (e.g., continuum mechanics ones) becoming involved in the subject or willing to incorporate texture effects in advanced simulations, for example, by computational homogenization. It is by no means a complete review, so many important works and formulation aspects, are omitted. Advanced long reviews taylored for researchers in materials science are available elsewhere [20][21][22].
Dislocations are crystal defects which result, for example, from the absence of some atoms in the otherwise regular arrangement, see Figure 2c,d. Microscopically, slip is a result of the glide of many dislocations, facilitated by the stress field generated around them. This glide does not happen simultaneously in all the grain, but occurs in a worm-like movement, see Figure 3. The glide motion of a dislocation means that the dislocation moves in the plane which contains both its line and the Burgers vector b (Figure 2c). Dislocations and their associated local stress fields interact, blocking or further facilitating their movement. The accumulation of hundreds or even thousands of these displacements form a slip band. These slip bands can be easily observed on the surfaces of polished specimens which have been subjected to plastic deformation. Dislocations, slip bands and slip displacements are small with respect to the size of the grains, so even when the grain microstructure is considered, this plastic slip can be treated as homogeneous from a macroscopic standpoint [25]. The dislocation mechanisms can be homogenized in the form of statistical dislocation populations, and therefore, can be embedded in state evolution constitutive equations in a crystal plasticity theory. For a certain crystal, a slip system is determined by the normal direction of the slip plane m g (which we can call the plane "direction") and the slip direction s g in the crystal coordinate system. In Figure 2a,b, the plane formed by red atoms is a typical slip plane in a FCC crystal. At the atomic level, the forming of slip requires some breaking and reforming of dislocation-stressed interatomic bonds. This can be interpreted as a plane of atoms sliding over another parallel one on the slip plane, with the sliding direction as the slip direction. This is a volume-preserving shearing process that is typically not affected by the hydrostatic stress [26]. Slip planes m g are planes with the highest packing or atom density, whereas the slip directions s g are usually directions in these slip planes in which the atoms have the closest packing, see again Figure 2. The dyad s g ⊗ m g constitutes the slip system andγ g is the slip (shearing) rate. Sliding may occur in the positive +s g and negative −s g directions (hence a FCC crystal may have 4 × 3 = 12 slip systems).
Other important quantities of a slip system are the resolved shear stress (RSS) and the critical resolved shear stress (CRSS). The former is the shear stress in the plane and direction of the slip system, and the latter is the critical value for which sliding occurs. For most crystals, the CRSS values of a certain slip system for both "+" and "−" directions can be assumed to be the same [25].
The CRSS of a slip system depends on many factors, among them previous slips and the crystal structure itself, because they play an important role due to the difficulty in breaking and reforming interatomic bonds. Also, an increase in the temperature (increasing vibration) or decrease in the strain rate (facilitating adaptation) increase the thermal activation of atoms and allow dislocations to overcome obstacles more easily. This results in an effective lower CRSS [27][28][29]. Of course, the modification of the material composition by adding other elements increase in general obstacles (breaking crystal regularity) and could interrupt the glide of the dislocations as the case of BCC iron micropillar [30]. Dislocations pile-up occurs as a cluster of dislocations which are unable to move past the grain boundary. Each successive dislocation will apply a repulsive force to the dislocation incident with the grain boundary. These repulsive forces act as a driving force to help overcoming the energetic barrier for diffusion across the boundary, such that additional pile up causes dislocation diffusion across the grain boundary, allowing further plastic deformation in the material. Decreasing grain size decreases the statistical amount of possible pile up at the boundary, increasing the amount of applied stress necessary to move a dislocation across a grain boundary. The higher the applied stress needed to move the dislocation, the higher the macroscopic yield stress [27]. This effect is well known as the Hall-Petch effect/law [31,32]. Of course, given a grain size, these effects deppend on the amount of statistical dislocations and other deffects.
In essence, slip plane m g , slip direction s g and the CRSS τ c g are crucial features of a slip system. According to the Schmid law, they determine which slip systems are "active" (sliding) [33]. The crystal structure is anisotropic, and the elasticity and plastic slip are also anisotropic. Furthermore, different slip systems can have different CRSS, which further determine the anisotropy of the material. This holds especially if the material is a single crystal or a polycrystal with a strong texture. Only when the material is a polycrystal that has a large number of different crystallographic orientations, the material can be considered as statistically isotropic on average; and von Mises plasticity is adequate. If the crystallographic orientation has a defined statistical distribution, continuum anisotropic models using Hill's or Barlat's yield functions with anisotropic or isotropic elasticity may also be adequate in some cases. However, note that in continuum plasticity tracking the macroscopic evolution of such anisotropy is not simple [34], which is a good reason to resort to crystal plasticity.
Rate-dependent crystal plasticity modeling is all about computing the shear amount of the activated slip systems from which the overall plastic deformation follows. The key point becomes how to describe/model the evolution of the CRSS of slips. The "hardening" evolution of the CRSS of slips is caused microscopically by the overcoming obstacles during slip gliding and from the interaction with other dislocations, as well as from the generation and annihilation of dislocations. As mentioned, temperature and strain rate are related to the thermal activation, so they are also related to the easiness for a dislocation to glide [27]. Different types of obstacles (as grain boundaries, solutes, and precipitates or other phases) may hinder gliding [30]. Of course, some factors changing the dislocation density or which may influence the interactions among dislocations may contribute to work hardening. An important one is the already pre-existing dislocation density. As well known, this could be reduced by heat treatments like annealing. Active secondary slip systems can cause multiplication of dislocations. Slip transfer at grain boundaries and cross slip can cause annihilation of dislocations. Stress concentrations (mainly at grain boundaries) can even cause the nucleation of dislocations. In practice, the different crossinteractions between dislocations which increase the CRSS are accounted for by means of the latent hardening [35][36][37].
As explained below, to prescribe the evolution of the CRSS of the slip system, different models have been proposed in the literature by relating its value and evolution to different variables during deformation. One widely used category follows a "phenomenological" approach [5,[38][39][40]. They relate the evolution of the CRSS with the shear rate of all the activated slip systems. Another popular type are the "physical" models, which are dislocation-density based. This type of model relates the evolution of the CRSS with the evolution of the dislocation density of all the slip systems, e.g., the Taylor law [41], the early models [42,43]. This approach is extended to account for more detailed factors that influence hardening, e.g., the annihilation of dislocations [44]; the latent hardening [45][46][47]; the grain size effect [48]; the slip transfer at the grain boundary [49]; the influence of strain gradient [50]; the incorporation of phase transformation [51].
Dislocation slip mechanisms are not the only ones present in crystals producing permanent deformations. Aside, there are other types of mechanisms as displacive (difusionless) transformations and twinning. The former is a rearrangement of the atoms of a crystal structure to form another structure without a relevant change in volume. The most important case appears in carbon steels (e.g., martensitic transformations). Another plastic mechanism is twinning deformation, which produces mirrored structures. The treatment of these additional mechanisms are out of the scope of this review.

Kinematics
As mentioned, the structure of a crystalline solid at the µm scale may be understood, and hence it is often modelled, as an aggregate of volume elements (finite elements), or lattice blocks (voxels). This constitutes the reference "undeformed" lattice, which contains dislocations (hence microscopically deformed); Figure 2. At the continuum scale (e.g., at a stress integration point in a finite element, computed from the reference coordinates X, the updated coordinates x = X + u, and displacement field u as F = ∂x/∂X = I + ∇ X u with ∇ X being the reference-coordinates gradient), the deformation can be considered as composition of two basic mechanisms that take place simultaneously: (1) slip between lattice blocks (as an homogenization of dislocation slips), followed by (2) elastic lattice distortions which include (3) lattice rotations (local or "elastic" and rigid-body) that also arise under a general deformation gradient F. Rotations produce observable phenomena of geometric origin, such as the geometric softening of crystals due to their reorientation [52] (similar to structural P-δ effects). In consequence, at large deformations, Kröner and Lee [53] introduced a multiplicative decomposition of F into plastic slips and crystal "elastic" deformations [54,55]. Whereas F is a compatible deformation gradient (X-integrable from the potential u(X) or F −1 is x−integrable from u(x)), it is decomposed into two (local, incompatible) deformation components as: where the plastic deformation gradient F p is irreversible and transforms the reference configuration into a fictitious local state, termed intermediate configuration. This ("elastically unloaded") intermediate configuration is "stresses-free" and is considered to maintain an "undeformed" slipped lattice. In turn, the elastic deformation gradient F e is the deformation considered reversible, and includes rigid body rotations; namely after an overimposed rotation Q we have F + := QF = QF e F p = F + e F p . The possible indeterminacy of the intermediate configuration respect to arbitrary rotations of the type F e Q T QF p has involved several discussions (see e.g., [56][57][58]). In summary, specially in the case of crystal plasticity, F p is continuously time-integrated, so it is always uniquely defined. Remarkably, the conservative/dissipative character of the two gradients give an important physical difference in their polar decompositions: whereas F e = R e U e is meaningful in the sense that R e are path-independent rotations which in principle do not enter the constitutive equations and U e are path-independent stretches that may be used in a stored energy as state variables, F p = R p U p is physically meaningless because the crystallographic slip does not change the orientation of the crystal lattice even though F p include mathematically a rotational part; stretch and rotations are related by the isoclinic flow and must be considered in an incremental setting. Indeed, this is an argument to build the formulation in terms of elastic measurements, avoiding also weak-invariance issues [59,60].
As mentioned, whereas F is compatible, the gradients F e and F p are incompatible, meaning that for general non-homogeneous deformations we havē where dξ is the transformation in the intermediate configuration of the material differential dX, and Γ is the transformation of a close circuit C in the reference configuration, i.e., C dX = 0. Even though the gradients are to be interpreted in a continuum sense (at a scale much larger than that of the dislocations),b represents the resultant Burgers vector (net Burgers vector of the dislocations tansversing Γ); it may be conceptually interpreted as the Burgers vector b (the sign may be reversed) when one dislocation is present; i.e.,b has a modulus proportional to the amount (density) of dislocations; see an ilustration of the concept in Figure 4 motivated in the typical Burgers circuit, and a discussion of the continuum/crystallographic views in [61]. To the last integral in Equation (2), the generalized Stokes theorem may be applied where S Γ is a surface with support on Γ, dl is a line differential of Γ, and ds is a surface differential of S Γ . The operation * is any operation consistent with the physical object A (vector or tensor of any order), and ∇ indicates that the operator operates on A. Then where the last identity is inmmediate using index notation. In dislocation density based models, the Kröner-Nye tensor α is handy [62,63], usually defined as the tensor field α(X) Of course, because the deformation gradient F is compatible, so c F e F p · dX = c dx = 0, the same integral may be obtained from the spatial configuration and For the case of infinitessimal deformations, splitting the local displacement into elastic and plastic local (incompatible) parts (u = e e + u p ), we immediately obtain  The multiplicative decomposition is accepted in crystal plasticity, but at a single integration point in a finite element program, the decomposition assumes that the represented volume is fully transversed by dislocations. However, in reality, this assumption is not fulfilled; in general, there is still a dislocation content between neighboring material points if both experience different slip system activity from nonhomogeneous deformations, even if a full traverse by dislocations or a (statistical) homogeneous population is considered for individual material points. Much research is ongoing in the seek of the best formulation to take into account these effects [64][65][66][67][68][69][70][71][72][73].
The constitutive law can be obtained from the Helmholtz free energy density expressed in the intermediate configuration as proposed by [8,69,74]. The free energy density can be decomposed into the elastic and plastic parts [8]. The elastic part is recoverable through elastic strain unloading, whereas the plastic part needs some plastic deformation to recover it (accounts for example, for kinematic hardening).
where q is a set of internal variables. The second Piola-Kirchhoff stress tensor in the intermediate configuration then is obtained, for example, using Green-Lagrange strains, as being E e the Green-Lagrange elastic strain tensor and S |e its work-conjugate stress tensor (both living in the intermediate configuration). Note that the strain energy may be written in terms of any arbitrary deformation measure, and the derivative of the energy respect to that measure gives the work-conjugate stress measure in the same configuration. If the elastic strain energy Ψ e is defined by a quadratic potential of the Green-Lagrange elastic strains, S |e can be obtained as a linear function of E e , where C is the fourth-order elastic tangent moduli tensor of the single crystal.
Considering the decomposition in Equation (1), the velocity gradient L can similarly be decomposed: being L e and L p the elastic and plastic velocity gradients. It is noted that the total deformation gradient is decomposed multiplicatively, whereas the decomposition of total velocity gradient is additive. Furthermore, both L and L e are defined with respect to the current configuration. The plastic velocity gradient L p , however, lives in the intermediate configuration, so it requires a mapping (push forward) from the intermediate to the current configuration, as shown by Equation (9).
Rice [54] assumed that for multislips, the plastic velocity gradient L p is the sum of all crystallographic slip rates: where γ g is the internal variable that accounts for the accumulated plastic slip in each glide (slip) system g; vectors m g and s g are, respectively, unit normal vector to the slip plane and the unit vector in the slip direction that define the same glide system g in the reference configuration. An important aspect here is that the addition in Equation (10) is not exact in discrete algorithms, because the order in which the different slips take place changes the final result. The reason is that after one slip has taken place, L p for the next one would be in a different configuration; see Figure 5. However, the approximation is considered acceptable in the literature given that steps are typically small. Indeed, recently we have proposed a different framework in terms of corrector rates [75], which results immediately from the chain rule applied to F e (F, where ctḞ e is the continuum corrector rate of the elastic gradient, and trḞ e is the continuum trial rate of the elastic gradient; see further details about the elastic corrector rates framework in [75][76][77][78]. In this case, Equation (9) may be written as where tr L e is the elastic rate with the plastic flow frozen (the concept of a partial derivative), and ct L e is the elastic rate with the external power frozen. Then, as an alternative to Equation (10), the following addition takes place in the same trial-frozen configuration: As a key difference, this poses an evolution on the elastic deformation gradient instead of on the plastic one, facilitating by construction the weak-invariance property.

Phenomenological Models
In crystal plasticity, the phenomenological models can be considered as the first approach adapted from the framework of continuum mechanics, and are dominant in the literature. This approach uses mostly the resolved shear stress (RSS) τ g on a glide system g, as material state variable. With S |e denoting the second Piola-Kirchhoff stress tensor in the intermediate configuration, the RSS on glide system g, is typically defined as: where F T e F e S |e is the non-symmetric Mandel stress tensor. In the case of infinitesimal elastic strains which typically hold in crystal plasticity, Equation (14) can be approximated as: An alternative is to use the recent formulation based of elastic corrector rates, which do not involve the Mandel stress tensor, but the symmetric generalized Kirchhoff one [75,[79][80][81]. This formulation has general validity (not only for infinitesimal elastic strains). In general, the stress is defined as the derivative of the stored energy respect any convenient strain measure, for example, the logarithmic (Hencky) elastic strains H e := 1 2 log(F T e F e ) in the intermediate configuration, namely where D is the plastic dissipation and ctḢ e is the corrector rate ofḢ e (the one with the external power frozen). Plastic slip occurs in this glide system when τ g reaches a critical level τ c g (CRSS) that is typically considered as a function of the total accumulated plastic slip γ = |γ|dt and its rateγ = ∑ γ g , This CRSS value can be interpreted as the hardened yield stress of the slip system. In inviscid plasticityγ > 0 just when τ g > τ c g . In rate dependent plasticity, the slip rate in each glide system is then formulated as a function of the RSS τ g and the CRSS τ c g (e.g., as their ratio),γ Much work was performed to establish the specific form of Equation (18), e.g., [5,54,55,82,83], who have formulated a power-law function for slip rate as a common choice (other options are possible):γ beingγ 0 and m material parameters that quantify the reference shear rate and the rate sensitivity of the slip mechanism, respectively. τ c g is also referred to as strength or hardness of the slip system. In the above equationγ g can take both signs, while in the original work of Pierce et al. [82] the positive and negative slip directions are considered. The sensitivity exponent and the reference shear rate may include the influence of the temperature, for example, as [84,85] m = kT where k is the Boltzmann constant, T the absolute temperature, µ the shear modulus, τ f is the average strength of the forest, d is the distance to overcome an obstacle, ρ m is the density of the mobile dislocations (see below), ∆G is the stored Gibbs free energy when a dislocation bypasses an obstacle, and f D is the Debye characteristic frequency of the crystal, given by where N/V is the number of atoms per volume and v s is the speed of sound. Typical values are [85] For the hardening evolution of the CRSS τ c g in each glide system, a hardening law is typically defined. A formulation proposed by Peirce et al. [82] for isotropic type hardening is often used. It includes the explicit contributions of the occurring slip in the same glide system (self-hardening) and of the cross-effect from slips on all other glide systems (latent hardening):τ Here h gj is a matrix of hardening coefficients, fitted/determined empirically to capture the micro-mechanical interaction between the different slip systems [82]: The material parameter q differentiates latent hardening and self hardening of the glide systems; typically in an interval [1, 1.4]. For FCC single crystals, a typical choice for the hardening function h(γ) is [82]: being h 0 the initial hardening rate. τ 0 and τ s denote initial and saturation values for the CRSS. We note that Equations (19), (22) and (24) are not the only ones proposed in the literature. Some researchers [86][87][88][89][90][91] modify Equation (19) to take into account the cyclic deformation by introducing a phenomenological backstress term to formulate the kinematic hardening. Other authors [92,93] use classical hardening laws from continuum plasticity such as Voce's hardening law to avoid vanishing hardening rates though the introduction of a limiting slope.

Physics-Based Models
In contrast to the phenomenological models, the physics-based ones use physical (or physically motivated) measures, as well as principles and important physical insight from microstructure information and its evolution. Examples of that information are dislocation density, precipitate morphology, grain size and shape, etc. Nevertheless, only a few of these aspects have been introduced into physics-based models. In following subsections, we present more in detail some predominant models.
(a) Dislocation density-based models There are 10 6 to 10 12 dislocations in every cm 2 of material. The amount is variable and strongly depends on the processing history, in particular deformation and thermal treatments, possibly on microstructure [27,94]. Whereas dislocations are mostly responsible for plastic deformation, crystal plasticity models are usually developed at a different scale. Then, the different scales are considered throughout the formulations in terms of statistically stored dislocation (SSD) densities in the crystal plasticity models. At a given time, there are two types of SSD: mobile (glissile) and immobile (sessile), with respective densities ρ m g and ρ im g ; in total ρ g = ρ m g + ρ im g . These densities represent the equivalent length of dislocations (Burgers vectors) along a direction per unit volume, e.g., with units [m/m 3 ] = [m −2 ]. Mobile dislocations are the ones producing effective slips, but during deformation a change of status of dislocations is possible [95,96] which means a transfer between ρ m g and ρ im g . This decomposition has been proposed to develop a different treatment needed for both components, because mobile dislocations are the source of plastic deformation and ductility; whereas immobile dislocations block motion, strengthening or hardening the material. The laws of evolution for mobile dislocations require a coupling with the ones for immobile dislocations. The evolution laws for dislocation density is typically based in the following Kocks-Mecking model [95,97,98]: This law includes both storage and recovery effects (in a similar layout to the Armstrong-Frederick rule in hardening) that determine the variation of the dislocation density with plastic deformation. The first addend in the brackets accounts for the storage of dislocations 1/(b g l(ρ g )); the dynamic recovery by the term 2γ c ρ g is associated with the annihilation of stored dislocations. The function l(ρ g ) is interpreted as immobilization of mobile dislocations at obstacles after having travelled a characteristic (statistical average) distance. b g is the Burgers vector modulus of slip system g. The term γ c is the critical annihilation distance. Some modifications of Equation (25) are possible; refer, for example, to [99,100].
Further evolution laws for metals of different crystal structures is proposed in [101] and enhanced by subsequent works [102][103][104][105]. An example is the following which is a linear superposition of various mobile and immobile mechanisms: where ρ im sat is an experimentally determined saturation for immobile dislocation densities, The coefficients h s g , h m0 g , h im0 G , h m1 g , h im1 g and h r g depend on the degree of activity of the specific dislocation mechanisms g.
Regarding the evolution of the slip rate, in dislocation-based models the equation proposed by Orowan [106] is typically used instead of Equation (19) as it establishes the connection between the plastic slip rate on a given glide systemγ g and the mobile dislocation density, i.e., it creates a bridge to convert a continuum mechanical term into the physics of dislocations:γ g = ρ m g b g v g (28) where ρ m g is the mobile dislocation density, b g is the Burgers vector modulus of g and v g is the average velocity of mobile dislocations that depends on the details of the evolving dislocation mechanism. During slip, depending on whether thermal activation can contribute to get over the obstacle, the barriers to dislocation movement can be referred to as temperature-dependent or independent.
Therefore, the average dislocation velocity is also temperature-dependent [38,68,99,[107][108][109]. Generally, the average dislocation velocity v g can be expressed by the typical thermodynamic relation [110]: where λ g is the average distance between the obstacles in the glide system g (compare with Equations (19) and (20)). During slip, due to the crystalline structure, there is a local stress field associated to the distortion around each dislocation core. The Gibbs free energy in presence of a general array of obstacles in the slip plane can be formulated as [107] being ∆F the activation free energy necessary to overcome the obstacles without the aid of an applied shear stress, τ t the strength of the barrier at 0 K, p and q the parameters that are function of distance propagated by the dislocation, and < · > is the Macaulay brackets. Furthermore, it is well known that a crystal always contains a number of immobile dislocations, so if a mobile dislocation keeps moving, it interacts with those immobile dislocations and also with the local stress fields around those.
Immobile dislocations therefore represent obstacles to mobile dislocations, and their density is a key factor determining the threshold plastic stress (CRSS, τ c g ) required for dislocation motion within a glide system. In addition, the value of τ c g also depends on a variety of factors such as: point defects, crystal structure, temperature, and other internal variables [27]. Taylor [111] determined this threshold stress for a single slip to occur: where µ is the shear modulus. For several slip interactions and to include the influence of the temperature, Elkhodary et al. [112] proposed a modification of Equation (31): being τ 0 g the initial slip resistance on glide system g, nss the number of slip systems and a gi the Taylor coefficient that represents the strength of interaction between glide systems. T, T 0 are the current and reference temperature, respectively. ξ is the thermal softening exponent.

(b) Geometrically necessary dislocations (GND) models
One of the best known frameworks in crystal plasticity for FE analysis is the one from Kalidindi, Bronkhorst and Anand [39]. Here, the constitutive laws can be totally expressed from the loading history at a stress integration point and are defined as local models. However, experiments show that the mechanical behavior of single crystals is size-dependent and the classical continuum models of plasticity cannot reproduce the details of the behaviour of the material at such small scales, important in nanoindentation [113][114][115][116][117], bend-ing tests for small-scale beam [118][119][120][121][122][123] and micropillar compression [124]. Hall [31] and Petch [32] introduced the grain size dependence of the flow stress through an empirical relation. Basically, smaller grains result in larger yield stresses. Thereafter, numerous studies have demonstrated that this strengthening effect is due to the fact that in the surroundings of the grain boundaries there is a higher volume fraction of heterogeneous plastic deformation [62,[125][126][127], hence more apparent in smaller grains. Typically, the nonhomogeneuos plastic deformation takes place and it may affect strain gradients and orientation near the boundary. These gradients are associated to the geometrically necessary dislocations (GNDs) [125]. The GND concept appears in nonhomogeneous deformations and is incorporated in the constitutive framework by Nye's dislocation density tensor explained above, see Equation (5). For a given slip system g, given that the Schmid dyad remains constant in the intermediate configuration, we have for that systemḞ T = F T p L T p =γ g F T p · m g ⊗ s g anḋ where ρ GND g is the vector of geometrically necessary dislocations for slip system g, proportional to the Burgers vector. This vector can be decomposed in screw and edge dislocations. The screw dislocation is in the direction of the slip direction, whereas the edge dislocations may be decomposed in the direction perpendicular to the slip plane m g and in the direction in that plane which is perpendicular to the slip direction, i.e., in direction t g := m g × s g . These density rates arė which can also be written as [73,128]: The integration of GNDs into the constitutive model can be introduced by different ways in crystal plasticity. In some models, GNDs are added to the statistically stored dislocations (SSDs) and also contribute to hardening by the Taylor hardening relation [129,130]. In the case that the density of GNDs in each system is explicitly resolved, the critical resolved shear stress τ g can be formulated as [131]: where ρ SSD g is the SSDs density in system g; a SSD kg and a GND kg are coefficients that define the latent hardening interactions among the different glide systems. To obtain the evolution of GNDs, it is necessary to determine the Nye tensor or the slip rate at the integration point level. This task is typically done through an extrapolation of the internal variables from the integration points to nodes in each element [132,133], or alternatively by computing the gradients based on recovery techniques [134]. Another approach to include GNDs in high order models. The most popular model of this type was proposed in the works [65,72,135], and thereafter extended in many works [67,[136][137][138][139][140][141]. Here, the plastic slip in each system γ g is treated as a kinematic variable in the constitutive model. Another tensorial quantity defined as geometric-dislocation tensor G is proposed to measure GNDs and can be expressed in terms of pure dislocations.
Here ρ κ is a signed density and l κ ⊗ s κ is a dislocation dyad, where l κ is the line direction (either s k or m k × s k ). The relevant aspect in the model is that it focuses on the work for each independent kinematic process that is energy-balanced, and therefore, micro-forces conjugated with slip and its gradient are properly considered in the principle of virtual power. Moreover, for constitutive dependencies on the geometric dislocation tensor, G is decomposed in a typical elastic-strain energy and a specific so-called defect energy Ψ G . A simple defect energy is defined in a quadratic form of the type [65]: where λ is a scalar material parameter. The micro-force τ g gives the following viscoplastic yield conditions (Equation (7.19)) [65]: Phenom. slip hardening + additional hardening due to the energy stored by the GNDs (39) where H(γ g ) is a function characterizing rate-dependence, σ g a slip resistance, S g = s g ⊗ m g is the Schmid tensor, T |G = ∂Ψ G /∂G is defined as the thermodynamic defect stress. In Equation (39), there are two different hardening mechanisms. The first addend is purely phenomenological without backstress effects, whereas the second one is purely energetic from restrictions imposed by the thermodynamical framework, and furthermore, it represents a backstress on the slip system g.

(c) Continuum Dislocation Dynamic Models
The GND models presented above are in general based on the Kröner-Nye dislocation density tensor that can be used for modelling size effects. However, those models present a fundamental weakness when only GNDs are considered: after averaging, the dislocation density quantities only account for the GNDs, while all moving dislocations (redundant and macroscopically cancelling SSDs) must be taken into account. In this sense, GND models typically need some phenomenological assumptions in order to account for the contribution of the geometrically redundant dislocations for the deformation process. Such GND models only can provide a closed and kinematically consistent description of the dislocation dynamics if all dislocations share the same line direction [142]. In order to overcome this weakness, recent density-based theories of dislocation dynamics are developed by different researchers [142][143][144][145][146][147][148][149][150] which contribute to a new family of physicsbased plasticity model named continuum dislocation dynamics (CDD).
Here, as an example of such theories, we focus on the more recent CDD model introduced by Hochrainer [142,147,148,150]. This model is based on the description of connected dislocation lines by a generalized dislocation density tensor, i.e., modifies the definition of the dislocation density tensor. As mentioned above, models based on the Kröner-Nye tensor normally average over the line directions of all dislocations contained in a volume element. However, in the CDD models the dislocation density tensor can distinguish a priori dislocations by their line direction before any averaging is introduced, which permits the description of the kinematics of very general systems of dislocations. The modification of the dislocation density tensor can be achieved by defining it in a higher-order configuration space [142,147,148,150]. For the configuration space R × R × S (where S is the orientation space [0, 2π]), the modified dislocation density tensor takes the form where (r, φ) defines a point in R × R × S, ρ (r,φ) is the scalar dislocation density at this point (r, φ), b is the Burgers vector, and L (r,φ) defines the generalized line direction in the higher-order configuration space that depends on the canonical spatial line direction l (φ) = (cos φ, sin φ) and the curvature k (r,φ) . The evolution of this tensor can be obtained by the following relation [142,147,148] ∂ t α I I (r,φ) = −(V (r,φ) × α I I (r,φ) ) × ∇ (41) being V (r,φ) the vector of the generalized velocity in R × R × S, that is perpendicular to the generalized line direction. By ∂ t (•) we mean partial derivative of (•) respect to t. This generalized velocity consists of spatial component where v(r, φ) is the scalar velocity of dislocation segments at point (r, φ) and ν(φ) = (sin φ, cos φ), and an angular component θ ((r, φ)) that indicates the rotation velocity of these segments. The dislocation density ρ (r,φ) and the curvature k(r, φ) then can be formulated from Equation (41) by assuming only glide dislocation movements [147] ∂ t ρ = −(∇ · (ρv) + ∂ φ (ρθ)) + ρvk (42) This approach allows to address physically relevant microplasticity problems with simple deformation geometries [142,147,148], however the main problem found with this approach is the large computational cost, since there is a large number of degrees of freedom required for conducting simulations in a higher dimensional state space. In order to reduce the dimensionality, a Fourier expansion can be used for the dislocation density and curvature functions, and the expansion coefficients then are used as variables of a simplified CDD model [148]. The total dislocation density ρ t is defined as the zeroth-order expansion coefficient of the dislocation density function ρ (r,φ) , The total geometrically necessary dislocation density is defined as ρ G = (κ 1 ) 2 + (κ 2 ) 2 , where κ 1 and κ 2 are the first order terms of the Fourier expansion of ρ (r,φ) in the material point (r, φ), which are the components of the dislocation density vector κ (r) = (κ 1 , κ 2 , 0), From κ (r) , the classical Kröner-Nye tensor is recovered as α = κ (r) ⊗ b. Under the assumption that the scalar dislocation velocity and curvature fields are independent of the line orientation, the evolution laws for the total dislocation density ρ t and the dislocation density vector κ are given by where κ ⊥ = (κ 2 , −κ 1 , 0), m is the glide plane normal, and k is the average curvature. A phenomenological evolution equation is normally adapted for the average curvature [148].

Numerical Implementation Aspects
The finite element method (FEM) is the widely preferred method for computational simulation of solids, and that includes crystal plasticity. Finite element programs typically have a driver subroutine for material models, and commercial ones allow users to implement their own material models via user subroutines. In a typical analysis, the load is applied incrementally in "time" steps, and at each time increment the global equilibrium of the structure is reached by means of an iterative process using a global nonlinear solver (typically based on global Newton-Raphson methods). Once the problem is solved at time t, the crystal model must provide two important quantities at each material point for the solution for the next time increment: (i) the stress at time t + ∆t and (ii) the tangent moduli tensor C = dS/dE also at time t + ∆t (S is the second Piola-Kirchhoff stress ten-sor and E = 1 2 (F T F − I) is now the Green-Lagrange strain tensor, both in the reference, undeformed configuration).

Integration Algorithm
Usually, an elastic-predictor/plastic-corrector algorithm is used for the integration of the constitutive equation. In principle, all data at previous time t is known and the strain increments at time t + ∆t are given by the element subroutine calling the material one, so one can start "elastically"-predicting any of the involved quantities. Subsequently, the prediction is corrected by solving the nonlinear equations in residual form, typically using a Newton-Raphson scheme. The second Piola-Kirchhoff stress S [39], the slip ratė γ g [82], the elastic deformation gradient F e [151], the plastic deformation gradient F p [152] or the plastic velocity gradient L p [153] are suggested as starting point in the numerical implementation. An example of elastic predictor/return-mapping algorithm for single crystal models is summarized in Algorithm 1. In the case of viscous-type rate constitutive equations, since no yield function is employed, no check for plastic slip is needed (there is always a plastic slip, even if it is very small). However, convergence problems are common if the rate sensitivity parameter m (see Equation (19)) is too small (as to approach the rate-independent solution), so it is common to increase this parameter artificially just to overcome these numerical problems. On the other hand, in the rate-independent case, non-uniqueness of solution issues may arise (the final stress is unique, but not the active set of slip systems).

Algorithm 1 Typical structure of an elastic-predictor/plastic-corrector return-mapping algorithm
Require: Given the state variables at time t and the total strain at time t + ∆t. Solve iteratively the residues of the variables used in prediction t+∆t R = 0 using the Newton-Raphson method iii. Update the state variables. Compute the elastoplastic tangent matrix during this phase.

4: EXIT
By using the Newton-Raphson method to obtain the solution, it is necessary to compute the local Jacobian matrix and its inverse. The dimension of this matrix is equal to the dimension of the residues vector, i.e the number of independent variables (and equations) that is used as predictor, for example, six variables for S, eight for F p , etc. This dimension can be a relevant number if the slip rates are chosen (for example, 48 × 48 for BCC crystals). Taking into account that this system is solved several times until local convergence for each global iteration and at each integration point, it requires some efforts in reducing the number of active slip systems [8], because this task may be the most computational demanding in solving the structure.

Type of Elements
In most crystal plasticity simulations, tri-linear elements are used, because they are easily generated and are inexpensive from a material standpoint (typically 2 × 2 × 2 integration points are employed. However, these elements can not describe strain gradients within one element due to the fact that this type of elements use linear interpolation functions for the displacements. Furthermore, if they employ standard formulation, they are well-known to lock for the incompressible case typical in plasticity. Linear elements are especially not sufficient to capture accurately strong strain gradients. Higher-order elements should be used in such cases with either reduced integration or a formulation to alleviate locking. Typical elements include mixed u/p formulations, incompatible models formulations, or reduced integration. Details on these formulations can be found in [2]. In cases with GNDs-based models, the situation becomes more complicated when those models include strain gradients. The formulations for standard elements are only continuous in the displacements (C 0 -continuity) and gradients are often undefined. Enhanced element formulations have been proposed to be used to overcome this problem [67,154]. However, in the case of complex loadings the complexity in defining boundary conditions increases for such element formulations. Furthermore, these element formulations are computationally very costly, and lack of robustness at large strains, so they are often applied to two dimensional problems at small strains. Therefore, the standard elements are still used and the necessary strain gradients are derived based on recovery techniques [134].

Searching the Active Set of Slip Systems
The question whether a specific slip system is activated is decided by a yield criterion mathematically expressed in terms of one inequality per slip system. It essentially states that plastic slip occurs if the RSS reaches the CRSS subject to the additional constraint that only slip rates in the positive and negative direction of the RSS are permitted. The overall constitutive material response can be determined once the set of active slip systems is known. However, it is not possible, in general, to know which systems will be active for the final converged deformation. Noteworthy, for multisurface plasticity models such as the Tresca and Mohr-Coulomb models there are robust algorithms to search the active set of yield surfaces. Those algorithms have been formulated based on the geometric representation in the principal stress space. Nonetheless, in the present crystal plasticity case, its application is more complex because the set of active slip systems that satisfies the discrete plastic consistency may not be unique and many combinations of plastic multipliers may exits for a given set of active slip systems, which provide the same incremental plastic deformation gradient [5,155]. If the set of Schmid tensors S g = s g ⊗ m g is linearly dependent in the space of deviatoric tensors, non-uniqueness in the solution is possible. In 3D there are 8 linearly independent tensors that generate the deviatoric tensors space. If it happens that the total number of active slip systems is larger than eight, then the Jacobian matrix of the return-mapping system of equations becomes singular. Therefore, an additional condition has to be formulated to be able to select the set of active slip systems which appears in reality. Typically a Moore-Penrose (minimum norm) solution is chosen, although other options are possible [156].

Rate-Dependent Approach
One of the preferred solutions to overcome the problems related to non-uniqueness are viscoplastic (rate-dependent) formulations. These have been introduced firstly in the works [82,83,157], in which the slip rate is given employing a typical power law function as shown in Equation (19). This approach is then popularly used not only in continuum crystal plasticity models [86][87][88][89][90][91]93] but also in the homogenization of polycrystals [158][159][160][161][162][163][164] or in multiscale modelling of polycrystals [165][166][167][168]. As highlighted by [82], this rate-dependent model allows a more accurate large strain formulation for polycrystals. In particular, it can obtain a good prediction of texture and specially obtain its dependence on both strain-hardening/latent-hardening and strain rate sensitivity. But much of the choice of the rate-dependent formulation is to avoid the non-uniqueness problem of rate-independent models, a formulation which can be recovered by using reasonably low values for the ratesensitivity constant (i.e., a high exponent 1/m in the slip law of Equation (19)). However, in this case, the resulting set of equations that are integrated can become extremely stiff as the rate-sensitivity constant reaches a very small value, so a robust formulation for those cases is not simple. These issues have been addressed and tackled in the work of Cuitiño and Ortiz [8].

Rate-Independent Approach
Borja and Wren [155] proposed a completely different approach. Restricted to the infinitesimal theory, a fairly robust algorithm (the ultimate algorithm) for the selection of the active slip systems was developed within the stress-update procedure for rate-independent single crystal models. The basic idea for determining the set of active slip systems is to start the application of the predictor/return-mapping scheme with some trial set of active slip systems, and then, to generate a sequence of sets of active slip systems as the Newton-Raphson iterations for obtaining the solution are applied. At the end of this rootfinding process, the algorithm will converge to a state when the discrete complementary condition is satisfied. The set of active slip systems is the converged set and is unique. Furthermore, to avoid the possibility that the Jacobian matrix can be singular, Borja and Wren proposed the triangular factorization of the Jacobian matrix after the redundant equations are eliminated. This algorithm has been extended recently for the finite strain framework by Borja and Rahmani [169]. Miehe also used this algorithm in the finite strain context [156,170], but introducing the concept of exponential map integrators of [171] into the field of single crystal simulation.

Engineering Forming Problems
The above formulations have been employed in a wide number of problems where, for example, texture and its evolution may be important. These have been compared to experimental data and to continuum models to show comparatively enhanced predictions. A typical example is a cylindrical cup deep drawing adopted to investigate the texture effects on the localization of the deformations and in earing [172]. In order to do it, the texture components were incorporated into the crystal plasticity model by introducing the measured orientation distribution function (ODF) at each integration point [173]. Three typical crystal orientations {111} 110 , {113} 110 and {001} 110 were selected from the dominant textures in three steel sheets (mild steel Deep-Drawing Quality DDQ dual-phase steel-DP600Ze and high-strength steel CP800). The results obtained from simulations are compared to experimental results. Figure 6 shows the computed results compared with one obtained in experiments for the high-strength steel (CP800) material. The predicted and measured thickness strain distributions (very important in deep drawing), along the rolling direction are similar (see Figure 6c). The thinning of the material appears at the hemispherical bottom region, where a high strain localization can also be observed at the bottom region. Figure 6a,b show the localized necking precluding fracture in numerical and experimental results, respectively. The expected four ears at 45 • to the rolling direction are observed. With the fitted simplified texture, the texture component-based crystal model give already results very close to the experimental ones, demonstrating the value of this type of formulations. Furthermore, the results obtained in the study [172] also confirmed that the more the γ fiber texture components, the better is the drawability of bcc steel sheet. The γ fiber also gives rise to four ears at 0 and 90 • to rolling direction. Figure 6. Comparison of the cup drawing between numerical simulation and experiment for CP800 material: (a) predicted thickness strain distribution and earing at punch travel 10 mm; (b) photographs of experiment specimen; (c) thickness strain distribution along the rolling direction. Reproduced from [172], under permission.
The same simulation benchmark was also employed in [174] but using a new concept for representing textures by mapping small sets of mathematically compact spherical Gaussian texture components on the integration points [175][176][177]. Figure 7 shows the comparison between the computed and experimental results for cube-textured aluminium sheet after cup drawing. In Figure 7a, the result obtained using the texture componentbased crystal model is compared with the experiment and also with continuum-based elastoplasticity obtained by use of a Hill-48 yield surface. Two texture components and a random scattering background component are used for fitting the texture component. These results show a better fit of the texture component-based crystal model than that obtained with the classical Hill-based model. Figure 7b,c show simulation results for different approximations using a volume fraction and using two rolling texture components with a random texture component.
Using the elasto-plastic self-consistent homogenization scheme in the finite element framework, Zecevic and Knezevic [178] explored its potential in simulation and prediction of springback of cup drawing from an AA6022-T4 sheet, see Figure 8; only a quarter is modelled due to symmetries. The springback is obtained by eliminating the contact between the blank and other parts when the punch reaches a travel distance of 12 mm. The von Mises stress contours at the end of drawing and springback are shown in Figure 8a,b, while the dimensional changes are depicted in Figure 8c,d.

Texture Evolution
Obviously one of the simulations that can be performed with crystal plasticity models is the analysis of texture evolution and, hence, the evolution of anisotropy. In Figure 9 we show the texture evolution in a copper polycrystal cubic specimen in a tensile test, see full details about the constitutive equations (similar to the avove-presented phenomenological ones) and the material parameters employed in the simulation in [75]. For this simulation we used 512 elements with an initial random texture with 50 orientations. In the same figure, the initial mesh and the deformed mesh is shown. In the middle and lower rows of the figure, the pole figures showing the evolution of texture for directions (111) (middle row) and (100) (lower row) are shown. This evolution has been computed in Figure 9 using both the classical Kalidindi et al. framework and the novel framework employing elastic corrector rates, being the results very similar for this case (see again details in [75]).
After Before  Figure 9. Evolution of texture in copper under a tensile test. Upper row contains the initial and final mesh, where colors correspond to different orientation. Middle row is shows the texture pole plots evolution for direction (111) and the lower row for direction (100). Adapted from [75], under permission. A Matlab toolbox for pole figures, and documentation, is available in https://mtex-toolbox.github.io/.

Impact Test
Taylor impact tests were originally devised to obtain the dynamic yield strength of a material at moderate strain rates, but now such tests are used frequently to validate the constitutive model for the simulation of plastic deformation. Figure 10a shows simulations of Ta cylinder by using the viscoplastic self-consistent (VPSC) model in the finite element framework. The simulation was performed by [179] and compared with the experimental results reported in [180,181]. Figure 10a shows the equivalent plastic strain rate contours that develop in the cylinder during the simulation. It can be seen that the dramatic change in the strain rate is observed over the course of simulation. At the beginning of the impact, the highest equivalent plastic strain rate develops at the center of the foot of the cylinder. As the impact progresses, the high strain rate zone moves away from the foot to the other side of the cylinder as a consequence of the plastic wave propagation. The side profile and footprint obtained in the simulation after the impact are compared with the experimental ones and depicted in Figure 10b. The highly heterogeneous plastic deformation is observed at the center of the foot. The deformation in lateral directions is non-uniform as shown in the ellipsoidal shape of the footprint.

Rolling Test
Using the multiscale framework, a typical rolling test is normally carried out to explore the mechanical behavior of polycrystals subjected to non-homogeneous stress states and examine heterogeneous texture development through the thickness of the sample during the process [165,182]. The rolling simulation that has been done in [182] is here recalled. The material of rolling plate is uranium α-U (see [182] for detailed mechanical and texture properties). The initial texture was represented using 1600 individual orientations. Each orientation was assigned to two elements and distributed over their eight integration points. The sample was rolled between two cylinders leading to 60% reduction in thickness during five steps. The geometry of the rolled plate at each step is depicted in Figure 11a, together with the corresponding contour of the equivalent plastic strain. Due to the highly anisotropic nature of the uranium α-U, the strain is developed heterogeneously through the thickness. Similarly the heterogeneous texture that can be appreciated in Figure 11b.
According to [182], the determination of the texture gradients can help with understanding recrystallization kinetics and the final recrystallized microstruture. Consequently, the texture gradients can affect recrystallization by altering crystallographic orientation, grain boundary and dislocation densities. Figure 11. Simulation results of rolling test for the uranium α-U: (a) equivalent plastic strain contours at different step; (b) comparison between measured and predicted evolution at a polycrystalline material point near the bottom surface of the plate. Reproduced from [182], under permission.

Conclusions
In the 21th century, computational crystal plasticity modelling has achieved a high degree of maturity. Thanks to the current computational power and to the algorithms with increased efficiency as multiscale ones based on the FFT, it is becoming an alternative to classical continuum plasticity for engineering problems. Crystal plasticity simulations allow for incorporating many ingredients of the structure of the material as grain size, dislocations, texture, etc. In this paper we have made an overview of different aspects of the present techniques employed in computational finite element simulations of crystal plasticity, as well as some applications of such simulations.