Analytical integration of the forces induced by dislocations on a surface element

An analytical formulation of the nodal forces induced by a dislocation segment on a surface element is presented. The determination of such nodal forces is a critical step when associating dislocation dynamics simulations with continuum approaches to simulate the plastic behaviour of finite domains. The nodal force calculation starts from the infinite-domain stress field of a dislocation and involves a triple integration over the dislocation ensemble and over the surface element at the domain boundary. In the case of arbitrary oriented straight segments of dislocations and a linear rectangular surface element, the solution is derived by means of a sequence of integrations by parts that present specific recurrence relations. The use of the non-singular expression for the infinite-domain stress field ensures that the traction field is finite everywhere even at the dislocation core. A specific solution is provided for virtual semi-infinite segments that can be used to enforce global mechanical equilibrium in the infinite domain. The proposed model for nodal forces is fully analytical, exact and very efficient computationally. A discussion on how to adapt the proposed methodology to more complex shape functions and surface element geometry is presented at the end of the paper.


Introduction
Over the past several decades, discrete dislocation dynamics (DD) methods have been developed into robust tools for simulating the evolution of dislocation structures in crystalline materials [1][2][3][4][5]. With the recent interest of plasticity at interfaces and in finite crystals, DD simulation methods have been applied to simulate the behaviour of finite domains. However, classical DD methods rely on Green's functions to quantify the elastic interaction between elements of the dislocation microstructure in an infinite medium. If not modified, these Green's functions, derived for infinite homogeneous media and closed dislocation loops, are not particularly well suited to plasticity problems associated in finite domains. Many analytical expressions have been derived over the years to complement the 'classical' dislocation stress expressions. Stress expressions have been derived for dislocations in half spaces [6][7][8].
Head [9] investigated the case of a straight dislocation contained in a 'bimetallic' medium with various constraints at the boundary between the two solids. Khraishi and Zbib [10] proposed a solution to ensure the traction free conditions associated with free surfaces in DD simulations. A distribution of image dislocation loops is employed to annul tractions on selected collocation points. Chu et al [11] reduced the problem of the stress field associated with an polygonal dislocation loop, contained in an anisotropic elastic semi-infinite medium, to line integrals that can be evaluated straightforwardly. Dundurs and Mura [12,13] have investigated the interaction between a straight dislocation and a circular elastic inhomogeneity with different elastic constants than the surrounding matrix. Recently, Wang and Sudak [14] considered the interaction of a screw dislocation with non-spherical inhomogeneities by means of a complex variable approach and expansions in terms of Faber series.
For more general boundary-value problems and realistic dislocation microstructures, Finite Element (FE) and Boundary Element (BE) methods have been coupled with DD [15][16][17][18][19][20][21][22] to impose desired traction-displacement boundary conditions on simulations of finite domains and to calculate the forces on the dislocations contained within the domain. Many of these hybrid methods follow the general approach proposed by van der Giessen and Needleman [15] based on the superposition principle and linear elasticity. The stress field of dislocations in a bounded domain is decomposed into the contribution of the classical stress field due to dislocation ensembles in the infinite medium, and a correction stress field that enforces the prescribed boundary conditions on the domain limits. The correction field is resolved by FE or BE methods for an elastic boundary-value problem and requires the determination of the residual tractions [23] induced by the dislocation ensemble (obtained from the classical infinite-domain stress field) on the boundary of the domain.
Two issues associated with such hybrid methods have been recently addressed in the literature [23,24]. The first issue occurs when dislocations terminate at the boundary of the domain, the stress field calculated from the classical line integral over such open dislocation loop may not be divergence free [23,24], and the integration of the residual tractions on the entire domain boundary may not cancel out [23]. In other words, mechanical equilibrium is not satisfied. The divergence of the infinite-domain stress field can be removed by closing the dislocation network with the help of virtual segments added outside from the bounded domain [21,23,24], or through the use of rational stress expressions [25]. As the exact virtual segment arrangement has little influence on the final solution of the boundary-value problem, some authors [21,23] prefer to close the dislocation network by simply extending segments, that end at the boundary, to infinity. The second issue has to do with the fact that the classical infinite-domain stress field is singular at the dislocation core. This poses some numerical issues especially when evaluating the residual tractions due to segments in contact with the domain boundary. To circumvent this difficulty, a number of authors [20,21,23] now use the non-singular expressions proposed by Cai et al [26] to express the infinite-domain stress field.
The evaluation of the residual tractions on the domain boundary, which is a quantity that is passed on to FE or BE simulations, constitutes a critical step of such mixed DD-continuum approaches. This procedure, which requires the integration of the continuous infinite-domain stress field due to the dislocation ensemble on surface elements and the proper transformation into a force at FE or BE nodes, has been overlooked until now. To date, analytical expressions have not been developed for the forces at FE or BE nodes due to a dislocation line (hereafter referred simply to as nodal forces). Therefore the nodal forces on surface elements are evaluated by means of a numerical integration on the surface of the element and the quality of this integration will dramatically affect the precision of the FE or BE solve [24,27].
The numerical cost associated with the evaluation of nodal forces on FE or BE is another aspect to consider. For a start, the cost of a numerical integration grows linearly with the number of quadrature points (QPs) required to achieve a given accuracy. Then, this numerical procedure is performed for every dislocation segments constituting the microstructure and every surface elements delimiting the domain. The evaluation of nodal forces on FE or BE is thus often considered to be one of the most expensive procedures when combining DD and continuum approaches. In the case of the interaction of two dislocation segments, the separation distance between the two dislocations determines the spatial variation of the stress field due to one dislocation onto the other. As a consequence, the number of QPs required to achieve a given accuracy has to be dramatically increased when the two dislocations get closer [4]. A similar behaviour can be expected when the dislocation ensemble approaches a surface. It seems difficult to control integration errors and what the appropriate number of QPs should be.
In this study, we show that the nodal forces induced by a dislocation segment of arbitrary orientation on a surface element can be analytically expressed for linear elastic isotropic solids. The model involves a sum of triple integrals (along the dislocation line and the surface element) that are solved by means of a sequence of integrations by parts given as recurrence relations. The proposed model exhibits the following advantages.
• The use of the non-singular continuum theory of dislocation [26] yields to stress and traction fields that are defined and finite everywhere in the infinite domain, even at the dislocation cores.
• The model is fully analytical and exact. The proposed formalism based on recurrence relations makes it very easy to implement. It is also very efficient computationally.
• A specific solution has been derived for semi-infinite segments. The determination of the asymptotic behaviour of the nodal forces at infinity is done analytically. This type of segments can be used to extend open dislocation loops in the infinite domain.
• The adopted formalism and geometry are chosen to match most DD and FEM implementations.
The second objective of this paper is to investigate the integration errors associated in the numerical evaluation of the nodal forces. This study will be performed on a simplified surface element geometry, and a discussion as to how the technique can be adapted to more complex surface elements required by arbitrarily shaped finite domains will be offered at the end of the paper.

Problem formulation
Consider a linear elastic body V containing a dislocation ensemble as illustrated in figure 1.
The body is subject to traction boundary conditions T = T 0 on surface S T and displacement boundary conditions: u = u 0 on surface S u . Following the formulation proposed in [15] to impose traction-displacement boundary conditions on a DD simulation within a finite domain, the total displacement and stress fields are written as The (∞) fields are elastic fields associated with the dislocation network in an infinite medium. The (∞) fields are themselves a superposition of the fields of each individual dislocation. These fields are governed by the standard equations of linear elasticity. The (ˆ) fields or 'image fields' are superimposed to enforce boundary conditions. To put it differently, a finite domain simulation is equivalent to having an infinite domain with an internal surface where a prescribed traction or displacement condition is applied. The boundary conditions at the boundary S read [15] where n is the outer unit normal vector to S. This method, based on the superposition principle [15], is therefore limited to small elastic strains. A discussion on the advantages and drawbacks of this method can be found in [28]. Solving a boundary-value problem in the example of free surface boundary conditions requires the following steps [15,23].
• In these simulations, one must ensure that the dislocation microstructure provided by DD is valid in the infinite system (associated to a divergence-free stress field) as it will be used as a starting point for the finite domain simulation. This requires the introduction of virtual segments: segments that exist in the infinite domain but not in the finite domain. Semi-infinite virtual segments can be used to close loops between segments that intersect the surface. The association of the dislocation microstructure contained within the domain and the virtual segments in the infinite domain is sometimes called 'augmented dislocation configuration' [23]. • The unbalanced tractions due to the dislocation ensemble are evaluated at the domain boundary. This calculation concerns dislocations contained within the domain and virtual segments in the infinite domain. The sign of the resulting traction field is reverted and is appropriately transformed into a set of nodal forces on nodes that are at the boundary of the domain. • Then, the classical system of elastic equations is solved by FE or BE methods for the image displacement field, from which the strain tensor is determined and Hooke's law gives the image stress field. • The image stress field is finally applied to the dislocation network and DD will simulate the dislocation microstructure evolution. Calculations in step 2 of the surface element nodal forces induced by the dislocation microstructure are not straightforward. The stress field σ ∞ (x) due to the dislocation ensemble is defined and finite everywhere in the infinite domain [26]. The resulting traction field σ ∞ (x) · n(x) is defined and finite as well at the domain boundary. Therefore, the traction field due to the dislocation ensemble can be seen as a distributed load on the surface (in that case loads have units of pressure). To match the discrete formalism of FE or BE methods, the traction field, or distributed load, must be appropriately transformed into a set of nodal forces [29]. This is usually accomplished such that the mechanical work done by the distributed load over the element is equivalent to the mechanical work performed by the nodal loads. The equivalent nodal force on node n is thus written [29] where dS is an infinitesimal surface element and S (E) is the surface area of the surface element. N (n) (x) are the so-called shape functions, which assure partition of the traction field among the nodes of the surface element.
To avoid the core singularity associated with the classical Volterra dislocation, the nonsingular formulation proposed by Cai et al [26] is used to express the stress field σ ∞ (x) of the dislocation ensemble in the infinite domain. In an homogeneous and infinite linear elastic medium, the stress field associated to a dislocation loop can be expressed as a contour integral along the loop [13]: where C ij kl is the elastic stiffness matrix, lnh the permutation operator, b the Burgers vector and x is the coordinate that spans the dislocation. G kp (x − x ) is the Green's functions of elasticity [13]. G kp (x − x ) is defined as the displacement component in the x k direction at point x due to a body force applied along the x p direction at point x . The singularity of the conventional stress field of dislocations is due to the fact that the Burgers vector distribution is taken as a delta function. Cai et al replaced the sharp delta distribution by a wider-spread distribution. As a realist Burgers vector distribution would be material dependent and has yet to be determined, i.e. from atomistic calculations, Cai et al proposed a simple but practical isotropic spread mainly contained in a radius a about every points of the dislocation 5 . In the case of isotropic elasticity, the Green's function now takes the following form: 5 The Burgers vector distribution is actually not bounded in space, but it decreases as fast as 1/R 7 and therefore becomes quickly negligible far from the dislocation. This will be discussed in detail in a forthcoming paper.
x 1 x 2 x x p q x 3 x 4 x 5 Figure 2. Definition of the considered geometry and associated variables. The surface element, of rectangular shape, is delimited by four nodes of Cartesian coordinates: x 3 , x 4 , x 5 , x 6 . p and q are unit vectors used to define any points x on the surface element. The dislocation segment is described by a line vector t and Burgers vector b, and is spanned by coordinate x .
where R = ||x − x ||, and µ and ν are the isotropic shear modulus and the Poisson's ratio, respectively. R a is the convolution of R with the Burgers vector distribution. R a is simply R 2 k + a 2 in the case of the isotropic distribution proposed in [26].
In DD, the dislocation microstructure is commonly discretized into straight segments. The contour integral of equation (6) becomes for a straight segment of non-singular dislocation (omitting the ∞ superscript for now) [26]: σ (12) The superscript 12 refers to the end nodes delimiting the dislocation segment. x 1 , x 2 are Cartesian coordinates delimiting the dislocation segment. This equation is presented in a form independent of any coordinate system ('dyadic' form). The analytical solutions of these integrals are given in [26] and are recalled in appendix C. In principle, a can be defined from atomistic calculations to match the core size of the dislocation structures. In practice, a is typically small (and positive), of the order of one b magnitude. When a equals zero, the well-known singular stress expressions are recovered.
In this study, we show that the nodal force F (n) in equation (5) that is required to solve the boundary-value problem can be analytically evaluated using the non-singular stress expression σ ∞ (x) of equation (8) provided by Cai et al [26]. However, as this study represents a first attempt of this kind, we first consider a simple surface element geometry and simple shape functions. The surface element has a rectangular shape and is delimited by its Cartesian coordinates x 3 , x 4 , x 5 and x 6 as shown in figure 2. x is the coordinate that spans the surface of the element. Recalling that x is the coordinate that spans the dislocation segment delimited by x 1 , x 2 . Vector R = ||x − x || can also be written as R = yt + rp + sq, with t the dislocation line vector and p = x 4 −x 3 ||x 4 −x 3 || and q = x 5 −x 3 ||x 5 −x 3 || (see figure 2 and appendix A). For a rectangular surface element, p · q is obviously null. The decomposition of R = yt + rp + sq is valid for dislocation segments that are non-parallel to the surface element: t · (p × q) must be different from zero. A different decomposition of R could be formulated for the case when the dislocation is parallel to the surface element. However, this represents a lot of effort to derive one specific solution. A more practical solution is proposed instead. The nodal force for a segment parallel to the surface can be approximated from the present work by averaging the nodal forces obtained for a set of situations where the dislocation is slightly tilted from the parallel case.
Functions N (n) (x) are chosen as linear shape functions. Taking node 6 as an example, N (6) is defined as where r 1 , r 2 and s 1 , s 2 are the lower and upper limits of the r and s scalars, respectively. Inserting equation (9) into equation (5), the nodal force at node 6 becomes Finally, replacing the stress field σ (12) by its non-singular expression of equation (8) and recalling that R = yt + rp + sq leads to It becomes clear that finding an analytical solution for the nodal force means solving all the triple integrals that are introduced in equation (11). This is the subject of the next section.
3. Analytical expressions for the nodal force due to a dislocation segment on a surface element

The finite segment
In this section, we first present the analytical integration of the stress field induced by a segment of finite length. The segment can be part of closed loops that are entirely contained inside the body or in contact with the surface element for opened dislocation geometries. The configuration for the dislocation segment and the surface element as depicted in figure 2 is reused here. The nodal force expression as derived in equation (11) is a sum of 44 triple integrals of the form dr ds dy (12) with l integer = 3 or 5 and i, j, k integers ∈ [0, 3].
We show next that these triple integrals can be obtained by a sequence of integrations by parts that present specific recurrence relations. First, we recall that R a = R 2 k + a 2 and R 2 k = R · R = r 2 + s 2 + y 2 + 2cry + 2 dsy where c = p · t and d = q · t. Then, the following partial derivatives of r i s j y k R l a are introduced and will be the starting point of the future integrations by parts: ∂ ∂s ∂ ∂y Three successive integrations over r, s, y of these partial derivatives lead to (omitting integral bounds for clarity) a dr ds dy integrals appear in the right-hand side of equations (16), (17) and (18). Therefore, these equations also constitute a system of three linear equations with three unknowns: the integrals of type r i s j y k R l+2 a dr ds dy. Solving the system of equations (16)-(18) by substitution leads to Here, we chose to always solve the system of equations (16)- (18) for the unknown r i s j y k+1 R l+2 a dr ds dy first. In addition to these three equations, a fourth recurrence relations can be obtained from the definition of R 2 a = r 2 + s 2 + y 2 + 2cry + 2dsy + a 2 and is as follows: The recurrences relations (19) through (22)  dr ds that are seen in the right-hand side of these equations. It can be noted that monomials are present in front of the double integrals. These double integrals need to be found as well and this can be achieved by a procedure identical to the one proposed for triple integrals. That is to say, the partial derivatives of equations (13)- (15) are integrated twice over r and s or r and y or s and y to obtain recurrence relations for double integrals. This second sequence of integrations by parts introduces, this time, single integrals. The single integrals are also obtained by recurrence relations derived from one single integration of the partial derivatives of equations (13)- (15). The development is detailed in appendix A.
After the complete sequence of integrations by parts, only a few integrals cannot be expressed by any of the recurrence relations proposed and need to be formally expressed. They correspond to three single integrals, three double integrals and one triple integral without any polynomial in their numerator: This set of integrals constitutes seed functions, from which all other integrals are constructed. Analytical solutions were successively obtained for all the linear and the double integrals but not for the triple integral of the form 1 R 3 a dr ds dy. Any first two integrations of 1 R 3 a leads to a convoluted arctangent function. Since this function is continuous and well behaved (if a > 0), its antiderivative must exist. However, we were not able to find an analytical expression for it, even with the help of the main softwares for symbolic calculation. Because of the recurrence relations, this integral appears in many other triple integrals. Fortunately, the terms that are function of 1 R 3 a dr ds dy all cancel out perfectly in the nodal force expressions.
Since the nodal force is a vector, proving that F (n) is independent of 1 R 3 a dr ds dy is equivalent to show that every components of this vector are independent of 1 R 3 a dr ds dy. To facilitate the demonstration, the basis made of unit vectors p, q and n = p × q is used to express the components of the nodal force components. By developing the triple integrals recurrence expressions in the definition of the nodal force given by equation (11), one can show that all terms that are function of 1 R 3 a dr ds dy exactly cancel out for the three nodal force components. The details of this demonstration are given in appendix B. As a consequence, there is no need to numerically integrate 1 R 3 a dr ds dy, and its contribution is set to 0. The present model for nodal forces is therefore fully analytical and exact.
In practice, the various steps presented above are done in a different order. The six seed functions of equations (23) are calculated first, and then are used to build the single integrals. The single integrals are themselves used to built the double integrals and the double integrals are used to build the triple integrals that are the objects actually required in the nodal force expression of equation (11).

The semi-infinite segment
This semi-infinite segment can be used to extend open dislocation loops in the infinite domain. The configuration under consideration is now as follows. Node x 1 of the dislocation segment is put in contact with the surface of the element and the other end node x 2 is taken to +∞, as depicted in figure 3. x 4 x 5 Figure 3. Definition of the considered geometry and associated variables. The surface element, of rectangular shape, is delimited by four nodes of Cartesian coordinates: x 3 , x 4 , x 5 , x 6 . p and q are unit vectors used to define any points x on the surface element. The dislocation segment is described by a line vector t and Burgers vector b, and is spanned by coordinate x .
Dealing with a semi-infinite dislocation is equivalent to having one of the integral bounds, here y 2 , set to +∞ in the stress expression of equation (8) or the nodal force on a surface element of in equation (11). However, attention must be paid when dealing with such infinite integrals, and one has to verify that the integration result will be finite. So, we first express analytically the stress field due to a semi-infinite dislocation starting with the non-singular equation (8) and taking x 2 (and y 2 ) to +∞: Limits of the antiderivatives solution of the integrals shown in the previous equation are obtained and detailed in appendix C. The resulting stress field σ (1+∞) for a semi-infinite dislocation is defined and finite everywhere in the infinite domain, the associated traction field σ (1+∞) · n is defined and finite for any surface element at the boundary of the body. Consequently, as the integrand ([σ (1+∞) · n] × N (n) ) in the nodal force expression of equation (11) is continuous and finite, the nodal force expression as given in equation (11) is still valid and constitutes the starting point for the derivation of expressions specific to the semi-infinite segment. We recall that equation (11) is composed of a sum of 44 triple integrals of the form r i s j y k R l a dr ds dy. Let us call F (r, s, y) the antiderivative of f (r, s, y) = r i s j y k R l a along y. In the case of a semi-infinite segment, they can be written The second term in the right-hand side of equation (25), corresponds to the 'finite end node' of the dislocation segment. Therefore, use can be made of the same integrations by parts and analytical functions as in the case of the finite segment. Now considering the first term of the right-hand side of equation (25) and recalling that the limit of a (converging) sum is the sum of the limits, we simply need to find the limits of the few analytical seed function in equations (23), while preserving the recurrence relations to build the semi-infinite nodal traction. This does not introduce any new integration constant of the type C(r, s) as the integration path for y 1 and y 2 remain essentially identical. Therefore, only limits lim y→+∞ for the set of six seed functions that are analytically solved in equation (23) need to be formally expressed. Attention has been paid to the determination of the asymptotic behaviour at infinity of the seed functions of equation (23). Because these functions are convoluted functions of y, they are systematically approximated by a Taylor series of the order of n associated with an error of the form O( 1 y m ) with m a positive integer. The order of the Taylor series is chosen so that the error exponent m is always bigger than 1. To illustrate this point, we consider the seed function

Comparison with a numerical evaluation of nodal forces
In this section, the proposed analytical expressions for the nodal force induced by a dislocation segment on a surface element is compared with a conventional numerical integration for two purposes. First, our analytical work constitutes a clear reference to evaluate how fast a numerical integration converges as the number of QPs is increased. Second, the numerical cost of the two methods will be quantified in order to evaluate when it is preferable to use one technique over the other.
The numerical quadrature is performed with the help of two Gauss-Legendre schemes [30,31], the nodal force at node n is approximated as follows: where r i and s j are abscissa locations where the function [σ (12) is to be evaluated and w i and w j are weighting functions. M is the number of function evaluations, the total number of QPs is M 2 . Relations to determine abscissas r i , s j and weights w i , w j from the desired number of function evaluations M can be found in [30]. When N (n) are linear shape functions, equation (26) becomes, i.e. for node 6, In the last equation, σ (12) corresponds to the stress field induced by a finite or semi-infinite dislocation segment. The double Gauss-Legendre integration (simply referred to as 'Gaussian  quadrature' later in the text) is thus carried out only on the surface element and not along the dislocation line.
The configuration considered to achieve the comparison is as follows. The dislocation segment has a size of 400 b (||b||) which is a typical discretization length for a density a 10 12 m −2 . The surface element is taken as 4000 × 4000 b 2 . Interestingly and for a given number of QP, only two parameters seem to affect in practice the numerical error in the Gaussian quadrature: the distance D between the centre of the dislocation and the centre of the surface, and the dislocation character. Indeed, the numerical error is smaller for screw dislocation as it is for edge, probably because the stress field of screws is simpler. Therefore, for the purpose of the comparison, we will be varying the dislocation/surface distance, and the dislocation chosen to be of edge character.
The relative error associated with the Gaussian quadrature is shown in figure 4 as function of the number of QP for various dislocation-surface spacings. This relative error is calculated in taking the result of our analytical nodal force expressions as reference values. Every curve in figure 4 exhibits the same behaviour, with a 'slow' then 'fast' decrease of the relative error as the number of QP used in the Gauss quadrature increases. Then the numerical error reaches a plateau at about 10 −12 , which is close to the precision limit allowed by the double precision floating point that is used for our numeric values (12 digits after decimal point). The rate at which the relative error decreases with the number of QP is strongly dependent upon the dislocation-surface separation distance. For instance, the relative error in the Gauss quadrature is still about 1% with 10 6 QP when the dislocation is in contact with surface, whereas only 49 QP are sufficient to reach the maximum numerical precision of 10 −12 when the separation distance is now of 10 000 b.
A similar evaluation of the numerical error associated with the Gaussian quadrature is now performed in the case of a semi-infinite segment for which the finite end node x 1 is placed on one of the surface node, i.e. x 6 , the other end node of the segment x 2 is sent to +∞. It can be noted that this configuration would constitute the most challenging geometry to properly described with the classical singular theory of dislocation. The line direction t coincides with the surface normal, and the dislocation exhibits again an edge character. The integrand used in the Gauss quadrature [σ (1+∞) (r i , s j )·n] N (n) (r i , s j ) is now a function of the stress field for Figure 5. Relative error between a numerical integration of the nodal forces associated with a piercing semi-infinite dislocation segment and the analytical integration proposed in this work. The dislocation, in contact with the surface in x 6 , is normal to the surface and has an edge character. The surface element has a square shape and a size of 400×400 b 2 . The error is given as a function of the number of QPs used. a semi-infinite segment σ (1+∞) given in appendix C. As before, the double Gauss-Legendre integration is carried out only over the surface element and not along the dislocation line.
The relative error in the nodal force calculation obtained by a Gauss quadrature for the semi-infinite segment is displayed in figure 5. The curve exhibits a behaviour similar to the one obtained for finite segments. As expected, the relative error decreases as the number of QP is increased. The negative slope is at first relatively small and increases with the number of QP until the error plateaus at a value close to 10 −12 at about 4 × 10 4 QP. Surprisingly, the Gaussian quadrature for the piercing semi-infinite segment converges towards our analytical solution much faster than the piercing finite segment shown in figure 4 (curve in red).
To compare the computational cost of these two approaches on equal footing, an independent cost reference is introduced. It is chosen to be the cost associated to one single stress at a point σ (12) evaluation using equation (8). For this calculation, the segment length (finite or semi-infinite) and position x 1 , x 2 will be the same as those used for nodal force calculations; the point at which the stress is evaluated is now fixed in the centre of the surface element.
Figures 6 present the cost comparison for, respectively, a finite and semi-infinite segment. Our analytical model is 14 times more computationally expensive than a single stress at point calculation for a finite segment, and only 4 times more expensive in the case of a semiinfinite dislocation. As expected, the computational cost associated with the Gauss quadrature increases considerably with the number of QP independently of the distance from the surface; it rapidly becomes several orders of magnitude more computationally expensive than a single stress at point calculation. Finally, the domain where a Gauss integration is less expensive than the proposed analytical model is very narrow and corresponds to configurations for which less than 9 QP are sufficient to evaluate the traction.

Discussion and concluding remarks
In this study, we propose an analytical expression for the nodal forces associated with the tractions generated by straight dislocation segments on a rectangular surface element. The determination of these nodal forces is a critical step in hybrid methods associating DD and FE or BE techniques to solve boundary-value problems. The starting point of the nodal force derivation is the non-singular stress at a point formulas proposed in [26] that present the advantage to be defined and finite everywhere in the infinite domain even at the core of dislocations. In the case of linear isotropic elasticity and simple linear shape functions, nodal forces are derived as sums of triple integrals of the form r i s j y k R l a dr ds dy that need to be obtained. Partial derivatives of r i s j y k R l a dr ds dy with respect to r, s and y are the starting point of a sequence of integrations by parts to solve the required triple integrals. The proposed model for nodal forces is analytical solely, which makes it exact and very computationally efficient. The solution of the required triple integrals, obtained by recurrence relations, is also straightforward to implement in a program.
Despite a formulation that can seem complicated, based on a hundred or so integrals, the analytical surface integration proposed in this work is very advantageous from a computational point of view. In addition to this benefit in term of performance, it is clear from figures that the analytical integration of stress fields is the only solution capable of providing a precise traction measure when dislocations get close to the surface element. When considering a realistic microstructure, the closest segments can represent the biggest contribution to the total traction on a surface element. It is clear from figures 4 and 5 that the error in a conventional numerical integration would be difficult to control as it depends on the dislocation-surface distance and on the character. Adapting the number of QP targeting a given precision can be easy to achieve for simple dislocation configuration, but is probably more difficult when the dislocation microstructure evolve and discretization becomes finer. For segments very far from the surface element, a numerical integration could be used as an alternative, as the number of necessary QPs becomes minimum. However, for segments that are that far away, it might be beneficial not to perform an N 2 calculation at all but use instead Taylor expansions of the stress following a fast multipole method framework [32,33].
Since this work is our first attempt to solve the surface integral problem, the shape functions and the surface geometry considered remain simple and below general FEM standards. We discuss next how to adapt the present analytical work to consider quadratic shape functions and various surface elements. First, let us consider a eight-node rectangular surface elements. This type of elements contains four mid-side nodes x 7 , x 8 , x 9 , x 10 in addition to the four corner nodes already present in a simple linear rectangular element. The nodal forces N (n) become for a corner node, i.e. x 6 [29], N (6) = 1 4 (1 + r)(1 + s)(r + s − 1), and for a mid-side node, i.e. x 9 (centre, top), In these equations, r and s ∈ [−1, 1] and a proper rescaling would be required for a general configuration. Replacing the linear shape functions by these quadratic functions in the nodal force expression of equation (11) simply leads to higher polynomials of r and s. In other words, the determination of the nodal force requires that we solve triple integrals of the form r i s j y k R l a dr ds dy (28) l and k remain unchanged, with l = 3 or 5 and k ∈ [0, 3] Few new integrals are required and are associated with higher exponents i and j for r and s, but the recurrence relations that were used to obtained all the other r i s j y k R l a dr ds dy integrals are still valid and can be used for these additional triple integrals required for the high-node rectangular element. Also, because the recurrence equations express integrals with polynomials of a given order r i s j y k as functions of integrals with polynomials of low order, solving these new triple integrals will not introduce any new seed functions.
Going from a rectangular surface element towards a parallelogram element involves a bit more modification. Since the definition of p = x 4 −x 3 ||x 4 −x 3 || and q = x 5 −x 3 ||x 5 −x 3 || are preserved, p · q now differs from 0. R a becomes (y 2 + r 2 + s 2 + 2cry + 2 dsy + 2(p · q)rs + a 2 ) 1/2 , where an additional cross-term 2(p · q)rs can be seen. Additional terms are also obtained in the partial derivatives of equations (15) . Assuming that all other integrals are known from previous integrations by parts, these three equations constitute again a linear system of equations with three unknowns that can be easily solved by substitution. Nodal forces on a parallelogram surface element can be found following the same treatment as the rectangular element, but every recurrence relations must be reworked to account for the additional crossterms. This procedure should not introduce any new seed function, however, because the definition of R a has changed (but not R), the analytical solutions obtained for the seed function need to be reworked as well.
Integration of the traction field over a general triangular surface element introduces other modifications to the formulation proposed in this work. First, taking a linear triangular element as an example, the shape functions N (n) slightly differ from the linear shape functions used for the linear rectangular element [29]: N (n) = A + Br + Cs, where A, B and C are now constants that depends on the geometry of the triangle. Second, since p and q are taken as unit vectors aligned with two sides of the triangular element, p · q = 0 like in the case of the parallelogram. R a includes again an additional cross term 2(p · q)rs. Third, the integral bounds r 1 , r 2 and s 1 , s 2 are not constant any more and are linear functions of s and r, respectively: r 1 = f (s), r 2 = g(s), s 1 = h(r) and s 2 = j (r) (in the most general case). The order of the sequence of integration over r and s now becomes important. Consequently, two different analytical solutions are now required for every seed function depending on the lower and upper bounds of the integral (over r or s) that is performed first. All other integrals are still build by integration by parts and starting from the seed functions, although the order of integration needs to be memorized. It must be noted that any quadrilateral element can be built from the combination of two triangular elements.
The feasibility of applying similar methodology to obtain nodal forces on a non-planar surface elements is finally discussed. In the case of non-planar surface element, the surface normal n now varies spatially. One of the conditions to solve analytically the nodal forces was to provide a rather simple expression for vector R. For some specific types of non-planar surface element, it may still be possible to decompose or approximate R by a sum of elements, and recurrence relations may even be obtained. It is expected that the number of recurrence relations and seed functions to obtain will greatly increase with the complexity of the geometry of the surface elements under investigation. Eventually, some of the seed integrals may require an approximation by a numerical integration procedure. The cost associated with the evaluation of seed functions is a small cost to pay in comparison to the gain associated with the analytical work done thanks to the recurrence relations. Therefore, the proposed methodology has a good chance to be still advantageous computationally and in terms of precision when compared to a numerical integration of the nodal forces over the surface element.

Acknowledgments
This work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Appendix A. Nodal forces induced by a finite segment on a rectangular linear element of surface
In this section, we provide analytical expressions for the nodal forces due to a finite segment on a planar surface element. The configuration under investigation is depicted in figure 2. The dislocation segment has a Burgers vector b and a line vector t and is associated with an infinite domain stress field σ ∞ (x). The variable x spans the segment between x 1 and x 2 . We consider a surface element that is one of the element describing the boundary of a finite domain. The rectangular surface element has a constant normal n over its surface and is delimited by nodes x 3 , x 4 , x 5 and x 6 . Variable x spans the surface element.
The continuous traction field σ ∞ (x) · n, due to the dislocation segment, is transformed into a set of nodal forces F (n) on node n of the surface element as follows: where S (E) is the surface of the element and dS is an infinitesimal surface element. N (n) (x) are the so-called shape functions associated with node n. In the case of a linear surface element, N (n) is defined as (i.e. of node 6) where r and s are scalars that span the surface of the element. r 1 , r 2 and s 1 , s 2 are lower and upper bounds of r and s, respectively. The infinite-domain stress field σ ∞ (x) is expressed following the non-singular framework proposed by Cai et al [26]. In the case of linear elastic isotropy, the stress field σ ∞ (x) is as follows: the vector R can be written in a coordinate-independent manner: Inserting the expression of the stress field of equation (A.3) and the shape function of equation (A.2) into the nodal force expression leads to The nodal force becomes when using the definition of R of equation (A.6): with A = (r 2 −r 1 )(s 2 −s 1 ). We can further rearrange the nodal force expression by introducing the following vectors and integrals: The nodal force is thus a function of 44 triple integrals of the form We show next that these triple integrals can be obtained by a sequence of integrations by parts that present specific recurrence relations. To this aim, let us introduce the following partial derivatives of r i s j y k R l a that constitute the starting point of the future integrations by parts: dr ds dy are shown in these equations. Consequently, these equations also constitute a system of three linear equations with three unknowns. The solution of this system is easily obtained by substitution: where we are now using the short notation for the integrals and D ij l , E ikl and F jkl are double integrals defined latter. A fourth recurrence relations is obtained from the definition of R 2 a = r 2 + s 2 + y 2 + 2cry + 2dsy + a 2 : The double integrals D ij l , E ikl and F jkl are defined as Assuming that the single integrals on the left-hand side and all the other double integrals are known from previous integration by parts, these two equations constitute a set of recurrence relations to express r i+1 s j R l+2 a dr ds and r i s j +1 R l+2 a dr ds as functions of integrals of rank 'l' and rank 'l + 2' with low i and j . These two equations are also a system of two equations with two unknowns that can be solve as follows: where A il and B jl (and C kl required latter) are single integrals. An additional recurrence relation is again obtained from the definition of R 2 a = r 2 + s 2 + y 2 + 2cry + 2dsy + a 2 : Single integrals A il , B jl and C kl are defined as In that case, exponents i, j , k and l now take the following values:

A.3. Recurrence relations for single integrals
Finally, one integration of equations (A.17)-(A. 19) gives directly a set of recurrence relations for linear integrals: Once more, the definition of R a gives a last set of recurrence relations: These integrals constitute a set of seed functions from which all other integrals are built. These integrals are simple enough to be analytically solved as follows: We were unable to solve analytically the last seed integral of the form 1 R 3 a dr ds dy. Indeed, any first two integrations lead to a convoluted arctangent or hyperbolic arctangent function. For example, in the case of two successive integrations over r and s: These arctangent functions are continuous, finite and well-behaved. Consequently, by virtue of the first fundamental theorem of calculus, an antiderivative should exist. However, we were unable to find any analytical antiderivatives for these arctangent functions. A numerical integration procedure could easily provide an approximate solution for this integral, but this is actually not required in practice as all terms that are dependent of 1 R 3 a dr ds dy exactly cancel out in the final nodal force expression (see appendix B). Consequently, H 000−3 is simply set to zero.

Appendix B. Proof that the present model is fully analytical
In appendix A, we showed that the nodal forces due to a finite segment on a rectangular surface element, that involve triple integrals (along the surface and the dislocation), can be analytically solved using a set of recurrence relations and seven seed integrals. One of the seed integrals of the type H 000−3 = 1 R 3 a dr ds dy cannot be analytically solved. We show next that all the terms that are function of this integral actually cancel out in the final nodal force summation.
Let us recall that the nodal force on node n due to a finite segment on a rectangular surface element as depicted in figure 2 can be expressed in a compact form as where I ijkl are constant vectors that are functions of the dislocation Burgers vector b and line vector t and the surface element geometry: p, q and normal n. J ij kl are defined as sums of r i s j y k R 3 a dr ds dy triple integrals including H 000m3 . Because the nodal force is a vector, proving that F (n) is independent of H 000m3 is equivalent to show that the three components of the vector are independent of H 000m3 . The choice concerning the basis of vectors used to define the nodal force components is open. The use of the unit vectors p, q and n = p × q as a basis will reduce the length of the demonstration. To start, we consider the component F (n) ·p of the nodal force along the p axis.
First, let us show the list of triple integral where L ij kl gather all other terms appearing in J ij kl integrals that are not function of H 000−3 ; and r i = r 1 or = r 2 and s j = s 1 or = s 2 depending on the surface node under consideration. Second, the component of the nodal force F (n) ·p along the p axis is written as Replacing vectors I ij kl by their expressions leads to It yields after some simplifications and rearrangements, It can be easily seen from the definitions of J ij kl that the H 000−3 terms cancel out two by two in the first six lines. F (n) ·p reduces to We now focus onto the last three lines that are still function of H 000−3 . The following dot products of b are introduced: Replacing the mixed products [(b × t) · p], [(b × t) · q] and [(b × t) · n] by functions of f , g and h into F (n) ·p leads to where M F·p contain all terms of F (n) ·p that are not functions of H 000−3 . The last J ij kl integrals are now developed: Recalling that e 2 = 1 − c 2 − d 2 , the final sum for the s i terms is The component along p of the nodal force is therefore independent of H 000−3 . This independence can be obtained in a similar way for the two other components F · q and F · n of the nodal force. Consequently, the nodal force is independent of the integral H 000−3 .

Appendix C. Stress at a point induced by a semi-infinite straight segment
In this section, non-singular stress expressions for a semi-infinite dislocation in an isotropic medium are given. The dislocation segment is described by a Burgers vector b and a line vector t and is delimited by nodes x 1 and x 2 . Following the formulations given in [26], the stress at a point x associated to a finite segment is σ (12) Replacing R by its definition into the stress field expression leads to σ (12)  where (C.23)

Appendix D. Limits of the seed functions for the semi-infinite segment
In this section, we provide limits of the seed functions that are required for the nodal force due to a semi-infinite segment. The configuration under investigation is depicted in figure 3, where x 1 remain a finite end note and x 2 is sent to infinity. All the recurrence relations derivated for the finite segment can be reused for the semi-infinite segment, only the set of seed functions will be different for y 1 or y 2 .
As the seed functions are convoluted function of y, the asymptotic behaviour at infinity of these functions is determined by means of Taylor series expansions. For a given seed function, Taylor series expansions are employed until the functions of y that are obtained cannot be decomposed any further. The order of the Taylor series expansions is chosen with care so that the error follows O( 1 y ). As can be seen from the recurrence relations given in appendix A, the seed functions are often multiplied by a monomial of r, s or y to build other integrals. If these integrals are single or double integrals, they can themselves be used to build other integrals and other monomials can be introduced. The complete list of monomials for every seed functions is as follows: Monomials of y γ are obtained in the case of A 0−1 , B 0−1 and D 00−3 . Consequently, the order of the Taylor series expansions for these seed functions must be increased when they are used with monomials of y γ to maintain an error of the type O 1 y .

D.1. Taylor expansions of the seed functions
We introduce the following variables: R · p = r + cy R · q = s + dy R · t = y + cr + ds.
For practical reasons, we focus on configurations associated with y → −∞. The limits of the seed functions in the case when y → +∞ can be derived by symmetry from the case when y → −∞. We pose Obviously, x goes to 0 − when y goes −∞. The following classical Taylor series expansions will be used: By using up to three of the previous Taylor series expansions, the seed functions are approximated as follows:  with C = f s − gr (D.14)
Some of these limits are functions of polynomials of y or logarithm ln(|y|) that go to infinity when y 2 goes to infinity. Many of the terms that go to infinity are only functions of y, y and r or y and s. These terms will cancel out when the seed integrals are calculated using the integral bounds r 1 , r 2 and s 1 s 2 . Few terms in these limits contain monomials of the three variables r, s and y. However, for reasons of symmetry, they exactly cancel out in the final nodal force evaluation. Consequently, the nodal force of the semi-infinite segment is finite and every terms that are functions of y 2 could be removed or set to any finite values in the limits expressions.