{'fraction_non_alphanumeric': 0.034700163634613534, 'fraction_numerical': 0.015641543940706517, 'mean_word_length': 4.738470035901685, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing forces on interface elements exerted by dislocations in an elastically anisotropic crystalline material

Driven by the growing interest in numerical simulations of dislocation-interface interactions in general crystalline materials with elastic anisotropy, we develop algorithms for the integration of interface tractions needed to couple dislocation dynamics with a finite element or boundary element solver. The dislocation stress fields in elastically anisotropic media are made analytically accessible through the spherical harmonics expansion of the derivative of Green's function, and analytical expressions for the forces on interface elements are derived by analytically integrating the spherical harmonics series recursively. Compared with numerical integration by Gaussian quadrature, the newly developed analytical algorithm for interface traction integration is highly beneficial in terms of both computation precision and speed.


Introduction
It is well known that a dislocation approaching a free surface experiences an attractive force, so called image force [1]. In multiphase materials or single-phase elastically anisotropic polycrystalline materials, a dislocation near a phase or grain boundary is subject to a similar force due to the change in elasticity across the interface, which can be either attractive or repulsive [2][3][4]. As the magnitudes of these forces are inversely proportional to the dislocation-interface distance, the effect of free surfaces on dislocation motion and multiplication becomes significant in submicrometre-sized crystals [5][6][7][8], and the elastic interactions between dislocations and phase/grain boundaries have a stronger impact on the mechanical properties of the nanostructured materials [9][10][11][12] than those of their coarse-grained counterparts.
For free surfaces and phase/grain boundaries of elastically anisotropic half spaces, such virtual forces (negative energy gradients) can be determined through the image force theorem of Barnett and Lothe [13][14][15][16][17]. For finite domains of more complex shapes, this problem of elastic interactions between dislocations and interfaces can only be solved through the coupling of a dislocation dynamics (DD) code and a finite element (FE) or boundary element (BE) solver [18][19][20][21][22][23][24][25][26][27][28]. In such cases, the stress field of the dislocation in an infinite elastically homogeneous medium is used to calculate the traction on the surface bounding the elastic solid, and then correction fields are added to impose the proper boundary conditions on the domain. For free surfaces this results in an imposition of an equal and opposite surface traction such that the net result is a zero traction condition, and for phase/grain boundaries there is a traction balance and displacement continuity condition that must be imposed. This work is focused on the first part of DD-FE/BE simulations of dislocation-interface interactions, i.e. determination of forces on interface elements due to tractions imposed by dislocations in an infinite elastically anisotropic medium.
Interface tractions exerted by dislocation stress fields must be integrated into nodal forces in a FE or BE solver. These nodal forces are surface integrals of the traction field T (force per unit area) over the individual interface FE or BE elements [20,24,28], where F (n) is the force on a FE or BE node n of an interface element S exerted by dislocation stress field σ (x), n is the interface normal, and N (n) s (x) are the FE or BE shape functions. In this paper, we will refer to these nodal forces as traction forces. The traction force on a FE or BE node is analogous to the interaction force due to dislocation stress field on a dislocation node in nodal based (one-dimensional FE) dislocation dynamics models [20,[29][30][31]. The dislocation interaction forces are line integrals of the Peach-Koehler force f P K (force per unit length) along the individual dislocation segments [20,31], where f (n) is the force on a node n of a dislocation segment L exerted by dislocation stress field σ (x), with b, t, and N (n) L (x) being the Burgers vector, line direction, and shape function of the dislocation segment, respectively. While the dislocation interaction force calculation is the core of a DD model, the interface traction force evaluation is the main link in DD-FE/BE finite domain simulations. However, the compromise between accuracy and efficiency of the numerical integrations of dislocation interaction forces and surface traction forces has been the bottleneck for large-scale DD simulations [31] and DD-FE/BE simulations of the elastic interactions between dislocations and free surfaces [20,24,[32][33][34]. Due to the limitations of numerical integrations, an alternative analytical integration of interface traction forces is highly desirable.
For isotropic elastic media, the stress field of a dislocation loop can be calculated analytically through line integrations of the derivatives of Green's function along piecewise straight dislocation segments [35,36]. Arsenlis et al. [31] developed analytical expressions for dislocation interaction force calculations that involves double integrations along the individual pairs of interacting dislocation segments, i.e. the first line integration along the source segment to get its stress field and the second line integration along the receiving segment to obtain the interaction force. These analytical expressions for dislocation interaction force calculations have already been used in many DD simulations, e.g. discovery of ternary dislocation junctions in body-centered cubic metals [37], interpretation of the size-dependent strength for micrometersized crystals [38], determination of low-angle grain boundary penetration resistances [39,40], observation of strain localization via defect-free channels in highly irradiated materials [41], and revealing the mechanisms of dynamic recovery during high temperature creep of single-crystal superalloys [42]. Queyreau et al. [43] have recently formulated analytical expressions to calculate surface traction forces induced by stress field of a dislocation in isotropic elastic media for rectangular surface elements, which give the precise solutions to triple integrals, i.e. one integral along the dislocation segments for the stress field and then a double integral over the surface element to obtain the traction force.
In anisotropic elasticity theory of dislocations, the analytical expression for the stress field of an arbitrary dislocation loop does not exist [44,45]. The dislocation stress field calculation based on Stroh's formulism [35,46] requires numerical integrations [47,48] or solving an eigenvalue problem for a six by six matrix [48] to obtain the associated matrices and their angular derivatives. The dislocation stress field computations using Mura's formula [49] has relied on direct numerical integrations of the derivatives of Green's function [50]. Recently, Aubry and Arsenlis [51] used a truncated spherical harmonics expansion to approximate the derivatives of Green's function, and have analytically integrated the spherical harmonics series to calculate dislocation stress field (single integrals) and interaction force (double integrals) for straight dislocation segments. In this work, we use the same spherical harmonics expansion of the derivatives of Green's function formulated in the previous work [51], and analytically integrate the associated spherical harmonics series to determine the interface traction force (triple integrals) in anisotropic elastic media for quadrilateral surface elements.

Method
In a homogeneous infinite linear elastic solid, the stress field of a dislocation loop can be expressed in terms of a contour integral along the loop [49], where C ijkl is the elastic stiffness tensor, is the permutation tensor, b is the Burgers vector of the dislocation loop, and ∂G vp /∂x d is the derivative of the Green's function G vp (x − x ), which is defined as the displacement in the x v -direction at point x in response to a unit point force in the x p -direction applied at point x . The Green's function in an anisotropic elastic medium has been obtained as a single integral [44], where M −1 vp (ξ) = vsm prw (ξξ) sr (ξξ) mw 2 lgn (ξξ) 1l (ξξ) 2g (ξξ) 3n , with the notation (ξξ) ij = ξ k C kijl ξ l . R is the norm of vector R = x − x , i.e. R = R . T is the direction of vector R, i.e. T = R/R. ξ is a unit vector that varies in the plane ξ · T = 0,ê x andê y are two orthogonal unit vectors in the same plane, and the angle betweenê x and ξ is ψ, Fig. 1. The corresponding integral expression for the derivative of Green's function is [52], where . See the overview of Bacon et al. [44] for more details.
The derivative of the Green's function is a product of a part depending only on 1/R 2 and an angular part g depending only on the direction T where (θ, φ) are the spherical coordinates of T .
There is no analytical expression for g vpd , but the function g vpd (T ) is suitable for decomposition in spherical harmonics.
The expansion coefficients g lm are independent of T (θ, φ), and are defined as The spherical harmonics Y m l are defined as the complex functions where P m l are the associated Legendre polynomials. To be consistent with the definition of elastic stiffness tensor C ijkl , we rewrite Y m l in the Cartesian coordinate system (e 1 , e 2 , e 3 ), i.e. x = T · e 1 , y = T · e 2 , and z = T · e 3 , in the form of The function g can then be defined as Using the spherical harmonics series expansion of its angular part and defining e 12 = e 1 +ie 2 with x + iy = R R · e 12 and z = R R · e 3 , the derivative of the Green's function can be evaluated by This definition involves a quotient of terms that are a function of R, which depends only on two variables m and l − 2k. It can be simplified further to obtain where S qm vpd is a sum of products composed of Q m l (k) and g lm vpd . The Green's function and its derivatives depend only on odd powers of 1/R. This property means that in the spherical harmonics expansion, the non-zero terms correspond to the odd powers of 1/R.
Within linear elasticity, the total stress field of a dislocation network is a superposition of the stress fields of individual dislocation segments in the network. For a straight dislocation segment linking two nodes at Cartesian coordinates x 1 and x 2 , respectively, its stress field in an anisotropic elastic medium can then be expressed by Consider a quadrilateral element delimited by its Cartesian coordinates x 3 , x 4 , x 5 , and x 6 as shown in Fig. 2, where x is the coordinate that spans within the element, and x is the coordinate that spans along the dislocation segment delimited by x 1 and x 2 . The vector R = x − x can be rewritten as Using n = p×q p×q and dS = p × q drds, with a linear shape function in the form of a four-term polynomial N (n) (x) = a 0 + a 1 r + a 2 s + a 3 rs, the traction force F (n) can be expressed as Combining Eq. (14) and Eq. (15) with dx = −tdy, the interface nodal force due to dislocation traction in an anisotropic elastic medium can be determined through The unknown parts of Eq. (16) are the following series of integrals: x 3 x 4 x 5 x y r s Figure 2: Geometry and associated variables for interface traction integration where i + j = k − 2, and the common part of the integrands will hereafter be referred as I ijk , i.e.
The main effort of this work lies in solving these series of triple integrals. We find that with the analytical solutions of a few seed integrals, all other integrals in the spherical harmonics expansion can be calculated analytically through recurrence relations. In Fig. 3, the required triple integrals for the traction force calculation are indicated with a blue color, and the rest triple integrals are necessary for the recursive integration to reach the required triple integrals of the same expansion order q and the next expansion order q + 1. As illustrated by the blue dashed lines in Fig. 3, the double integrals at the intermediate levels are needed for the calculations of the triple integrals at the highest levels, and the single integrals at the lowest levels must be used to calculate the double integrals at the intermediate levels.
The required recurrence relations are constructed using the following partial derivatives: where α = t · e 12 , β = p · e 12 , γ = q · e 12 , φ = t · e 3 , θ = p · e 3 , and ψ = q · e 3 . The first three recurrence relations Eq. (26) where H y ijk = The single integrals J y ijk , J r ijk , and J s ijk must be formerly calculated using the third set of recurrence relations given below.
The first two recurrence relations Eq. (46) and Eq. (47) for the single integrals J y ijk are obtained by integration over y of the partial derivative Eq. (17) and multiplying both sides of the integration equation by α and φ, respectively. The third recurrence relation for the single integrals J y ijk Eq. (48) is obtained by integration over y of the partial derivative Eq. (23). Analogously, the recurrence relations for the single integrals J r ijk (from Eq. (49) to Eq. (51)) and J s ijk (from Eq. (52) to Eq. (54)) are obtained using the corresponding partial derivatives.
As illustrated in Fig. 3, a number of seed integrals have to be first calculated before starting the recursive integrations. The seed single integrals J y 00−1 , J r 00−1 , and J s 00−1 can be directly calculated using the corresponding analytical solutions Eq. (55), Eq. (56), and Eq. (57), respectively.
The rest single integrals are all calculated using the recurrence relations Eq. (46)  H 003 = 2 F 003 = 2 E 003 = 2 and H r ijk+2 from the available double integrals H ijk and single integrals J y ijk and J r ijk at a lower k index, see Fig. 3. The last two recurrence relations Eq. (34) and Eq. (35) are used to calculate the double integrals H ijk at higher i and j indices, respectively. The seed triple integral K 001 can be directly calculated from the seed double integrals H 001 , F 001 , and E 001 using Eq. (64), which is obtained by successive integrations over y, r, and s of the partial derivatives Eq. (20), Eq. (21), and Eq. (22), and then summation of these three integration equations.
The remaining triple integrals are all calculated using the recurrence relations Eq. (26) to Eq. (30). The first three recurrences Eq. (26) to Eq. (28) are used to calculated the double integrals K y ijk+2 , K r ijk+2 , and K s ijk+2 from the available triple integrals K ijk and the double integrals H ijk , F ijk , and E ijk at a lower k index, see Fig. 3. The last two recurrence relations Eq. (29) and Eq. (30) are used to calculate the triple integrals K ijk at higher i and j indices, respectively.
The triple integrals K rs ijk can be directly calculated from K ijk , K r ijk , H r ijk , F ijk , and E r ijk using the relations Eq. (65) to Eq. (67). The first equation is obtained by successive integrations over y, r, and s of the partial derivative Eq. (21). The last two equations are obtained by multiplying the partial differential equations Eq. (17) and Eq. (19) by r on both sides, and successive integrations over y, r, and s of these two partial differential equations, respectively.

Results
To evaluate the accuracy and efficiency of our analytical traction force calculation, we perform q max (spherical harmonics expansion order) convergence tests, and compare with Gaussian quadrature numerical integrations. The infinite series of spherical harmonics expansions must be truncated in practice. Figure 4 shows how the relative error of the traction force calculation evolves as the spherical harmonics expansion order q max increases. Similar to the dislocation interaction force calculations of Aubry and Arsenlis [51], the traction force calculations converge faster for materials with lower elastic anisotropy ratios, which means high orders of spherical harmonics expansion are needed for materials with high elastic anisotropy.
How the computation cost increases with the order of spherical harmonics expansion is presented in Fig. 5. The cost of our traction force calculation grows quadratically as the stress field calculation using the analytical expressions given in Ref. [51]. For the traction force calculation, using q max = 20 is twenty times more expensive than using q max = 1. The computation cost ratio between the traction force and stress field calculations is rather insensitive to the spherical harmonics expansion order, and changes from seven to six when q max increases from one to twenty.   We now compare our analytical interface traction force integrations with Gaussian quadrature numerical integrations using the analytical stress field expressions of Ref. [51]. As both the traction force and stress field calculations use the same spherical harmonics expansions of the derivatives of Green's function, the comparison of analytical and numerical integrations for a fixed q max is roughly the same when q max is changed from one to twenty. Figure 6 depicts how the relative error between analytical and numerical integrations varies as the number of Gauss quadrature points increases in the numerical integrations. Similar to the traction force calculations of Queyreau et al. [43] for isotropic elastic media, the relative error decreases faster for larger dislocation-interface distances as the number of Gauss points increases in the numerical integrations. While the analytical solution is used as a reference to assess the error in numerical integrations, the comparison with numerical integrations using a large number of Gauss points can verify the correct implementation of the analytical traction force integration. Figure 7 illustrates the computation cost comparison of analytical and numerical traction force integrations. The analytical integration becomes more efficient when the number of Gauss points exceeds eight in the numerical integration. With eight Gauss points, the error of the numerical integration is above one percent for a large dislocation-interface distance of 1000 b, and above ten percent for smaller dislocation-interface distances of 100 b and 10 b, Fig. 6. Keep in mind that the traction force calculation is only for one dislocation segment and one interface element, and the computation error can escalate in DD-FE/BE simulations with larger numbers of dislocation segments and interface elements over many correlated time steps. As proposed by Weygand et al. [20], a minimum of one hundred integration points has to be used for numerical integrations of surface traction forces. As shown in Fig. 7, compared with a numerical integra- tion using one hundred integration points, our analytical integration is more than one order of magnitude faster in speed.
Despite the obvious advantages of the analytical interface traction force integration, the algorithm implementation into a specific DD-FE/BE code can bring a considerable overhead to a short research project. In such a case, using simpler analytical stress field expressions and standard Gaussian quadrature numerical integrations may be preferred. Figure 8 shows how the computation error and cost of numerical integration varies with the number of Gauss points and the order of spherical harmonics expansion used in the stress field expressions. The numerical integration error is evaluated with respect to an analytical traction force integration using q max = 20 for Ni with a dislocation-interface distance of 100 b. As can be seen in Fig. 8, there are optimized combinations of integration point number and spherical harmonics expansion order to achieve a desired level of accuracy, and it is often more efficient to increase the order of spherical harmonics expansion than the number of Gauss integration points.

Concluding remarks
Using spherical harmonics expansions of the derivatives of Green's function, we constructed the expressions for the interface traction force exerted by dislocation stress field in anisotropic elastic media, and develop hierarchical recurrence relations to integrate the spherical harmonics series to calculate the traction force. It is found that all the triples integrals associated with the spherical harmonics are functions of a few analytically solvable seed integrals. Compared with numerical integrations of the traction force, our analytical integrations have gained substantially  in computation precision and speed. This development of analytical interface traction integrations can impart accuracy and efficiency to DD-FE/BE simulations of the elastic interactions between dislocations and interfaces in general elastically anisotropic crystalline materials.