On the computational solution of vector-density based continuum dislocation dynamics models: a comparison of two plastic distortion and stress update algorithms
Peng Lin, Vignesh Vivekanandan, Kyle Starkey, Benjamin Anglin, Clint Geller, Anter El-Azab
OOn the computational solution of vector-density based continuum dislocation dynamics models: a comparison of two plastic distortion and stress update algorithms
Peng Lin , Vignesh Vivekanandan , Kyle Starkey , Benjamin Anglin , Clint Geller , Anter El-Azab Purdue University, West Lafayette, IN 47907, USA Naval Nuclear Laboratory, West Mifflin, PA 15122, USA
Abstract
Continuum dislocation dynamics models of mesoscale plasticity consist of dislocation transport-reaction equations coupled with crystal mechanics equations. The coupling between these two sets of equations is such that dislocation transport gives rise to the evolution of plastic distortion (strain), while the evolution of the latter fixes the stress from which the dislocation velocity field is found via a mobility law. Earlier solutions of these equations employed a staggered solution scheme for the two sets of equations in which the plastic distortion was updated via time integration of its rate, as found from Orowan’s law. In this work, we show that such a direct time integration scheme can suffer from accumulation of numerical errors. We introduce an alternative scheme based on field dislocation mechanics that ensures consistency between the plastic distortion and the dislocation content in the crystal. The new scheme is based on calculating the compatible and incompatible parts of the plastic distortion separately, and the incompatible part is calculated from the current dislocation density field. Stress field and dislocation transport calculations were implemented within a finite element based discretization of the governing equations, with the crystal mechanics part solved by a conventional Galerkin method and the dislocation transport equations by the least squares method. A simple test is first performed to show the accuracy of the two schemes for updating the plastic distortion, which shows that the solution method based on field dislocation mechanics is more accurate. This method then was used to simulate an austenitic steel crystal under uniaxial loading and multiple slip conditions. By considering dislocation interactions caused by junctions, a hardening rate similar to discrete dislocation dynamics simulation results was obtained. The simulations show that dislocations exhibit some self-organized structures as the strain is increased. age | 2
Keywords:
Continuum Dislocation Dynamics; Field Dislocation Mechanics; Plastic Distortion; Cross Slip; Junction Reaction; Finite Element Method age | 3 Introduction
The continuum dislocation dynamics approach to crystal deformation has gained a great deal of attention lately (El-Azab, 2000; Groma et al., 2003; Hochrainer, 2016; Hochrainer et al., 2014; Monavari et al., 2016; Monavari and Zaiser, 2018; Sudmanns et al., 2019; Xia and El-Azab, 2015; Yefimov et al., 2004). This approach is attractive due to the fact that one can use it to integrate the mechanics and physics of dislocations into continuum mechanics descriptions of metal deformation in much the same manner as with discrete dislocation dynamics methods (Cui et al., 2019, 2014; Devincre et al., 2008, 2006; El-Awady, 2015; Liu et al., 2017; Sills et al., 2018; Stricker et al., 2018; Stricker and Weygand, 2015). While the latter methods proved to be powerful in tackling important mesoscale plasticity problems (Cui et al., 2014; Devincre et al., 2008; Liu et al., 2009; Mishra and Alankar, 2019; Stricker and Weygand, 2015), they remain computationally prohibitive from the viewpoint of bulk deformation studies due to the expense of calculating long-range interactions. Such computational difficulties are perceived to be surmountable with continuum dislocation dynamics, since dislocation-dislocation interactions are computed indirectly via an eigenstrain approach (Xia and El-Azab, 2015; Yefimov et al., 2004). A typical continuum dislocation dynamics formalism consists of two coupled sets of governing equations (Xia and El-Azab, 2015; Yefimov et al., 2004). The first includes the conventional crystal mechanics equations cast in the form of stress equilibrium and deformation kinematics, with the plastic strain to be fixed from the dislocation motion, and the second comprises the equations governing the transport of and the reactions between dislocations. These two sub-problems are two-way coupled in that the plastic strain is fixed from the dislocation fluxes via Orowan’s law, while dislocation motion is dictated by the local internal stress via a dislocation mobility law. A very important aspect of the continuum dislocation dynamics model is the representation of dislocations by multiple density fields distinguished based upon their Burgers vectors and slip planes. Such a representation enables the accurate calculation of both dislocation motion and reactions (El-Azab, 2000; Groma et al., 2003; Hochrainer et al., 2014; Monavari et al., 2016; Monavari and Zaiser, 2018; Sudmanns et al., 2019; Xia and El-Azab, 2015). The aforementioned representation of dislocations in continuum dislocation dynamics is an extension of the earlier work of Kröner (Kröner, 1958), and Nye (Nye, 1953), who established a theoretical framework that can inform a generic relationship between the dislocation microstructure, the plastic distortion, and the associated internal stress fields in plastically age | 4 deformed crystals. In such a framework, dislocation density is measured by a dislocation density tensor, , which is defined by the curl of the plastic distortion, p = − , with the plastic distortion, p , and its elastic counterpart, e , serving as the fundamental components of the deformation of crystals, in the sense that their sum forms the gradient of a compatible displacement field, u . The time evolution of the dislocation density tensor is found from ( ) v = , which was formulated by Mura (Mura, 1963) and Kosevich (Kosevich, 1965). This formula, however, can only be applied to families of dislocations of the same Burgers vector on a relatively small scale since the dislocation velocity v is only meaningful at that level. In recent years, several attempts have been made to refine an average, statistical description of the evolution of dislocation microstructures. Pioneering work was done for straight, parallel dislocations in two dimensions (2D) by Groma (Groma, 1997), Zaiser et al (Zaiser et al., 2001), and Groma et al (Groma et al., 2003). The joint evolution of statistically stored and geometric dislocation densities can be captured by these models; see also the recent work of Kooiman et al (Kooiman et al., 2014) and Finel et al (Cottura et al., 2016). It is quite challenging to extend such 2D approaches to 3D systems where dislocations are interconnected, curved lines moving perpendicular to their line direction. Different approaches have been established to represent 3D dislocation microstructures, for example, by solving for the evolution of the scalar density ( ) , t x and orientation function ( ) , t x (Sedláek et al., 2007) or the screw dislocation density, ( ) screw , t x , and the edge dislocation density, ( ) edge , t x (Arsenlis et al., 2004; Grilli et al., 2018; Leung et al., 2015). Another approach proposed by Hochrainer et al (Hochrainer et al., 2007) is based on picturing the 3D curved dislocation lines in a higher dimensional space containing line orientation as an additional dimension, so that the dislocation density can carry additional information about the line direction and curvature. Simplified variants of this theory have been formulated, which consider only low-order moments of the dislocation orientation distribution (Sandfeld et al., 2011; Sandfeld and Zaiser, 2015). A further development of this theory is achieved through a hierarchy of evolution equations of the so-called alignment tensors, and the tensors contain information on the directional distribution of dislocation density and dislocation curvature (Hochrainer, 2015; Monavari and Zaiser, 2018). Generally speaking, the implementation of CDD is challenging since the numerical solution of the age | 5 transport equations always suffers some level of dispersion and diffusion; see for example (Sandfeld et al., 2015; Schulz et al., 2019). In this paper, a different description of dislocation evolution that was recently formulated by Xia and El-Azab (Xia and El-Azab, 2015) is adopted. In this description, dislocations on various slip systems are represented by vector fields ( ) , with the superscript being the slip system index. The dislocation tensor is decomposed according to the contributions of the dislocation vectors ( ) from various slip systems weighted by their corresponding crystallographic Burgers vectors, ( ) b . Here the vector fields ( ) locally represent dislocation bundles with the same line direction. In this case, the direction of the dislocation velocity vector ( ) v for each slip system is uniquely determined as the direction perpendicular to the dislocation line direction. To complete a full description of crystal plasticity, the stress field applied to the dislocations is needed to fix the dislocation velocity. Different approaches have been adopted in the literature to calculate the long-range interactions arising from dislocation networks. One approach is to develop an analytical expression for the elastic stress between two infinitesimally short dislocation segments based on fundamental dislocation theories (Balluffi, 2016; Hirth and Lothe, 1982). Another approach is to calculate the plastic distortion p (or the plastic strain p ) and use it as an eigenstrain with which to solve the stress equilibrium equation (Evers et al., 2004; Lin et al., 2016, 2015; Xia and El-Azab, 2015; Yefimov et al., 2004). The latter approach is more commonly used in crystal plasticity, because the stress imposed by boundary conditions can also be obtained equivalently by solving the stress equilibrium equation. The evolution of plastic distortion itself is evaluated using Orowan’s law, p ( ) ( ) = − v . When the dislocation density field is adequately smooth in space, updating plastic distortion by direct use of Orowan’s law provides accurate results. This is the case for most crystal plasticity problems of interest at the mesoscale. However, in this paper, we show that, when dislocation density fluctuates rapidly in space, as in models that aim to capture dislocation patterns, numerical errors in p can accumulate, resulting in a significant discrepancy between p − and . We further show that such errors tend to propagate into the predicted stress-strain behavior and dislocation patterns. age | 6 In the current work, we use the field dislocation mechanics approach (Acharya, 2004; Acharya and Roy, 2006) in conjunction with our vector density-based dislocation dynamics formulation. The main idea here is to express p as the sum of two components: p = − z , with z and being compatible and incompatible components of p , respectively. In this way, the incompatible part is solved for directly from the current dislocation configuration represented by . It is thus guaranteed that the incompatible part of the plastic distortion p tensor is consistent with the dislocation density tensor , as further explained in the next section. In Section 2, we introduce the basic formulation of continuum dislocation dynamics in terms of the dislocation vector density and its evolution equations. In Section 3, the algorithms for updating plastic distortion are presented, including the direct time integration algorithm and the field dislocation mechanics formulation. In Section 4, we present the numerical implementation scheme. In Section 5, we present the results produced by the two plastic distortion/stress update algorithms for the case of an austenitic steel crystal under uniaxial loading and multiple slip conditions. We provide brief concluding remarks in the final Section 6. Continuum dislocation dynamics model
The concept of distortion is fundamental to the mechanics of deformed crystals (Kröner, 1958). In the case of small deformation, the distortion of a crystal, , corresponds to the displacement gradient, u . This distortion has two components, elastic, e , and plastic, p . We may thus write: e p = = + u . (1) Due to the presence of dislocations, e and p are incompatible fields and thus do not correspond to displacement fields individually. The continuum dislocation dynamics approach is based on the premise that the dislocation system in a deforming crystal can be modeled as a continuous field. To solidify the concepts being discussed, we first consider the case of a single active slip system with all dislocations having the same Burgers vector. With dislocations characterized by two vectors, the line direction and Burgers vector, a dislocation field is often described in the continuum theory by a second order tensor field, , defined in terms of the plastic distortion by (Kröner, 1958): p = − . (2) age | 7 On the other hand, the evolution of the plastic distortion is caused by the motion of dislocations at rate p obtained by Orowan’s law, p − v = , (3) where v is the velocity field of dislocations. We remark here that the above equation is valid for a single slip situation and at sufficiently small spatial resolution permitting the use of Orowan’s law (3). Taking the time derivative of equation (2) and combining it with equation (3), the equation for the time rate of change of the dislocation density tensor is formed (Mura, 1963), ( ) v = . (4) The above-stated requirement of sufficiently fine spatial resolution means that, on each slip system the dislocations have the same line direction at a given (continuum) point in space. The dislocation density field then can be represented by a vector field, , which carries both the number density and line orientation at every point in space. The dislocation density tensor field is related to the vector density according to equation (5) (Gurtin et al., 2007) b = , (5) where b is the Burgers vector of dislocations. In order to obtain the evolution equation of the vector density, , equation (5) is substituted into equation (4), and the Burgers vector is dropped, yielding ( ) v = , (6) as previously shown by Xia and El-Azab (Xia and El-Azab, 2015). Equation (6) is considered to be the main kinetic equation of dislocation dynamics in our model. It is itself the evolution equation of dislocations on each slip system. It describes the motion and bowing of a dislocation line, and hence the line length increases in a natural way. Also, it can easily be modified to incorporate cross slip by allowing the direction of the velocity vector to change from one slip plane to another. As pointed out in the introduction, that transport equation (6) is based on the line bundle idealization of dislocations (Xia and El-Azab, 2015). This idealization assumes that at a continuum point on a given slip systems, the density has only a single line direction. To satisfy this assumption, the mesh size should be comparable to the annihilation distance of dislocations of opposite sign, and as such the geometric cancellation of dislocations of opposite sign coincides with their physical annihilation. age | 8 The vector dislocation density is a combination of the scalar density and the unit tangent to the dislocation line , = . (7) The dislocation velocity is also represented by a vector field, written in the form v = v . (8) As the dislocation line moves perpendicular to its line direction, the dislocation velocity direction can be determined by the slip plane normal m and the line direction as follows: m = . (9) The scalar velocity v is determined by the resolved shear stress applied on dislocations. According to (Zaiser et al., 2001) and (Yefimov et al., 2004), the dislocation glide velocity is assumed to change linearly with the local resolved shear stress, ( ) ( ) sgn p bv B = − + (10) where b is the magnitude of Burgers vector and B is a drag coefficient. The function “sgn” returns the signature of its argument, denotes a Macaulay bracket, which returns its argument if it is positive and zero otherwise, is the stress representing lattice friction, and is the resolved shear stress corresponding to the Peach-Koehler glide force with magnitude b (Peach and Koehler, 1950). This expression for v includes the long-range interactions between dislocations, together with the stress arising from boundary conditions. Since dislocation reactions, e.g. junctions, can lock or impede the motion of dislocations, a short-range interaction stress is included via p . The details are discussed later. The dislocation transport equation (6) can be solved when the dislocation velocity field v is prescribed. As has just been mentioned, the velocity field comes from the resolved shear stress or the Peach-Koehler force according to equation (10), both of which require fixing the stress field, , which is obtained by solving the stress equilibrium equation. The latter equation is part of the crystal mechanics problem, cast below as an eigenstrain problem: age | 9 p sym in : ( ) in on on u = = − =
0c uu un t = , (11) in which c is the elastic tensor. In the above, the plastic distortion (or its symmetric part, the plastic strain) is considered as the eigenstrain. The numerical update of eigenstrain (eigendistortion) represents the main motivation of this work and it is explained in Section 3. In equation (11), is the crystal domain of the problem and its boundary is , with the latter partitioned into parts u over which the displacement boundary condition u is applied, and over which a traction boundary condition t is applied, where n is the unit outward normal to the boundary. By solving equation (11), the displacement and stress fields can be obtained from both imposed systems of boundary tractions and internal dislocations (eigenstrain). The Peach-Koehler force, PK f , acting on the dislocations in a stress field, , is given by (Hirth and Lothe, 1982) ( ) PK = f b , (12) where PK f is the Peach-Koehler force per unit length of dislocation, b is the Burgers vector, and is the local tangent of the dislocation bundle. Equation (12) can be rewritten to obtain the glide component and the climb component. Based on the definition of the velocity direction (equation (9)), we have m = . So ( ) ( ) ( ) ( ) PK = = f b m b m b m − . (13) Equation (13) shows the glide component and climb component of the Peach-Koehler force clearly, ( ) ( ) g c and b b = = f s m f s m − . (14) In the above, s is the unit slip direction / | | = s b b . Only the glide component of the motion is considered in the dislocation transport, which corresponds to a resolved shear stress of the form = s m , (15) to be used in equation (10). A generalization to a multi-slip situation requires equations (6), (10), (13) and (15) for the dislocation transport, dislocation velocity, Peach-Koehler force and resolved shear stress, age | 10 respectively, to be specialized to individual slip systems. This in turn requires the dislocation density, , velocity, v , slip plane normal, m , and slip direction, s , Burgers vector, b , resolved shear stress, , Peach-Koehler force, PK f , and its glide and climb components, and the short range stress, p , to be so specified. We will do so by adding superscripts, , , etc. As dislocations evolve, dislocation reactions and cross slip take place. In the current formulation, some reactions are explicitly considered such as the glissile junction formation and collinear annihilation. Others such as the immobile Lomer and Hirth junctions, which contribute part of the hardening behavior, are modeled via a short-range, Taylor type term, p , in the mobility law, Eq. (10). The form of this term is taken after (Devincre et al., 2006; Franciosi et al., 1980; Kubin et al., 2008): p b a = , (16) with being the shear modulus and a an interaction coefficient representing the average strength of the interaction between slip systems and . Dislocation cross slip, collinear annihilation, and glissile junction reactions, influence the balance of dislocations on various slip systems and thus lead to additional terms in the dislocation transport equation (Eq. (6)). The latter thus is modified by adding an extra term, cp , representing the contributions of these processes. That is, ( ) cp + v = . (17) Cross slip and reactions change the evolution of a dislocation network and lead to dislocation multiplication. The cross slip and collinear annihilation terms with the same Burgers vector will be coupled. Glissile junction reactions couple three slip systems, say, , , , where Burgers vectors satisfy + = b b b . A detailed formulation of these coupling terms can be found in an earlier work by the current authors (Lin and El-Azab, 2020). Plastic distortion and stress update algorithms
We now turn to the update algorithm for the plastic distortion tensor, p , which is required to calculate the long-range stress on the dislocation system subject to the prescribed boundary age | 11 conditions. As discussed earlier, p serves as an eigenstrain in the crystal mechanics problem described in equation (11). The dislocation dynamics problem yields the rate of the plastic distortion in terms of Orowan’s law, equation (3), which upon combining with equations (7) through (9) and considering a multi-slip situation, gives ( ) p v v b − = = b m s m s = , (18) with v b = being the slip rate or the time derivative of the local shear strain on slip system . We remark here that, in crystal plasticity models of small strain, the plastic distortion p is expressed in terms of crystal slips as follows: p = m s . (19) Algorithm 1: Update of p by direct time integration of p : From equations (18) and (19), it can be seen that, to update the plastic distortion tensor p , it is only necessary to calculate the crystal slip magnitudes for all slip systems by integrating the corresponding slip rates. That is, ( ) ( , ') ' t t dt = x x . (20) The plastic distortion can then be constructed by inserting all values from equation (20) into equation (19). As mentioned in the introduction, this algorithm results in an inconsistency between the dislocation density tensor, which is the measure of the elastic strain incompatibility, and the curl of the updated plastic distortion. This inconsistency arises due to accumulation of the numerical errors in the staggered solution schemes of the coupled dislocation transport/crystal mechanics equations. Such errors result in a stress field that is inconsistent with the dislocation tensor itself, the latter being the origin of stress. The same inconsistency also appears in differences between the slip system dislocation density and the curls of the ‘oriented’ plastic slip (shear strain) on individual slip systems. The latter is the plastic shear strain taken together with the unit normal to the slip plane for a given slip system. Algorithm 2: Update of p using the field dislocation mechanics approach : In order to obtain the stress field accurately, we adopt the field dislocation mechanics formulation (Acharya and Roy, age | 12
In field dislocation mechanics, the plastic distortion p is divided into a compatible part z and an incompatible part , p = − z . (21) Being the gradient of a vector field z , the compatible part is curl-free. The incompatible part, on the other hand is, divergence-free. The field problem governing the incompatible part, , can be derived from equations (2) and (21), in in on n = = = (22) The vector field z associated with the compatible part of p satisfies a boundary value problem derived from Eqs. (3) and (21), ( )( )( ) in on arbitrary value at one point in o = − = = z vn z n vz z − (23) The boundary conditions in equations (22) and (23) are generally applicable and have been discussed in detail by Acharya (Acharya and Roy, 2006). They are used to obtain a unique decomposition of the plastic distortion tensor. Because we concern ourselves with the total plastic distortion, we first assemble the dislocation density tensor, b = . We then use it in equation (22) and we substitute the Orowan’s term, v , into Eq (23). Upon solving the last two boundary value problems for the compatible and incompatible parts of the plastic distortion tensor, p , the latter can be assembled as per equation (21). We remark here that, although the compatible part still needs integration over time, the incompatible part will be calculated directly from the current dislocation density tensor regardless of the history of dislocation evolution. We will show the importance of calculating the incompatible part independently of the evolution history in a test example later. It should be pointed out that only the decomposition of plastic distortion of Acharya and Roy’s field dislocation mechanics formalism (Acharya and Roy, 2006) into compatible and incompatible parts is adopted here. Our dislocation transport equations (Eq. (17)) are different. While Acharya and Roy use a second-order dislocation density tensor to represent the dislocation field, a set of age | 13 slip-system level dislocation vector densities are used in our case. This enables us to compute the dislocation velocity from a Peach-Koehler force via a mobility law, and it enables us to explicitly incorporate cross slip and reactions along with their line direction dependence. Numerical Solution Scheme
As presented above, the overall problem of continuum dislocation dynamics can be divided into two parts, solving the stress equilibrium equation (11) for the displacement field u given the plastic distortion distribution, together with the dislocation kinetic equations (6) or (17) for the set of . With Algorithm 2 for the plastic distortion update, two additional sets of boundary value problems, summarized in equations (22) and (23), are to be solved for z and . The latter can be considered as subsidiary problems to be solved for the proper update of the plastic distortion and, in turn, the stress tensor. We employ a staggered scheme to solve the coupled crystal mechanics/dislocation kinetics problem. We first assign initial densities to all slip systems within the simulation domain with an initially known plastic distortion field p . We next solve equation (22) and (23) to obtain the corresponding plastic distortion p . Then the stress field can be obtained by using the stress equilibrium equation (11). Once the stress field is found, the resolved shear stress equation (15) and the dislocation velocity equation (10) can be used to solve the kinetic equation for dislocation density evolution, per equation (17). The new dislocation configuration is used for the next time step until the average strain reaches a desired value. The overall scheme is schematically shown in Figure 1. age | 14 Figure 1. The flowchart for the solution scheme. is the applied strain during loading and * is the final strain. We use a standard Galerkin finite element method to solve the stress equilibrium equation and the boundary value problem for the compatible part of plastic distortion, while, for ease of enforcing the divergence-free constraints, the least squares finite element method is used for the kinetic equation for dislocation density and the incompatible part of plastic distortion. Details of the finite element formulations of these equations are in Appendix. Results and discussions
To illustrate the improvement in the manner of updating plastic distortion by using field dislocation mechanics instead of direct time integration of a rate expression, we show a test example in which two dislocation loops expand with a prescribed velocity and then annihilate with each other. Then a bulk simulation for an FCC crystal (austenitic steel) with [001] loading, performed by the age | 15 continuum dislocation dynamics with field dislocation mechanics, is shown to illustrate the microstructure evolution during loading.
Expansion and annihilation of two dislocation loops
The two alternative algorithms for updating plastic distortion are compared in this example. Initially, two dislocation bundles in the form of loops were placed close together in a simulation domain of size 5µm×10µm×5.303µm, as shown in Figure 2. Dislocation line bundles are represented by a vector field initialized by a Gaussian distribution function over their cross section. The values in Figure 2 show the norm of dislocation density tensor α , which is found from the dislocation density vector ( ) k and Burger vector ( ) k b via Eq. (5). The x , y and z axes are taken to be along [110], [110] and [001], respectively, so that the dislocation loops are on a (111) slip plane. The Burgers vector is oriented along the y axis. A hybrid mesh with tetrahedral and pyramidal elements is used to replicate the FCC lattice structure, as described in (Xia and El-Azab, 2015). The mesh size here is mesh l = . A prescribed dislocation velocity v = 0.03 µm/ns is used to cause the two loops to expand and the time step t = . Periodic boundary conditions are used in this test example. Figure 2 Initial dislocation structure. Two dislocation loops are placed on a (111) slip plane. The evolution of the dislocation density, calculated by solving the transport equation (17), is depicted in Figure 3 in a view perpendicular to the slip plane. Aside from some diffusion, the evolution of the two dislocation loops looks reasonably accurate. At first, the two loops expand and then they interact with one another and begin to annihilate in places where opposing line age | 16 directions intersect. Finally, the two dislocation loops merge into one single dislocation loop that continues to expand. This simple result shows that the dislocation transport-reaction equation is solved accurately. In particular, by using the least squares finite element method and applying a divergence-free constraint on the dislocation density vector, the annihilation process is captured without the numerical dispersion that emerges whenever first-order hyperbolic partial differential equations are not carefully solved (Varadhan et al., 2006). It is to be noted here that, in this test problem, the dislocation transport equation (Eq. (6)) is solved using the same prescribed dislocation velocity for both algorithms. As such, the evolution of the dislocation field is the same, which facilitates the comparison of the plastic distortion fields resulting from the two update algorithms. Figure 3 Dislocation density evolution on the slip plane. The two loops expand and annihilate with each other to form a larger loop. As the dislocation density evolves, the associated plastic distortion will also change. As discussed in Section 3, two algorithms were used for comparison to update the plastic distortion. The plastic distortion calculated by direct time integration of the rate expression fixed by Orowan’s law is denoted as pI and the one found by our field dislocation mechanics approach is denoted as pII . The evolutions of these two quantities are shown in Figure 4 and Figure 5. (A Burgers vector age | 17 of 0.25 nm was used). The slipped area inside dislocation loops is shown in red, while the blue area shows unslipped area outside the dislocation loops. Figure 4 and Figure 5 can be compared with Figure 3 via the slipped area. The slipped area in Figure 3 is the area inside the dislocation loops. It can be seen in Figure 3 that the slipped area increases as the dislocation loops expand. After parts of the loops encounter each other and annihilate, the slipped area becomes connected. The slip area in Figure 3 is proportional to the plastic shear strain on the slip system considered in this simulation, and the latter gives the magnitude of the plastic distortion p . As such, the numerical algorithm that yields a plastic distortion distribution consistent with the slip area in Figure 3 is the more accurate one. Returning to Figure 4 and Figure 5, we notice that there are differences between the spatial distortion distributions of pI and pII . For example, at time steps 20 and 30 in Figure 4, fluctuations can be seen at the center where dislocations annihilate. It seems that the direct time integration algorithm accumulates error in calculating the plastic distortion in regions where dislocations of opposite direction annihilate. By using field dislocation mechanics (as shown in Figure 5), the evolution of the plastic distortion is similar, except that the fluctuations are highly reduced at the center. As shown later, the fluctuations in the field dislocation mechanics simulation only arise from the time integration of the compatible part, while the incompatible part is consistent with the current dislocation density tensor. This enhanced numerical stability has important implications for stress calculations. Figure 4 Evolution of the plastic distortion updated using the direct time integration algorithm. age | 18 Figure 5 Evolution of the plastic distortion, updated using the field dislocation mechanics algorithm. To see how accurately the plastic distortion is updated, we take the numerical curl of the plastic distortion field updated by the two algorithms, pI I = − and pII II = − , where I and II are the dislocation density tensors related to the corresponding plastic distortion field. Their evolutions are shown in Figure 6 and Figure 7. The dislocation density tensor , obtained from solving the transport equation, was already shown in Figure 3. Comparing I and II from Figure 6 and Figure 7, respectively, with the of Figure 3, we can see the differences in accuracy between the two algorithms. By using direct time integration, the dislocation density tensor I noticeably differs from . Namely, there are dislocations left in the annihilation region (Figure 6), which should not exist (Figure 3). This means that, although the dislocation evolution equation (17) is solved accurately in our model, the curl of the plastic distortion p may be inconsistent with the dislocation density tensor if the plastic distortion is updated by direct time integration of the rate expression derived from Orowan’s law. Numerical errors can accumulate during the numerical integration operation. Another difference between the dislocation density distributions of Figure 6 and Figure 3 is that there are two small dislocation loops with small density values in Figure 6 at which the initial loops are placed. Although an analytical expression of the initial plastic distortion pI is consistent with the analytical expression of the initial dislocation field , the age | 19 numerical curl operation from values at nodes of the mesh may lead to inconsistency in the initial values of the two fields. This inconsistency remains for the subsequent time steps. On the other hand, by using field dislocation mechanics, these problems are eliminated. The dislocation density tensor in Figure 7 is almost the same as in Figure 3. Since the incompatible part of pII is calculated directly from the current dislocation density tensor , numerical error is much smaller. Figure 6 The dislocation density tensor calculated from the plastic distortion updated by direct time integration. Figure 7 The dislocation density tensor calculated from plastic distortion, as updated by field dislocation mechanics. age | 20 In order to compare the accuracy of the update two algorithms quantitatively, it will be required to compare the errors in the plastic distortion obtained by the two algorithms relative to the plastic distortion computed by a third and a more accurate method. Generally speaking, such a third solution is not available. For the test examples presented here, however, that third method would be the analytical solution of the plastic distortion, which can in principle be found from knowledge of the density and the assumed velocity of dislocations. Such a comparison, however, would test together the accuracy of the solution of the transport equations together with the update algorithm, which is not the objective here. An alternate test of accuracy of the update algorithms is the error of the dislocation density computed as the curl of the plastic distortion relative to that which comes from the solution of the transport equation. Let us denote the latter by I and II for the two update algorithms. The relative error of interest, Error , can be defined in an integral sense over the domain of the solution by I(II)Error, I(II) || || d|| || d VV −= , (24) where is the volume of the simulation domain and || || is the norm operator for a second order tensor. The evolution of Error is shown in Figure 8. The figure shows that the direct time integration algorithm yields a much higher numerical error than that of the field dislocation mechanics approach. The rate of increase of Error is also higher in the case of direct time integration. Ideally speaking, the field dislocation mechanics approach should not accumulate errors. However, as the way the test is performed here is that we are comparing the gradients of the plastic distortion itself, which is obtained by numerical differentiation, i.e., by taking the numerical curl pI I = − and pII II = − , which introduce errors. age | 21 Figure 8 The relative error of dislocation density tensor (defined by Eq. (24)) calculated by the two algorithms. Although the dislocation velocity is prescribed and not calculated from the stress field in this test example, the stress fields I and II , obtained by using pI and pII as eigendistortions, can be compared to show how dislocation evolution would be affected if the stress field were used. The stress field is calculated using periodic boundary conditions, and the average stress is zero. The resolved shear stress values, I I = m s and II II = m s at time step 30 are shown in Figure 9. A Young’s modulus value of 112.5 GPa, a Poisson ratio of 0.34, and a Burgers vector of 0.25 nm were used in this stress calculation test. Figure 9 clearly shows that an inaccurate plastic distortion will lead to a significant error in the stress field. The spurious dislocations shown in Figure 6 cause the negative resolved shear stress at the center (left in Figure 9). The small dislocation loops in Figure 6 also perturb the stress field. On the other hand, the resolved shear stress calculated from field dislocation mechanics is more representative of the current dislocation state. It only has non-zero values around dislocation bundles, and it changes sign across the bundles. An accurate stress field is very important for showing correct patterns in bulk plasticity simulations at mesoscale, in which there are many dislocations and dislocation annihilation events occurring simultaneously. age | 22
Figure 9 The resolved shear stress, as alternately calculated from direct time integration (left) and field dislocation mechanics (right). As has been shown above, a direct time integration of the plastic distortion may lead to an inaccurate stress field. Accuracy can be improved by using field dislocation mechanics, in which the incompatible part of the plastic distortion is calculated directly from the current dislocation density tensor. However, the compatible part of the plastic distortion still needs time integration, see equation (23). A question that arises here is the extent to which the accuracy of this time integration affects the stress field. The stress field is not affected at all, if all mechanical boundaries are traction boundaries. Suppose the exact plastic distortion p has a compatible part z and an incompatible part , the compatible part of the plastic distortion pII calculated from field dislocation mechanics differs from the exact solution ( ) II = + z z z , while the incompatible part is exact II = . Then the corresponding displacement and stress fields satisfy the following two sets of equations, sym in : ( ) in on = = − +
0C u zn t = (25)
IIII II sym in : ( ) in on = = − − +
0C u z zn t = (26) If u is the solution of equation (25), it is easy to prove that II = + u u z (27) is the solution of equation (26). Substituting it back into equation (26), it is easy to obtain the following result, age | 23 II = . (28) Equations (27) and (28) imply that if the compatible part of the plastic distortion is calculated inaccurately, an error results in the displacement field, while the stress field is not affected. This fact explains why there are plastic distortion fluctuations at the center of Figure 5, but no stress fluctuation are present at the center of Figure 9 (right), which was generated from field dislocation mechanics calculations. All these errors derive from the calculation of the compatible part of the plastic deformation. So as long as an exact solution of the incompatible part of the plastic distortion is obtained, the stress field is exact. This fidelity can be guaranteed by using field dislocation mechanics, but not by using direct time integration. Also, although updating the compatible part of the plastic distortion involves time integration, this procedure incurs smaller errors than integrating the whole plastic distortion tensor over time, as per Figure 4 and Figure 5. It is noted that equations (27) and (28) hold only when all boundaries are traction boundaries. If a displacement boundary condition is applied, then the stress field calculated will be slightly different from the accurate value at that boundary, but it will remain accurate in the domain far from the displacement boundary. Bulk simulation with multiple slip systems
The purpose of this section is to present a larger test problem for the application of the update algorithms by running a 3D simulation with multiple slip in FCC crystal. In some simulations involving direct time integration, solution instability occurred after about 10 time steps, which is significantly before the target strain range is reached. That instability was manifested in the form of unbounded stress and dislocation density fluctuation in localized regions in the domain of solution. However, by using the field dislocation mechanics, the stress field is updated much more accurately and the solution remains stable well past a million time steps. As such, only the update algorithm based on field dislocation mechanics is used in the test presented in this section. The simulation domain used is a 5µm×5µm×5.303µm box, as shown in Figure 10(a). The x , y and z axes are selected to be along the [110], [110] and [001] crystallographic directions, respectively. The domain is discretized using a hybrid mesh with pyramidal and tetrahedral elements. The mesh size is chosen to be 62.5nm, so that the annihilation distance for an edge dislocation dipole is consistent with the property of the material (25 nm). Initial dislocation structures are set by placing dislocation loops in the domain with periodic boundaries, meaning age | 24 that if a dislocation exits the domain through one surface, it re-enters the domain from the opposite surface. In this simulation, 10 dislocation loops were placed randomly on each of the 12 slip systems of an FCC crystal, with radii ranging from 2µm to 6 µm, as shown in Figure 10(b). The initial total dislocation density over the domain is scaled to 1.5 µm -2 . Multiple dislocation interactions were considered. Cross slip was included with a coarse-grained rate found from discrete dislocation dynamics (DDD) simulations (Xia et al., 2016). Collinear annihilation and glissile junctions were considered in order to account for dislocation exchange between different slip systems. Lomer-Cottrell and Hirth locks were incorporated via a Taylor hardening term, with coefficients 0.084 (Lomer-Cottrell) and 0.051 (Hirth) (Madec, 2003). Figure 10 Initial dislocation structure. Each slip system has 10 dislocation loops placed randomly in the domain with radii ranging from 2 µm to 6 µm. (a) Total dislocation density; (b) Dislocation density on the first slip system. Material parameters were chosen to be those of steel (Déprés et al., 2004), as shown in Table 1. The crystal was loaded along the [001] direction at a strain rate of 20 s -1 . All boundaries are as periodic boundaries. An adaptive time step scheme was used, which was controlled by a Courant number max mesh c v t l = = , with max v being the maximum dislocation velocity in the domain. Typically, the calculated time increment for each step is about several nanoseconds. A relaxation step is performed prior to loading to release the stress field introduced by placing the initial dislocations into the simulation domain. It is performed by letting the initial dislocation microstructure evolve under zero applied stress until the maximum dislocation velocity in the is very close to zero, typically less than 10 -2 m/s. age | 25 Table 1 Material parameters used in the simulation (steel) (Déprés et al., 2004). Young’s modulus (GPa) Possion’s ratio Burgers vector (nm) Drag coefficient (Pa ∙ s) Lattice friction stress (MPa) 189 0.26 0.254 7.12×10 -6
3 The crystal is loaded to 1.5% strain. The stress-strain curve is shown in Figure 11. It starts in the elastic regime, which extends to 0.1% strain. A clear yield point of about 13 MPa is shown. After dislocations start to move, the crystal deformation progresses into the plastic regime. The hardening rate is roughly a constant after 1% strain. [001] loading is highly symmetric, with eight active slip systems. The reactions between these slip systems are the origin of the strain hardening. At 1.5% strain, the applied stress reaches about 17 MPa. The evolution of the dislocation density is shown in Figure 12. The total dislocation density increases with strain (Figure 12(a)). At 1.5% strain, the total dislocation density is about 8 µm -2 , about 5 times the initial dislocation density. There are several mechanisms accounting for dislocation multiplication in our model. Forest dislocations on other slip systems can act as pinning points. Therefore, according to the transport equation, dislocations will bow out, leading to an increase in dislocation line lengths. When double cross slip happens, dislocation density will increase (Hirth and Lothe, 1982; Hull and Bacon, 2011). Our model is able to handle such processes naturally, since the cross slip terms connecting every two collinear slip system pair are included in equation (17). Glissile junctions are another important source of dislocation multiplication, as the junction segment can bow out on its slip plane (Stricker et al., 2018; Stricker and Weygand, 2015) and ultimately detach. As depicted in Figure 12(b), the dislocation densities on the individual slip systems appear as two groups. The densities on the 8 active slip systems increase more or less similarly, with some fluctuations, while the dislocation densities on the 4 inactive slip systems increases at a significantly lower rate. From a strict crystal plasticity point of view, the density on the inactive systems should not increase at all. However, glissile junctions couple the “inactive” slip systems with the active ones, and cross slip also leads to density evolution. For example, slip system (111)[01 1] (numbered as slip system 1) and slip system (111)[101] (numbered as slip system 3) are active slip systems and they can form a glissile junction on slip system (111)[110] (numbered as slip system 11), which increases its dislocation density. Another reason for the finite nonzero increase in dislocation density on the inactive slip systems is that, although the Schmid factor is zero on theses slip systems, the local age | 26 resolved shear stress is not usually precisely zero everywhere, since it includes contributions not only from the external load, but also from the stress field arising from dislocations. This fact enables even those dislocations on the nominally inactive slip systems to move. Figure 11 Stress-strain curve. Figure 12 Evolution of the dislocation density. (a) Total dislocation density and (b) dislocation densities on individual slip systems. Figure 13 shows the dislocation microstructure evolution on a (111) slip plane. The location of the section is shown in Figure 13(a). Figure 13(b), (c) and (d) show the total dislocation density at different strains. The total dislocation density is defined as the sum of dislocation densities on all slip systems, ( ) total || || kk = , where ( ) k is the dislocation density vector on slip system k . The arcs in Figure 13(b) are the dislocation loops on that slip system, while the diffuse dots are dislocation loops on the other slip systems piercing the plane and serving as forest dislocations. age | 27 As the strain increases, dislocation density propagates over the domain through dislocation transport and dislocation reactions. Dislocations accumulate at some locations, forming a heterogeneous pattern. Figure 13(d) shows that dislocations are likely to accumulate in a pattern along three specific directions. One preferred direction is horizontal and the other two are oriented from the horizontal line. These directions are the intersections of the slip plane with the other three slip planes of the FCC crystal. These dislocation walls are not stationary. They evolve as the strain increases. The dislocation walls are at different locations at different strains. Figure 14 shows another cross section of the dislocation microstructure pattern. As it is not a slip plane, there are no arcs in the initial configuration and all dislocation loops pierce this plane, see Figure 14(b). The dislocation pattern in Figure 14 is different from that in Figure 13. Dislocation walls are formed at angles of from the horizontal line. They are shorter and less prominent than those in Figure 13. Figure 13 Dislocation microstructure evolution. (a) The location of the plane of view in the domain, which coincides with the (111) slip plane. (b) Initial total dislocation density. (c) age | 28 Total dislocation density at 0.6% strain. (d) Total dislocation density at 1.2% strain. The red lines show three directions along which dislocations accumulate the most. Figure 14 Dislocation microstructure evolution. (a) The location of the plane of view in the domain. (b) Initial total dislocation density. (c) Total dislocation density at 0.6% strain. (d) Total dislocation density at 1.2% strain. The red lines show three directions along which dislocations accumulate the most. Concluding remarks
Two algorithms for updating the plastic distortion tensor and internal stress tensor were presented and tested in this paper. The first algorithm performs a direct time integration of the plastic distortion rate computed from Orowan’s law, and the second utilizes the field dislocation mechanics formulation (Acharya and Roy, 2006; Brenner et al., 2014; Roy and Acharya, 2006, 2005). These two algorithms were introduced along with a full description of the vector-density based continuum dislocation dynamics approach, consisting of two sets of equations. One set of equations describes the crystal mechanics of the problem, with the kinematics and stress equilibrium parts, and the other set consists of the transport-reaction equations of continuum dislocation dynamics. The update algorithm based on field dislocation mechanics calculates the age | 29 incompatible part of the plastic distortion directly from the dislocation density field, regardless of the dislocation evolution history. While direct time integration of the plastic distortion rate accumulates numerical errors over time, the algorithm based on field dislocation mechanics was found to predict a plastic distortion field that is consistent with the instantaneous dislocation density field irrespective of time. However, this gain in accuracy incurs a computational cost. The use of field dislocation mechanics requires that one solves two more sets of equations for the compatible and incompatible parts of the plastic distortion. This method was explained. It consists of a staggered numerical scheme consisting of alternating steps of solving the stress equilibrium problem, followed by an update of the dislocation density effected by solving the dislocation transport-reaction equations. This sequence then is repeated. In this scheme, the first step yields the resolved shear stress required to compute the dislocation velocity, while the second step helps update the plastic strain required to solve the stress equilibrium problem. The stress field was solved by a standard Galerkin finite element method, and the dislocation dynamics calculations are performed using a least squares finite element method. A simple test example was performed that clearly shows that the update algorithm based on field dislocation mechanics predicts a more accurate plastic distortion distribution and stress field than the older time integration method. The relative influences of the accuracy of the compatible and incompatible parts of the plastic distortion on stress field fidelity also was discussed. The conclusion is that the errors developed in the stress field mainly derive from the manner in which the incompatible part of the plastic distortion is calculated, and these errors can be reduced using the improved update algorithm presented herein. This algorithm was used to simulate the behavior of FCC steel under uniaxial loading and multiple slip conditions. The expected trends for the mechanical behavior of a FCC metal at small strain were obtained. All known major dislocation multiplication mechanisms were included in the model. The evolution of the dislocation microstructures also was illustrated. In the continuum dislocation dynamics framework presented here, more mechanisms regarding dislocation reactions can be included to study the mechanical behavior of the material at microscale. For example, more accurate dislocation junction reaction rates can be included from DDD, and interactions of dislocation with point defects can be introduced. age | 30
Acknowledgements
The authors are grateful for the support from the Naval Nuclear Laboratory, operated by Fluor Marine Propulsion, LLC for the US Naval Reactors Program. A. El-Azab assisted with the formulation of field dislocation mechanics with support from the US Department of Energy, Office of Science, Division of Materials Sciences and Engineering, through award number DE-SC0017718. K. Starkey assisted with the deployment of the numerical approach and was supported by the National Science Foundation, Division of Civil, Mechanical, and Manufacturing Innovation (CMMI), through award number 1663311 at Purdue University.
Appendix: Finite element formulations of the equations
The Galerkin finite element method is commonly used to solve mechanical problems (Belytschko et al., 2013). By taking the weak form of equation (11) and discretizing it, a standard linear algebraic system is obtained for the nodal displacements, u u
K u P = , (A1) where T T11 21 31 1 2 3 , , ,..., , ,
I I I u u u u u u u = , with I being the number of nodes per element, is the collection of the nodal values of displacement in an element, and u K is the element stiffness matrix expressed in the form T d u K B C B V = , (A2) with C being the elastic stiffness matrix. B contains the derivatives of the shape functions of the element. u P is the load vector, which is contributed in the current case by the plastic distortion term in the stress equilibrium equation (11). Considering the two update algorithms explained in section 3, this load vector has the form: ( ) TT d with algorithm 1,d with algorithm 2, uu P B C E N VP B C B z N V = = − (A3) age | 31 where , z and are the corresponding node values of , z and . N and N are the matrices that contain shape functions of the element. E defines the slip normal and slip direction of each slip system. Equation (A1) is the linear algebraic system to be solved for displacement. The stiffness matrix u K of the discretized system remains unchanged at different time steps, so it is calculated and assembled only at the initial time step. The load vector u P changes due to the evolution of or z and , and it is updated at every step. Next, we turn attention to the dislocation transport problem. We use a 2D representation of the 3D dislocation vector field to reduce the number of degrees of freedom by considering a set of local coordinate systems with base vectors ˆ ˆ ˆ( , , ) e e e defined such that, for each slip system, ˆ e is in the direction of the Burgers vector, ˆ e is in the direction of the normal to the slip plane, and ˆ ˆ ˆ = e e e . With these relationships in hand, the vector and tensor quantities can be written in both global ( , , ) e e e and local ˆ ˆ ˆ( , , ) e e e coordinates. For example, ˆ ˆ ˆ ˆˆ and i i j j ij i j kl k l = = = = e e e e e e , (A4) with the components transforming using the matrix ˆ ij i j Q = e e according to ˆ ˆ and i ji j ij ki lj kl Q Q Q = = . (A5) In the local coordinate systems, however, the dislocation density and velocity on a given slip system have only two components, ˆ ˆ ˆ ˆ ˆ ˆˆ ˆ and v v = + = + e e v e e , (A6) so that the dislocation density evolution equation (6) can be written in its component form ˆ ˆ ˆ ˆ ˆ ˆˆ ˆ ˆ ˆ( ) and ( )ˆ ˆ v v v vx x = − = − − . (A7) The last equations can be rewritten in a matrix form, ˆ ˆˆ ˆ ˆˆ t A A A Ax x = + + , (A8) where age | 32 ˆ ˆˆ ˆ 0 01 0 ˆ ˆ, , , .ˆ ˆ0 1 0 0ˆ ˆˆ ˆ t v vx x v vA A A Av vv vx x − − = = = = − − (A9) Equation (A8) is discretized in time by a first order implicit Euler method, giving
11 10 1 21 2 2 2 ˆ ˆˆ ˆ ˆ ˆ n nt t
A t A t A t A Ax x + − + + + = − , (A10) with n + and n representing the unknown variable at time steps n +1 and n , respectively. The velocity in A , A and A is taken at time step n . If we write the unknown variable in terms of the nodal values ˆ Ii with the aid of shape functions of a finite element, which denote the i th component of I th node, ˆ ˆ ˆˆ ˆ ˆˆ ˆ ˆˆ ˆ ˆˆ ˆ ˆ, , ˆ ˆˆ ˆ ˆˆ ˆ ˆˆ ˆ ˆ M M MM M M
N B Bx x = = = (A11) where M is the number of nodes of the element and M MM MM M
N N NN N N NN N Nx x xB N N Nx x xN N Nx x xB N N Nx x x = = = (A12) For ease of reading, the superscript has been dropped in the nodal quantities. Eq. (A10) can be rewritten for the node values as follows: age | 33 ˆ ˆ n n L W + = , (A13) with ,,ˆ ˆ ˆ ˆ ˆ ˆ ˆ, , , , , , . tt I I L A N t A N t A B t A BW A N = − + + + = − = (A14) The system of transport equations for dislocation density evolution is hyperbolic. Some level of numerical diffusion and dispersion is always anticipated with the discretization of this kind of partial differential equations (Varadhan et al., 2006). A least squares finite element method (LSFEM) (Jiang, 2013) is used in our model to stabilize dispersion. The least squares residual corresponding to equation (A13) is given by ( ) ( )
Tn n n n
R L W L W + + = − − . (A15) A key attribute of dislocations is that they do not end inside a crystal. In continuum dislocation theory, this requirement is enforced by imposing the condition that the dislocation tensor is divergence free. That is, = . (A16) Theoretically, equation (A16) is satisfied in our model because of equation (2). However, in the numerical implementation, a divergence-free condition must be used as a constraint on the evolution of the field to ensure numerical accuracy. An additional term representing this constraint is thus added to the residual in equation (A15), which has the form cp R ch = − , (A17) where c is a control parameter, h is the mesh size, and cp is the coupling term due to dislocation reactions and cross slip, as defined in Eq. (17). Including cp in Eq. (A17) is based on the fact that when dislocation reaction happens, the dislocation line can not be considered as continuous and connected on a individual slip system anymore unless we take the reaction portion back. In a matrix form, the last residual becomes ( )
212 cp nd R ch B + = − , (A18) age | 34 where I Id
N N N N N NB x x x x x x = . (A19) Adding equation (A18) into equation (A15) and taking the first variation of the function to be zero, we have ( ) ( )
11 1 2 cp ˆ ˆ ˆ ˆ ˆ+ d 0
T nT Tn n n d d
L L W ch B B ++ + − − = . (A20) As n + is arbitrary, Eq. (A20) can be rewritten into the form ( ) ( )
12 12 cp ˆdˆ ˆ d
T T nd d nT Tn d d
L L ch B BL W ch B B + + + = + , (A21) which yields the algebraic system ˆ n K P + = , (A22) with ( ) ( ) d , ˆ ˆ d T Td d nT Tn d d
K L L ch B BP L W ch B B + = + = + . (A23) The solution of Eq. (A22) yields the dislocation density evolution. The stiffness matrix K and the load vector P depend on the problem information at time step n and are updated at every time step. When we use field dislocation mechanics to calculate the incompatible and compatible parts of the plastic distortion, Eqs. (22) and (23) are also solved by finite element methods. For the incompatible part of the plastic distortion, Eq. (22) can be written in a matrix form, age | 35 x xx xx xx x x − − = − . (A24) Using L to represent the linear operator, the equation can be expressed in the form, L = . (A25) By using the least squares finite element method to solve this partial differential equation, we finally reach an algebraic system of the form, K P = , (A26) where the stiffness matrix K and load vector P are given by TT dd K L LP L = = . (A27) For the compatible part of the plastic distortion, Eq. (23) can be rewritten for a time step t in the form ( ) ( ) ( ) ( ) ( ) ( ) n n n n t + + − + = z z v . (A28) A standard Galerkin finite element method is used to solve this partial differential equation, the weak form of which is ( ) ( ) ( ) ( ) ( ) ( ) ( ) d 0 n n n n n t + + + − + = z z v z . (A29) The corresponding algebraic system has the form nz z K z P + = , (A30) where the stiffness matrix z K and load vector z P are given by age | 36 ( ) TT d d z nz K B BP B B z p = = + . (A31) The B matrix is the same as in equation (A2) and p denotes ( ) ( ) ( ) n n t + v . References
Acharya, A., 2004. Constitutive analysis of finite deformation field dislocation mechanics. J. Mech. Phys. Solids 52, 301–316.
Acharya, A., Roy, A., 2006. Size effects and idealized dislocation microstructure at small scales:
Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part I. J. Mech. Phys. Solids 54, 1687–1710. Arsenlis, A., Parks, D.M., Becker, R., Bulatov, V. V., 2004. On the evolution of crystallographic dislocation density in non-homogeneously deforming crystals. J. Mech. Phys. Solids 52, 1213–1246. Balluffi, R.W., 2016. Introduction to Elasticity Theory for Crystal Defects. World Scientific Publishing Company. Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K., 2013. Nonlinear finite elements for continua and structures. John wiley & sons. Brenner, R., Beaudoin, A.J., Suquet, P., Acharya, A., 2014. Numerical implementation of static Field Dislocation Mechanics theory for periodic media. Philos. Mag. 94, 1764–1787. Cottura, M., Appolaire, B., Finel, A., Le Bouar, Y., 2016. Coupling the Phase Field Method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys. J. Mech. Phys. Solids 94, 473–489. Cui, Y., Po, G., Pellegrini, Y.P., Lazar, M., Ghoniem, N., 2019. Computational 3-dimensional dislocation elastodynamics. J. Mech. Phys. Solids 126, 20–51. Cui, Y.N., Lin, P., Liu, Z.L., Zhuang, Z., 2014. Theoretical and numerical investigations of single arm dislocation source controlled plastic flow in FCC micropillars. Int. J. Plast. 55, 279–292. age | 37
Deng, J., El-Azab, A., 2010. Temporal statistics and coarse graining of dislocation ensembles. Philos. Mag. 90, 3651–3678. Déprés, C., Robertson *, C.F., Fivel, M.C., 2004. Low-strain fatigue in AISI 316L steel surface grains: a three-dimensional discrete dislocation dynamics modelling of the early cycles I. Dislocation microstructures and mechanical behaviour. Philos. Mag. 84, 2257–2275. Devincre, B., Hoc, T., Kubin, L., 2008. Dislocation Mean Free Paths and Strain Hardening of Crystals. Science (80-. ). 320, 1745–1748. Devincre, B., Kubin, L., Hoc, T., 2006. Physical analyses of crystal plasticity by DD simulations. Scr. Mater. 54, 741–746. El-Awady, J.A., 2015. Unravelling the physics of size-dependent dislocation-mediated plasticity. Nat. Commun. 6. El-Azab, A., 2000. Statistical mechanics treatment of the evolution of dislocation distributions in single crystals. Phys. Rev. B - Condens. Matter Mater. Phys. 61, 11956–11966. Evers, L.P., Brekelmans, W.A.M., Geers, M.G.D., 2004. Non-local crystal plasticity model with intrinsic SSD and GND effects. J. Mech. Phys. Solids 52, 2379–2401. Franciosi, P., Berveiller, M., Zaoui, A., 1980. Latent hardening in copper and aluminium single crystals. Acta Metall. 28, 273–283. Grilli, N., Janssens, K.G.F., Nellessen, J., Sandlöbes, S., Raabe, D., 2018. Multiple slip dislocation patterning in a dislocation-based crystal plasticity finite element method. Int. J. Plast. 100, 104–121. Groma, I., 1997. Link between the microscopic and mesoscopic length-scale description of the collective behavior of dislocations. Phys. Rev. B - Condens. Matter Mater. Phys. 56, 5807–
Groma, I., Csikor, F.F., Zaiser, M., 2003. Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics. Acta Mater. 51, 1271–1281. Gurtin, M.E., Anand, L., Lele, S.P., 2007. Gradient single-crystal plasticity with free energy dependent on dislocation densities. J. Mech. Phys. Solids 55, 1853–1878. Hirth, J., Lothe, J., 1982. Theory of Dislocations. John Wiley & Sons, New York.
Hochrainer, T., 2016. Thermodynamically consistent continuum dislocation dynamics. J. Mech.
Phys. Solids 88, 12–22. age | 38
Hochrainer, T., 2015. Multipole expansion of continuum dislocations dynamics in terms of alignment tensors. Philos. Mag. 95, 1321–1367. Hochrainer, T., Sandfeld, S., Zaiser, M., Gumbsch, P., 2014. Continuum dislocation dynamics: Towards a physical theory of crystal plasticity. J. Mech. Phys. Solids 63, 167–178. Hochrainer, T., Zaiser, M., Gumbsch, P., 2007. A three-dimensional continuum theory of dislocation systems: Kinematics and mean-field formulation. Philos. Mag. 87, 1261–1282. Hull, D., Bacon, D.J., 2011. Introduction to dislocations. Butterworth-Heinemann. Jiang, B., 2013. The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics. Springer Science & Business Media. Kooiman, M., Hütter, M., Geers, M.G.D., 2014. Collective behaviour of dislocations in a finite medium. J. Stat. Mech. Theory Exp. 2014, P04028. Kosevich, A.M., 1965. DYNAMICAL THEORY OF DISLOCATIONS. Sov. Phys. Uspekhi 7, 837–854. Kröner, E., 1958. Kontinuums theorie der Versetzungen und Eigenspannungen. Springer. Kubin, L., Devincre, B., Hoc, T., 2008. Toward a physical model for strain hardening in fcc crystals. Mater. Sci. Eng. A 483–484, 19–24. Leung, H.S., Leung, P.S.S., Cheng, B., Ngan, A.H.W., 2015. A new dislocation-density-function dynamics scheme for computational crystal plasticity by explicit consideration of dislocation elastic interactions. Int. J. Plast. 67, 1–25. Lin, P., El-Azab, A., 2020. Implementation of annihilation and junction reactions in vector density-based continuum dislocation dynamics. Model. Simul. Mater. Sci. Eng. 28, 045003. Lin, P., Liu, Z., Cui, Y., Zhuang, Z., 2015. A stochastic crystal plasticity model with size- dependent and intermittent strain bursts characteristics at micron scale. Int. J. Solids Struct. based on coupled glide-climb model. Int. J. Plast. 92, 2–18.
Liu, Z.L., Liu, X.M., Zhuang, Z., You, X.C., 2009. A multi-scale computational model of crystal plasticity at submicron-to-nanometer scales. Int. J. Plast. 25, 1436–1455. age | 39
Madec, R., 2003. The Role of Collinear Interaction in Dislocation-Induced Hardening. Science (80-. ). 301, 1879–1882. Mishra, A., Alankar, A., 2019. Revisiting dislocation reactions and their role in uniaxial deformation of copper single crystal micro-pillars. Model. Simul. Mater. Sci. Eng. 27, 055010. Monavari, M., Sandfeld, S., Zaiser, M., 2016. Continuum representation of systems of dislocation lines: A general method for deriving closed-form evolution equations. J. Mech. Phys. Solids 95, 575–601. Monavari, M., Zaiser, M., 2018. Annihilation and sources in continuum dislocation dynamics. Mater. Theory 2, 3. Mura, T., 1963. Continuous distribution of moving dislocations. Philos. Mag. 89, 843–857. Nye, J.F., 1953. Some geometrical relations in dislocated crystals. Acta Metall. 1, 153–162. Peach, M., Koehler, J.S., 1950. The Forces Exerted on Dislocations and the Stress Fields Produced by Them. Phys. Rev. 80, 436–439. Roy, A., Acharya, A., 2005. Finite element approximation of field dislocation mechanics. J. Mech. Phys. Solids 53, 143–170. Roy, A., Acharya, A., 2006. Size effects and idealized dislocation microstructure at small scales: Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part II. J. Mech. Phys. Solids 54, 1711–1743. Schulz, K., Wagner, L., Wieners, C., 2019. A mesoscale continuum approach of dislocation dynamics and the approximation by a Runge-Kutta discontinuous Galerkin method. Int. J. Plast. 120, 248–261.
Sandfeld, S., Hochrainer, T., Zaiser, M., Gumbsch, P., 2011. Continuum modeling of dislocation plasticity: Theory, numerical implementation, and validation by discrete dislocation simulations. J. Mater. Res. 26, 623–632. Sandfeld, S., Thawinan, E., Wieners, C., 2015. A link between microstructure evolution and macroscopic response in elasto-plasticity: Formulation and numerical approximation of the higher-dimensional continuum dislocation dynamics theory. Int. J. Plast. 72, 1–20.
Sandfeld, S., Zaiser, M., 2015. Pattern formation in a minimal model of continuum dislocation plasticity. Model. Simul. Mater. Sci. Eng. 23, 065005. age | 40
Sedláek, R., Schwarz, C., Kratochvíl, J., Werner, E., 2007. Continuum theory of evolving dislocation fields. Philos. Mag. 87, 1225–1260. Sills, R.B., Bertin, N., Aghaei, A., Cai, W., 2018. Dislocation Networks and the Microstructural Origin of Strain Hardening. Phys. Rev. Lett. 121, 85501. Stricker, M., Sudmanns, M., Schulz, K., Hochrainer, T., Weygand, D., 2018. Dislocation multiplication in stage II deformation of fcc multi-slip single crystals. J. Mech. Phys. Solids 119, 319–333. Stricker, M., Weygand, D., 2015. Dislocation multiplication mechanisms - Glissile junctions and their role on the plastic deformation at the microscale. Acta Mater. 99, 130–139. Sudmanns, M., Stricker, M., Weygand, D., Hochrainer, T., Schulz, K., 2019. Dislocation multiplication by cross-slip and glissile reaction in a dislocation based continuum formulation of crystal plasticity. J. Mech. Phys. Solids 132, 103695. Varadhan, S.N., Beaudoin, A.J., Acharya, A., Fressengeas, C., 2006. Dislocation transport using an explicit Galerkin/least-squares formulation. Model. Simul. Mater. Sci. Eng. 14, 1245–1270. Xia, S., Belak, J., El-Azab, A., 2016. The discrete-continuum connection in dislocation dynamics: I. Time coarse graining of cross slip. Model. Simul. Mater. Sci. Eng. 24, 075007. Xia, S., El-Azab, A., 2015. Computational modelling of mesoscale dislocation patterning and plastic deformation of single crystals. Model. Simul. Mater. Sci. Eng. 23, 055009. Yefimov, S., Groma, I., Van der Giessen, E., 2004. A comparison of a statistical-mechanics based plasticity model with discrete dislocation plasticity calculations. J. Mech. Phys. Solids 52, 279–300.