Driving forces on dislocations: finite element analysis in the context of the non-singular dislocation theory

This work presents a regularized eigenstrain formulation around the slip plane of dislocations and the resultant non-singular solutions for various dislocation configurations. Moreover, we derive the generalized Eshelby stress tensor of the configurational force theory in the context of the proposed dislocation model. Based on the non-singular finite element solutions and the generalized configurational force formulation, we calculate the driving force on dislocations of various configurations, including single edge/screw dislocation, dislocation loop, interaction between a vacancy dislocation loop and an edge dislocation, as well as a dislocation cluster. The non-singular solutions and the driving force results are well benchmarked for different cases. The proposed formulation and the numerical scheme can be applied to any general dislocation configuration with complex geometry and loading conditions.


Introduction
Dislocations are one of the most influential defect types in crystalline materials. For instance, plasticity and hardening in crystalline materials are caused by the collective movement of dislocations, making the influence of dislocations essential for the mechanical properties [1,2]. Dislocations can interact with each other or with other types of defects, such as point defects and grain boundaries, and even functional structures like ferroelectric domain walls [3,4]. In this way, they are also expected to have a large impact on the functionalities on the mesoscopic scale. To understand the influential mechanisms of dislocations and hence to enable new routes of materials design, a reliable and flexible approach to evaluate the elastic fields induced by dislocations and to determine the driving force they experience under circumstance of external loading or other defects is indispensable.
In the field of materials science, the Peach-Koehler formulation [5] has been widely applied, particularly successful in the discrete dislocation dynamics [2,[6][7][8]. Given the stress field on dislocation line, e.g., by using the analytical infinite-domain solution [6,7] or the finite element solution [2,8], the driving force per line length can be determined through the Peach-Koehler formula. Though for simple situations with limited number of dislocations and of simple dislocation types the stress field can be determined, the accurate determination of the fields for general complex cases is non-trivial. Moreover, this formulation regards only the mechanical contribution to the driving force. In fact, other contributions to the potential drop and thus to the driving force on a defect, e.g., electrostatic driving force in the case charged dislocations, are also important. For this purpose, the Peach-Koehler formulation should be extended, such as the electroelastic Peach-Koehler force proposed by Agiasotou and Lazar [9].
The theory of configurational forces has been a valuable tool to study the driving force on defects in materials, and has become one important branch of the modern micromechanics. It can be traced back to the work in 1951 by Eshelby [10]. The second-order Eshelby stress tensor was derived from the negative gradient of the total potential with respect to the position of the defect. Then the configurational forces can be formulated as the divergence of the Eshelby stress tensor [11]. This concept demonstrates explicitly that, if a defect could move for a unit length from its current position, the driving force on the defect is in principle the maximum of the possible potential drop of the system. It also explains the alternative name used in the micromechanics theory, i.e., the configurational force. Since it stems from the variation of the system potential, the configurational force theory provides a general and efficient way to analyze driving forces on the defects in multi-physics scenarios. Alternatively, Lazar and Kirchner [12] derived the Eshelby stress tensor by making use of the concept of the dislocation density tensor, which is further extended for dislocated piezoelectric materials [9].
On the basis of the configurational force theory, driving forces on defects can be simply evaluated numerically using the calculated elastic fields. Mueller and Gross [11,[13][14][15] developed a general computational concept to calculate the driving forces on the defects in anisotropic materials, including point defects, dislocations, inclusions, material interfaces and so on. In numerical implementation, the configurational force balance is implemented into the finite element scheme to calculate the configurational forces on nodes. Then the driving force on the defect can be calculated through summing up the nodal configurational forces around the defect (in accordance with the one determined by the Peach-Koehler formula). The model has been successfully used to analyze the driving forces on point defects, material interfaces and cracks. This numerical concept on the basis of the configurational force theory allows studies on defects in problems involving complex external loads and boundary conditions [16] or for the coupled electro-elastic problems [17,18].
Nevertheless, there are few efforts applying the configurational force theory and the related numerical approach to study the driving force on dislocations. Baxevanakis and Giannakopoulos [16] presented a twodimensional finite element method to calculate the driving force on the edge dislocation. The edge dislocation is modeled using a thermal analogue combined with a non-singular variable core model developed by Lubarda and Markenscoff [19]. For the thermal analogue, the eigenstrain of the dislocation is treated as thermal expansion; therefore, it only contributes to the volume strain of the system. The model is used to study the interaction between the edge dislocation, crack, and inclusion. The driving force on the dislocation is calculated through the J-integral and compared with the Peach-Koehler force. However, the thermal analogue model can only be used to predict the climb force on the edge dislocation, because it does not include a shear eigenstrain component and the applied shear stress has no effect on the dislocation, as discussed by Kolednik et al. [20]. Kolednik et al. also tried to use the eigenstrain strip model to calculate the elastic fields and driving forces of dislocations, where the eigenstrain is distributed in a narrow rectangular strip. Although the driving force on the dislocation is correctly predicted, the prediction of the displacement field does not agree well with the analytical solution. Therefore, they proposed a "cut-displace-glue"procedure operating on the mesh to model a two-dimensional model for the edge dislocation. The model correctly predicted the elastic fields of the dislocation, as well as the driving forces on the dislocation. However, the "cut-displace-glue" procedure operating on the mesh makes it difficult to model curved dislocation lines or other complex dislocation configurations. Besides, none of these works discussed about the distribution of nodal configurational forces around the dislocation core due to complex elastic fields of dislocations.
In order to remove the singularity in the classical solutions of dislocations, one can consider the non-singular dislocation theory developed based on the non-local elasticity [21][22][23] or the gradient elasticity [24,25]. The non-singular dislocation theory based on the gradient elasticity introduces an extra characteristic parameter which could be estimated using atomistic calculations [26]. The non-singular dislocation theory based on non-local elasticity have been often used in finite element simulations [8,16]. Cai et al. [22] proposed a nonsingular continuum theory of dislocations. In particular, they spread the Burgers vector isotropically around every point on the dislocation line and obtained the stress field of dislocations with finite stress at the dislocation core. The spreading function is characterized by one single parameter namely the dislocation core width. In use of the Green function theory, the obtained non-singular solution can be supposed to obtain the fields of dislocation lines containing various segments. On the other hand, various eigenstrain distributions were used to calculate the non-singular stress fields of dislocations in the finite element simulation, such as the thermal analogue model and eigenstrain strip model discussed above [16,20]. Recently, the spreading function in Cai's non-singular theory was used to define the eigenstrain distribution, correctly predicting the stress fields of screw and edge dislocations [8]. Whether the spreading function can also be applied to derive the driving force on dislocations, remains an open question. The configurational force theory in the context of the non-singular dislocation theory based on the spreading of the Burgers vector should be revisited.
In this work, we first formulate explicitly the eigenstrain on the basis of the spreading function of the Burgers vector, and contrast the resultant numerical elastic fields of dislocations with those analytical results of Cai et al. [22] and the classical singular solutions. Then based on the generalized eigenstrain theory of dislocations [27], we derived the Eshelby stress tensor through the gradient of the potential energy of dislocated elastic bodies for both the singular and non-singular dislocation theory. From the obtained Eshelby stress tensor and nonsingular finite element solutions, we numerically evaluated the driving forces on dislocations. Benchmarks are obtained for an edge dislocation, a screw dislocation, a dislocation loop and the interaction between a vacancy loop and an edge dislocation. Both the elastic fields and the driving forces of the dislocations are verified. To demonstrate the flexibility of the model and its application potential for complex situations, a cluster of edge dislocations is simulated. We obtained both the critical passing stress and the driving force on each individual dislocation, which are otherwise difficult to calculate.

Singular representation
The configurational forces of dislocations are derived in the following using the basic ideas of Eshelby [10]. Consider a dislocated body that occupies a region V with the boundary S. Assuming the small strain approximation and quasi-static problem, the dislocated solid is governed by the mechanical equilibrium where () , j = ∂/∂ x j , σ i j is the stress tensor, and the body force is ignored. For the general anisotropic case, the constitutive relation for the stress σ i j and elastic strain ε e i j is given as where c i jkl is the stiffness tensor, ε e kl is the elastic strain, ε D kl is the eigenstrain of dislocations, and the total strain ε kl is defined as the symmetric part of the displacement gradient u k,l as Either displacement or traction boundary condition can be considered whereū i is the prescribed displacement on the boundary, n j is the normal vector of the boundary, and t i is the applied traction on the boundary. The associated eigenstrain ε D i j results from the plastic slip generated during the glide of the dislocation line. The eigenstrain tensor is defined as the symmetric form of a dyadic product of the Burgers vector b i and the normal vector n i of the slip plane D where δ(x − D) is the one-dimensional Dirac delta function in the normal direction of the slip plane D [27]. Index-notation and Einstein's summation convention for repeated indices are used here and in the following. Eq. (5) ensures that the eigenstrain is located merely on the slip plane, i.e., when x ∈ D. Consider the definition of the one-dimensional Dirac delta function [27] where x D is the location vector of the points on the slip plane D. Substituting above equation into Eq. (5), the eigenstrain can be rewritten in terms of convolution as where i j = b i n j + b j n i /2, and the symbol * means the convolution over the slip plane.
In the eigenstrain theory, by introducing the concept of dislocation density, the Burgers vector can be defined as [9] where C is the Burgers circuit, S is the Burgers surface bounded by C. β D i j = b i n j * δ(x) is the plastic distortion due to the dislocation, β i j = u i, j − β D i j is the elastic distortion, and u i, j is the total distortion. α i j is the dislocation density tensor, according to the Green theorem, it is defined as where jkl is the permutation tensor. The strain energy density for the dislocated body can be rewritten as where the symmetry of the stress tensor is used in the second equation. In the following, a Volterra dislocation with constant Burgers vector is considered. The gradient of the strain energy is given as By multiplying Eq. (9) with the permutation tensor k jl , one can reformulate the definition of the dislocation density tensor as Then, consider the mechanical equilibrium, the gradient of the strain energy is reformulated as Rearranging above equation yields the configurational force balance in the well-known local form where the generalized Eshelby stress tensor k j takes the following form and the generalized configurational body forces g k is given as The driving force per unit length on the dislocation can then be calculated by the surface integration of the divergence of the Eshelby stress tensor or the configurational body force. As an example, consider a straight dislocation line with dislocation density of where x C i is the location of the dislocation line,ξ l is the sense vector of the dislocation line, and σ 0 i j = σ i j x i =x C i is the stress field at dislocation core due to external loading. For more details regarding the physical interpretation of the Eshelby stress tensor and Peach-Koehler force, one can refer to Ref. [28,29].

Non-singular representation
The classical solution of dislocations discussed above has singularity at the core center. However, the field singularity at the dislocation core has no sound physical interpretation and leads to inconvenience in the numerical simulations. One general idea is to spread the eigenstrain distribution on the slip plane over the surrounding region (regularization of dislocation slip) [2]. The regularization of the eigenstrain distribution is crucial to the accuracy of the numerical solution of dislocations. Cai et al. [22] proposed a non-singular continuum representation of the Burgers vector, which ensures that outside the pre-defined dislocation core the derived stress fields agree well with the classic solutions.
Though the results in this subsection are not limited to it, we demonstrate the configurational force theory in the context of the non-singular dislocation theory. For an arbitrary point on the slip plane, spreading the Burgers vector around the point in the three-dimensional space as where f i (x) is the Burgers vector density function. The convolution of the spreading functionw(x) with itself defines the second spreading function where r = ||x||. The second spreading function w(x) is used to modify the distance function R = Then the dislocation core width parameter h is introduced in the formulation such that the modified distance function R h ensures the resulting non-singular solution closely resemble the singular solution of the classical theory outside the dislocation core. The analytical form ofw(x) is unknown, and Cai et al. [22] presented the following approximate expressioñ where h 1 = 0.9038h, h 2 = 0.5451h, and m = 0.6575 are constants.
In the framework of the non-singular continuum theory of dislocations, after spreading the Burgers vector around the slip, the Dirac delta function shown in the last subsection is replaced now by the spreading functioñ Comparison of Eq. (7) and Eq. (21) implies that the eigenstrain takes the similar form in both the singular and the non-singular theory. In the meantime, the first two equations in Eq. (8) are satisfied with the distribution of Burgers vector Eq. (18) in the non-singular dislocation theory. Therefore, following Eq. (7) to Eq. (16), the configurational force balance for the non-singular continuum theory can be derived by simply replacing The resulting configurational force balance takes the same form as Eq. (14), and the non-singular Eshelby stress tensor in the context of the non-singular continuum theory of dislocations is given as where β ns ik = u i,k − b i n k * w(x). From Eq. (22) it can be seen that the gradient of the distribution function does not show explicitly in the final form of the non-singular Eshelby stress tensor. Eq. (18) guarantees the magnitude of the regularized Burgers vector in the non-singular theory is the same as that of the singular theory, and Eq. (21) guarantees that the eigenstrain is obtained through the convolution over the slip plane similar to Eq. (7). Hence, the derivation of the non-singular Eshelby stress tensor (22) is mathematically consistent with the derivation of the singular Eshelby stress tensor (15) in the context of the non-singular continuum theory of dislocations developed by Cai et al. [22].

Numerical calculation of the regularized eigenstrain
In the conventional finite element simulation, the strain is calculated on integration points rather than on the nodes. The distribution of the eigenstrain on integration points is crucial to the accuracy of the numerical solution of dislocations. In the following, we consider the eigenstrain distribution calculated by using the concept of distributed Burgers vector [22], as introduced in Eq. (18). We demonstrate the calculation by the basic examples of the straight and circular dislocation line. For these dislocations, the normal vector will be a constant in local coordinate system. Therefore, the distribution of eigenstrain in Eq. (21) can be written as where The function W (x, h) defines the distribution of eigenstrain around the slip plane in the non-singular continuum theory of dislocations. In order to reduce the computational cost, the integration is truncated at a distance r c such that the distribution function becomes and where H (r c −r ) is the Heaviside step function. A small r c leads to under estimation of the eigenstrain. r c = 2h is used in this work and it leads to less than 5% error compared to the analytical solution (r c → ∞) [8].
The two-dimensional setup of the numerical integration ofw(r, h) is shown in Fig. 1a. The dislocation is located at x = {0, 0, 0} T and the dashed line represents the slip plane. It can be applied to the straight edge dislocation or straight screw dislocation with the dislocation line lying in the x 3 direction. The distribution of the Burgers vector along the x 3 direction remains constant. Therefore, the range of the integration domain for x D 3 is from −r c to r c . Then the surface integration of w(r, h) can be represented as where L represents the integration path along the x 1 direction on the slip plane, and H (r c − r ) is removed from the integration by changing the integration domain to be always inside −r c to r c around the integration point x. The analytical solution of the integration along the x D 3 -axis can be determined as where The integration along the x D 1 direction is done numerically, and the trapezoidal numerical integration is used in this work. In other words, where N is the number of integration points (N = 20 is used in this work), and X i is the X calculated at the ith integration point. Forx 1 ∈ [0, r c ], the slip plane is assumed to be extended outside the sample such that the integration domain is is only a function of the coordinate x i , therefore, it only has to be calculated once at the initial step. Figure 1b shows the numerical results of W (x, h), which also represents the distribution of the eigenstrain. The eigenstrain is constant along the slip plane, except for the region on the right side of the dislocation core. For illustration, only a circular dislocation loop will be considered in this work. The distribution function W (x, h) for a circular loop is axis-symmetric. It means that the analytical solution Eq. (28) can also be obtained in the circumferential direction, and the numerical integration along the radius direction is the same with a straight dislocation line. Therefore, the distribution function W (x, h) of the straight edge dislocation can also be used along the radius direction for the circular dislocation loop. For the more general application of the distribution of eigenstrain, readers can refer to Ref. [8].

Finite element formulation
Based on the finite element method, the weak form formulation is constructed for numerical implementation. Multiplying the mechanical equilibrium (1) with the test function η i and integration over the domain V V σ i j, j η i dV = 0.
Integrating by parts leads to the weak form formulation as According to the conventional linear finite element method, the test function and its gradient can be interpolated through the shape function N I as where η I i are the nodal values of the test function. Evaluating Eq. (31) over each element yields where V e is the volume of an element.
Similarly, multiplying the configurational force balance (14) with the test function η i and integrating over the domain V , the weak form formulation is obtained through integrating by parts as The weak form formulation of mechanical equilibrium is used for calculation of the stress σ i j and the displacement field u i of the system. Then the numerical solution of the elastic fields is used in the weak form formulation of the configurational force balance to obtain the nodal configurational forces. Since the elastic fields are calculated based on the non-singular continuum theory of dislocations, according to Eq. (34), the nodal configurational forces acting on a node I are calculated as where the assembly operation is performed over all n el elements adjacent to node I . Then the driving forces on the dislocation can be calculated by summing up all nodal configurational forces around the dislocation core as where n nod is the total number of nodes in the chosen integration volume. The user element for dislocations is developed and the finite element simulations are performed using the open source software MOOSE [30]. The numerical results are compared with the analytical solutions of dislocations and also benchmarked by the Peach-Koehler formula [1]. In the fifth example of a dislocation cluster we demonstrate the flexibility and the application potential of the proposed numerical scheme for determining the driving force in a complex scenario.

Edge dislocation
We start with a straight edge dislocation with the Burgers vector b = [b 0 0 0] with b 0 = 1, and the normal vector of the slip plane n = [0 1 0]. The edge dislocation is placed in the center of the sample. The straight edge dislocation is modeled as a plane strain problem. To assure that there are enough integration points for calculating the distribution of eigenstrain, it is suggested that two element should be covered in r c range. The minimum h used in this work is 1, thus, we choose the minimum element size of 0.5. The sample size is 200 × 200 and the mesh size is 400 × 400. In order to compare the numerical solution of a finite size with the analytical solution of a single dislocation in an infinite medium, one can either use a large-enough sample or apply asymptotically the traction t i = σ a i j n j on the surfaces based on the analytical stress solutions σ a i j . Since the computation is not very expensive, we apply the traction free boundary condition on a sufficiently large sample to reduce traction free boundary effects in this work. We studied the influence of the sample size on the relative error of the numerical solution compared to the analytical solution. For a square sample with a width of 50, the relative error of σ 11 at x 2 = ±10 is around −39%. When the width of the sample Fig. 2 Benchmarks for the edge dislocation. a-c The numerical and analytical stress fields (dimensionless) agree well outside the core region. d The displacement jump (dimensionless) is correctly predicted near the slip plane reaches 150, the relative error reduces to around −4%. Therefore, the sample size used in this numerical example is large enough to get good approximation of the stress field of dislocations near the dislocation core. The stress distribution near the dislocation core are shown in Fig. 2(a-c). Outside the dislocation core x i ∈ [−h, h], the numerical solution of the stress fields agrees well with both the singular analytical solution and Cai's non-singular solution [22]. Inside the dislocation core region, the numerical solution lies between the singular analytical solution and Cai's non-singular solution. The displacement u 1 on the left surface are shown in Fig. 2d. The displacement jump near the slip plane is equal to the Burgers vector, which indicates the validity of the numerical solutions.
The sample is subsequently subjected to uniform stretch in the x 1 direction, as shown in the left top corner of Fig. 3a. A uniform stress σ 0 = 1 is applied on left and right surfaces, the top and bottom surfaces are traction free. The simulated elastic fields are then used to calculate the driving force on the dislocation following the scheme outlined in subsection 3.2. A convergence study is performed to estimate the required integration area to obtain accurate results. A square integration area with a width of s is chosen, as shown in the left top corner of Fig. 3a. Figure 4 shows the driving forces calculated using Eq. (36). The result shows that the driving force is convergent when the sample width is larger than 6.
For the benchmark of the numerically calculated driving force, the Peach-Koehler force is used [5]. It is known that the driving force on a dislocation can be evaluated from the following equation where ki j is the permutation tensor,ξ j is the dislocation sense vector, and σ 0 il is the applied loading. For the uniformly stretched sample, according to Eq. (37), the driving force F i on dislocation is in the negative x 2 direction, as shown in Fig. 3a. The sample subjected to pure shear stress is also considered, and the loading is illustrated in the left top corner of Fig. 3b, where uniform shear stress τ 0 = 1 is applied on all  surfaces. For the sample under pure shear loading, according to Eq. (37), the driving force F i is in x 1 direction, as shown in Fig. 3b.
In the inset on the right bottom corner of Figs. 3a and 3b, the details of the distributed nodal configurational forces are shown. The configurational force is in principle defined as the gradient of the potential energy. Hence, Fig. 3a implies the dramatic drop of the potential energy in the compressive and tensile regions. Nevertheless, the drop on the tensile side takes over that of the compressive side. Thus the driving forces on the dislocation core is towards the tensile side. In other words, under the uniform stretch, the dislocation tends to move downwards. Similar discussion can be made on the inset of Fig. 3b, and one can conclude that the dislocation tends to move towards right under the given shear loading.
An integration area of 20×20 around the dislocation core is chosen to sum up all nodal configurational forces G i to calculate the driving forces F i . According to the convergence study discussed above, this integration area should give an accurate result to compare with the Peach-Koehler force. In Table 1 the numerically calculated driving force and the analytically evaluated Peach-Koehler force are contrasted. It can be seen that the numerical results agree well with the analytical one. The differences in both loading cases, defined as )/|F P K |, are less than 1%. We compare also the results for the dislocation core size h = 1 and h = 2. The numerical result implies that the numerically determined driving forces are not sensitive to this parameter.

Screw dislocation
In the next we consider a straight screw dislocation with the Burgers vector b = [0 0 b 0 ] with b 0 = 1, and the normal vector of the slip plane n = [0 1 0]. The sample size is 80 × 80 × 5 and the mesh size is 160 × 160 × 10. The periodic boundary condition for the displacement is applied along the x 3 direction (front and back surfaces).
The dislocation induced elastic fields are first simulated, assuming the traction free boundary condition. The stress fields are shown in Figs. 5(a-b). The displacement u 3 at the left surface is shown in Fig. 5c. Similar to the edge dislocation, good agreement is observed between the numerical solutions and the singular and non-singular analytical solutions outside the dislocation core region. The distribution of nodal configurational forces without external loads is shown in Fig. 5d. The nodal configurational forces decrease to 0 very quickly within the distance of 2h to the dislocation core.
According to Eq. (37), the Peach-Koehler force is nonzero when σ 13 and σ 23 are nonzero for the screw dislocation. Therefore, the pure shear stress is applied as τ 0 = σ i3 n i , (i = 1, 2), where n i is the normal vector of the corresponding surface. The shear stress is applied on the front and back surfaces, thus, the periodic boundary condition will not be applied for this simulation. In order to eliminate boundary effects, a thick sample along x 3 direction is used. The sample size is 50 × 50 × 50 and mesh size is 100 × 100 × 100. The straight screw dislocation line is coincident with the x 3 -axis in the center of the sample. The driving forces per unit length on the dislocation line are calculated by the total nodal configurational forces evaluated on the x 3 = 0 plane averaged along the dislocation line length. The comparison between the numerical results and the Peach-Koehler force is shown in Table 2. The difference remains below 1%. It shows that the driving forces on the screw dislocation are also not sensitive to the dislocation core size h, similar to the edge dislocation case.

Circular dislocation loop
In this subsection a circular dislocation loop is simulated. The dislocation loop is located at the x 2 = 0 plane, and the Burgers vector is b = [b 0 0 0] with b 0 = 1. The radius of the dislocation loop is 25. It exists in a cylinder with a radius and a height of 100 each. The finite element model is shown in Fig. 6. The mesh around the dislocation loop is refined with a mesh size around 0.5. Figure 6a shows the eigenstrain distribution function W (x, h) at the x 2 = 0 plane.
At the two intersection points of the dislocation loop with the plane x 3 = 0, it is of edge dislocation type, since the tangent of the dislocation line is perpendicular to the Burgers vector, whereas at the intersection points of the loop with the plane x 1 = 0 it is of screw dislocation type, since the tangent of the dislocation loop is parallel to the Burgers vector. The corresponding stress distributions confirm these features, as shown in Figs. 6(c-d), respectively. The driving forces on the dislocation loop are compared with the Cai's non-singular analytical model [22]. In the analytical model, the dislocation loop is approximated by 12 dislocation line segments connected by 12 nodes, as shown by the black line segments in Fig. 6a. The analytical driving forces on the 12 nodes are obtained by summing up the contributions from the Peach-Koehler force integrated over the segments. In the numerical model, the driving forces F i are calculated by summing up the nodal configurational forces in the dashed line box around the dislocation core shown in Fig. 6c, including the volume for θ ∈ [−π/12 + iπ/6, π/12 + iπ/6], (i = 0, 1, 2, ..., 11) in the circumferential direction.  The distribution of the 12 numerical driving forces F r along the radial direction is shown in Fig. 6a. The distributed nodal configurational forces G i are shown in Figs. 6(c-d). The results show that the driving forces F i are oriented towards the center of the dislocation loop, which means the dislocation loop tends to shrink without external loads. Fig. 6b shows the distribution of F r as a function of θ along the dislocation loop. The numerical results agree well with the analytical solutions, in both the cases of h = 1 and h = 2. The driving forces obtained with h = 1 are larger than that in the case h = 2. The reason is that the stress at dislocation core region is larger for h = 1. As result, the driving forces due to this stress field on the adjacent points will be larger. Due to the same reason, the element length along the circumferential direction can affect the accuracy of the numerically predicted driving forces. In this work, the element size on the dislocation loop is around 1.6. Fig. 6b also shows that the driving forces F r at the two intersection points of the screw dislocation type are larger than those at the two intersection points of the edge dislocation type. The Peach-Koehler forces on the loop due to the stress fields of the straight edge dislocation can be calculated as [1] where σ ns 11 is the non-singular stress field of the straight edge dislocation [22]. In Cai's non-singular continuum theory of dislocations, the non-singular stress σ ns 11 is not Cauchy stress [26], but rather an effective stress obtained through the convolution of thew(x) with the Cauchy stress (stress field of an "source" dislocation). The influence of the convolution will only be limited inside the dislocation core region, becausew(x) → 0 outside the core region. In finite element simulation, the stress field at an arbitrary point is resulting from superposition of the stress fields of "source" dislocations in linear elasticity. As a consequence, for a given dislocation core size h, the numerical stress fields will be different to σ ns i j inside the core region, as shown in the numerical results of the edge and screw dislocation above. The difference will further affect the determination of the interaction force between two very close dislocations in numerical simulations.
In the numerical simulation, the driving forces on the loop consists of two contributions. The first contribution is the driving forces due to the stress fields of the line segments of the loop, which can be calculated when the straight edge dislocation is not present, similar to the subsection 4.3. The results are shown in Fig. 7c. The second contribution is the driving forces due to the stress fields of the straight edge dislocation. This contribution can be compared with the analytical solution defined in the last equation (38). Since the driving forces on the loop due to the loop itself is always present, we simulated two problems: one with both the loop and the straight edge dislocation, and the other with only the vacancy loop. The difference between the driving forces on the loop of the two problems leads to the second contribution, which is comparable to the analytical solution (38). For numerical results, total nodal configurational forces on a point of the dislocation line are calculated by summing up all the nodal configurational forces on the cross section perpendicular to the dislocation line. Then the driving forces per unit length are calculated as the total nodal configurational forces divided by the element length in the circumferential direction.
When there is only the vacancy loop, Fig. 7c shows the nodal configurational forces are uniformly distributed along the dislocation loop in the circumferential direction. Figure 7d shows the nodal configurational forces Fig. 8 Interaction forces between the vacancy loop and straight edge dislocation. a Driving force F 2 , the attraction force (θ = 0) is larger than the repulsive force (θ = π). b Driving force F 3 , which tends to expand the loop on the straight edge dislocation increase slightly near the vacancy loop in the x 2 direction, which indicates the attraction between the vacancy loop and the straight edge dislocation. The driving force on the vacancy loop due to stress field of the straight edge dislocation are shown in Fig. 8. Figure 8a shows that the attractive forces F 2 at θ = 0 are larger than the repulsive forces at θ = π. In the meantime, Fig. 8b shows that the driving forces F 3 are positive in the radial direction. Therefore, the interaction force tends to expand the vacancy loop, as illustrated in Fig. 7a.
The comparisons between the analytical and the numerical solutions of driving forces are shown in Fig. 8. Minor differences between the analytical and the numerical solution can be found for the θ ≈ 0 and θ ≈ π, which is the bottom and the top of the loop, respectively. It may be due to the fact that for the bottom part, the vacancy loop is located very close to the dislocation core region and the stress fields are highly dependent on the dislocation core size h. At the top, the dislocation on the vacancy loop is far away from the straight edge dislocation; therefore, the stress fields of the straight edge dislocations are significantly influenced by the traction free boundary condition on a small sample. Increasing the sample size can reduce the difference between the numerical and the analytical results.

Taylor hardening
The proposed numerical scheme of determining the distributed driving force on the dislocation line is very general and can be used for complex situations. In the last example we show that the driving forces can be used to calculate the critical passing stress and determine the movement tendency of each individual dislocations in a dislocation cluster. The passing stress due to interaction between dislocations is crucial to understand the yield stress and strain hardening behavior of crystals. The Taylor hardening formula is widely used to predict these behaviors, which is defined as where τ c is the critical passing stress for passing of the oppositely signed edge dislocations in an array, as shown in Fig. 9a. b 0 = 1 is the magnitude of the Burgers vector b = [b 0 0 0]. ρ = 1/(2a 2 ) is the density of dislocations, where a = 5 is the distance between adjacent dislocations in the x 2 direction. It should be noticed that the value of a used in this numerical example is rather small and corresponds to high dislocation density. We choose the small a in order to achieve a computationally affordable but legitimate example, because increasing a implies larger simulation domain and thus higher computation cost. In fact, the constant α in the general law of the critical passing stress (39) is independent from the specific choice of the dislocation density and thus also independent from the distance a. In the case of varying slip plane distances in a complex dislocation cluster, the dislocation density is then related to the mean values of a. The analytical solution of the constant α for two oppositely signed edge dislocations (one dipole) is 0.08 for the current material. For a cluster of multiple dislocations or a dislocation cloud in general, this constant is unknown. For this complicated case, the proposed numerical scheme of calculating the driving force can be applied. In this numerical simulation, the sample size is assumed to be 50 × 50 and the mesh size is 100 × 100. The dislocation core size is h = 2. The relative position between the oppositely signed dislocation arrays is described by d. The distance between the same signed dislocations (lie in the same line in the horizontal direction) remains 2a and we vary the value of d to achieve different interaction forces between oppositely signed dislocation arrays. When the sample is under pure shear load τ 0 , the driving force on a dislocation can be calculated using the Peach-Koehler formula (37) as F P K 1 = τ 0 b 0 . In order to move this dislocation, the driving force on the dislocation due to external load should be larger than that due to the stress fields of the surrounding dislocations. Then a constant α 0 can be calculated at different d as where the driving force F 1 (d) is numerically calculated by summing up all nodal configurational forces in the integration area of 6 × 6, as shown in Fig. 9c. The constant α is determined by the maximum value of α 0 . Fig. 9b shows that the maximum α 0 is 0.175 at d = 2. The driving force on the center dislocation in the x 1 direction can be calculated as F 1 = σ 12 b 0 , where σ 12 is the stress at the core of the center dislocation.
With the Burgers vector in the x 1 direction, σ 12 is symmetric with respect to the x 2 -axis. Therefore, since one oppositely signed dislocation in a dipole contribute to α 0 = 0.08, the two oppositely signed dislocations close to the center dislocation (for a small d) already contribute to α 0 = 0.16. The numerical result α 0 = 0.175 shows that the other surrounding dislocations also contribute to larger passing stress. In reality, there are much more dislocation in a plastically deformed copper material than the simulated example, and the experimentally measured α is around 0.35 ± 0.15 [34]. For the unloaded case τ 0 = 0, the distribution of nodal configurational forces in Fig. 9c shows that the dislocation array tends to shrink in the x 2 direction while it expands along the x 1 direction. In order to show the tendency of the movement of dislocations, the maximum critical shear stress is applied on the sample, which means τ 0 = τ c | α=0. 5 . Fig. 9d shows that the external load largely changes the distribution of potential energy on the compressive and tensile site of the dislocation. The distribution of nodal configurational forces clearly shows that the oppositely signed edge dislocations tend to move towards each other. The total nodal configurational forces on the boundary should be equal to the total driving forces on the dislocation array in magnitude but in opposite direction. It shows that the total driving force on the array is along the x 1 direction, because there are more dislocations with positive Burgers vector in the dislocation array.
The results above show that the driving forces can be used to determine the critical passing stress of dislocations. More importantly, the tendency of the movement of each individual dislocation can be also clearly indicated by the nodal configurational forces. The methodology can be in principle used to study any arbitrary dislocation arrays with irregular distributions.

Conclusion
In this work, a finite element model is proposed to calculate the driving forces on dislocations. The model is based on the eigenstrain theory of dislocations, and the concept of spreading the Burgers vector around the slip plane from the non-singular continuum theory of dislocations is adopted to calculate the distribution of eigenstrain. The derivation of the Eshelby stress tensor through the gradient of the potential energy of dislocated elastic bodies is revisited based on the non-singular continuum theory of dislocations. It shows that the derivation of the non-singular Eshelby stress tensor based on the non-singular continuum dislocation theory is mathematically consistent with the derivation based on the singular dislocation theory. The model is implemented using the finite element method and applied to calculate the elastic fields and the driving forces of an edge dislocation, a screw dislocation, a circular dislocation loop, as well as the interaction force between a straight edge dislocation and a vacancy loop. For all the four cases, the elastic fields calculated numerically are compared with the corresponding analytical solutions. Moreover, the driving forces on dislocations are compared with the Peach-Koehler force. The numerical results show very good agreement with the analytical solutions. Different dislocation core sizes are considered in the simulation. It was found that although the stress fields inside the dislocation core are different for different dislocation core size, it does not significantly influence the driving forces on dislocations under external load. However, the driving forces on the dislocation loop and the interaction force between dislocations are highly dependent on the dislocation core size.
The proposed approach can be used to analyze complex dislocation configurations or clusters. It provides not only the critical information but also detailed movement tendency of each individual dislocation.