A Mesh-Independent Brute-Force Approach for Traction-Free Corrections in Dislocation Problems

Simulation of dislocation dynamics opens the opportunity for researchers and scientists to observe in-depth many plastic deformation phenomena. In 2D or 3D media, modeling of physical boundary conditions accurately is one of the keys to the success of dislocation dynamics (DD) simulations. The scope of analytical solutions is restricted and applies to specific configurations only. But in dynamics simulations, the dislocations’ shape and orientation change over time thus limiting the use of analytical solutions. The authors of this article present a mesh-based generalized numerical approach based on the collocation point method. The method is applicable to any number of dislocations of any shape/orientation and to different computational domain shapes. Sev-eral verifications of the method are provided and successful implementation of the method in 3D DD simulations have been incorporated. Also, the effect of free surfaces on the Peach-Koehler force has been computed. Lastly, the effect of free surfaces on the flow stress of the material has been studied. The results clearly showed a higher force with increased closeness to the free surface and with increased dislocation segment length. The simulations’ results also show a softening effect on the flow stress results due to the effect of the free surfaces.


Introduction
The dislocation is a type of defect and a source of stress in crystalline materials.
Dislocation theory explains plastic deformation in a material. The stress field solution in an infinite medium due to dislocation is a fundamental problem of linear elastoplasticity [1] [2]. Eshelby and Stroh [3] extended this stress field solution in the case of a semi-infinite medium by adding a correction term. The correction term is derived by introducing image dislocations to treat the boundary condition of the free surface, which is analogous to the electrostatic problem of a line charge in a medium of non-homogeneous dielectric constant [4]. To derive this correction term for different boundary conditions, there have been many methods developed by researchers in the literature.
The earliest work dates back to 1961 by [5]. In this work, the problem of a terminating dislocation at a free surface was considered. Here, the dislocation was straight with any angle of incident or Burgers vector. In the case of a dislocation loop below a zero-traction surface, [6] found the stress solution in the case of a perpendicular Burgers vector to the surface. Later, [7] [8] derived the displacement field solution for an infinitesimal dislocation loop of arbitrary Burgers vector and orientation, and in a half-space, which can be extended to find the elastic fields of a finite dislocation loop. For a linear dislocation segment perpendicular and parallel to a free surface, [9] [10] derived the correction term in an isotropic medium that was half-infinite. Comninou and Dundurs [11] derived the elastic field for an angular dislocation in a half-medium. This solution can be extended for an angular dislocation ending at the free surface in conjunction with the solution presented by [5]. In the case of a dislocation terminating at the free surface of a half-space, [12] derived an integral form for the field stress in an anisotropic medium. Using a line integral tracing the dislocation, [13] determined the stress solution for a dislocation in a half-medium. For a linear dislocation segment in an isotropic material, [14] found the stress field using a global coordinate system. Before that, [2] found the same stress field but using a coordinate system which was fixed to the dislocation segment.
Some of these developments enabled researchers to simulate dislocation dynamics in 2D and 3D space. Dislocation dynamics (DD) simulations brought to light many phenomena during plastic deformation, which are difficult, or sometimes impossible, to observe during real experiments. For three-dimensional dislocation simulations, the problem of dislocation interaction with surfaces or boundaries is always challenging. Nowadays, DD simulations are considered a valuable tool for plasticity research.
The first group to model three-dimensional dislocation motion and interaction in a simulation box was [15]. The effects of free surfaces were addressed using field correction methods in dislocation dynamics [16] [17]. Dislocation flux and traction boundary conditions were developed by [18] [19]. An image stress tensor to capture the effect of free surfaces was used by [20]. The effect of free surfaces on void strengthening was studied by [21]. Jing et al. (2005) [22] also studied the effect of voids on dislocation motion. Dislocation theory and molecular dynamics were used by [23] for sub-surface dislocation loops. Lastly, the A. B. Siddique, T. Khraishi development of methods to treat zero-traction boundaries for employment in three-dimensional dislocation simulations were presented by [24] [25] [26]. In the last works by [24] [25] [26], the developed methods revolved around a combination of image stresses and additional surface terms to correct for the proper physical boundary conditions. The corrections were needed since the employed dislocation stress fields in the simulations were for an infinite domain, whereas the actual computation was done in a finite domain.
In the references indicated earlier [24] [25] [26], correction terms to satisfy the traction-free surfaces were generated for dislocation segments using a surface-attached coordinate systems. They were then transformed to a global coordinate system in order to calculate the stress evolution problem. More recent treatments for the effect of free surfaces has been provided by [27] for a flat surface in 3D space and for a curved one in 2D space [28].
In the proposed method for this article, all the stresses caused by the real and image dislocations/correction terms are calculated using a global coordinate system such that no first-rank or second-rank tensor transformations are required. Moreover, unlike previous collocation-point methods by [24] [25] [26] discussed above, the current method can utilize a structured or unstructured mesh on the free surface with any mesh element type (rectangular, triangular, etc.). This is for the purposes of satisfying zero traction condition on the free surface. The triangular mesh can, for example, save on computational time (see below). Below, the development of the correction terms is performed in the Theory section. Then the numerical model is established, along with any considerations, in the Numerical Considerations section. Lastly, dynamic and static dislocation problems are presented in the Results section. The static case verifies the solution against a known analytical result for a simple case. The dynamic case implements the method in a 3D dislocation dynamics simulation code that can generate a stress-strain curve.

Theory
A dislocation inside a crystalline material causes stresses in it. Stresses at any field point P due a linear dislocation segment AB in an unbounded space or medium was presented by [14] as, and ( ) AB r αβ σ  is described in [14].
where the material is isotropic with a modulus for shear μ and a Poisson ratio indicated by ν. In addition, Here, ( ) , , a b c    is called a scalar triple product and is given by: is called symmetric tensor operator and is given by: The total stress P σ , in linear elasticity, at P (any field point) due to a number of dislocations γ can be written as, Here, the stress of dislocation l is given by l P → σ . The stress second-rank tensor is a 3 × 3 matrix, i.e. it has 9 components of stress.
The stress field formulation in Equation (2) is not valid for a semi-infinite or finite medium since it does not ensure zero traction on the free surfaces. The traction vector T n =   σ must be zero on the free surfaces of a semi-infinite or finite medium. The unit normal vector to the surface is indicated by n  , and the stress state at a surface point is given by σ . To match this physical condition at any point M on the free surface, a correction term can be added by re-writing Equation where, M corr σ is the correction term for stress. This term ensures traction is null on any free surface point. This correction term is a spatial function of the dislocations and the field point. For example, any two field points in a material may have different correction terms depending on the spatial coordinates of the field points. Also, the correction terms change if the sub-surface dislocation changes position, say during a dynamics simulation.
The Peach-Koehler force at the center c of dislocation q as follows: where, l qc → σ is the stress tensor at the center of dislocation q due to dislocation l. q b  , q t  and qc corr σ are the Burgers vector, line sense and the correction term (evaluated at the center) of dislocation q, respectively. Also the total stress tensor at any field point P (not just surface points) is A collocation point method was provided in [25] and [26] to capture corr σ . In this method, a number of collocation points on the surface in question are used to enforce the physical boundary conditions. As more points cover the surface, this numerical method is slated to match, in the limit of infinite collocation points, the correct analytical solution for the problem (if one existed of course).
The method introduced in [25] and [26] uses a local coordinate system to formulate corr σ . Specifically, the surface-attached xy-plane is the free surface.
The correction terms in this method represent a second-rank tensor and need to be transformed as such since they need to be summed with the stresses from the crystal dislocations (which are computed with reference to a global coordinate system using [14]). In addition to above, the coordinates or position vectors for points A and B (see Figure 1 and Figure 2) also require transformation as first-rank tensors (from global coordinates to local or suface-attached coordinate system). The method presented herein avoids such first-rank and second-rank transformations, in order to find the correction terms corr σ , as a global system is used for all computations (such as Equation (7) or Equation (9)).
Consider loops can be rectangular or triangular in shape. As shown in Figure 3, the collocation points represent the centers of the dislocation loops.
The condition for the surface traction translates to 0 yx yy yz σ σ σ = = = for the surfaces in Figure 3. The annulment of these stresses here happens on a set of N collocation points. These N points represent the centers of dislocation loops (all equilibriated sources of stress) meshing the surface.
This condition on the surface can now be expressed as, where, l i → σ and j i → σ represent the stress tensor at surface collocation point i due to dislocation l and generally-prismatic surface loop j, respectively. After simplification, the problem is formulated by rewriting Equation (10) as follows: where j i αβ σ → is the stress computed at point i due to loop j. Also, l i αβ σ → is the αβ component of stress due to dislocation l and calculated at collocation point i on the surface.
For a dislocation loop like DEFG in Figure 3(a), the stress is calculated at any (collocation) point as the contribution of four segments in the loop. The calculation here being done with respect to a global coordinate system (using [14]): is the stress of dislocation segment (1) of loop j (shown as DEFG in Figure 3 generally-prismatic triangular dislocation loop (see Figure 3(b)) in the global coordinate system one can write, When one analyzes the segment stress solution by [14], see eqation (2), it is apparent that is linear with respect to the Burgers vector components ( x b , y b and z b ). Therefore, in Equation (11), one can then insert: Here, the stress ( ) , , j i x y z αβ σ → is due to Loop j acting on collocation point i. Also, the x-component of the Burgers vector is labelled here j x b for loop j.
Lastly, the kernel term associated with this Burgers vector component is labelled It is a field function evaluated at collocation point i. In this article, the kernels are computed using a "brute-force" method. This method is a straight-forward and generalized method to compute kernels in Equation (14). The desired kernels can be computed by canceling the other Burgers vector components in Equation (14). For example, x K αβ can be calculated by putting 0 (14). In general, Take Equation (14) and insert it into Equations (11). The result is a system of linear algebraic Equations (replacing Equations (11)): Equations (16) has 3N equations. Such system requires increased solution time than the previously mentioned system of N equations which required the suspension of an image segment in air (e.g. [25]). However, the current method is more straight forward to implement. It also can be expanded to different types of meshes and saves time on the tensor transformations. For Equations (16), the sought after solution vector is composed of the different Burgers vectors for the N loops. Equations (16) can be re-written in matrix form as: The matrix containing the kernel terms is labelled here the kernel matrix. Its size is 3N × 3N (=9N 2 ). The Burgers vectors for the padded surface loops can be found by solving Equation (17). Once the Burgers vectors for the loops are found, the correction term corr σ due to all the surface loops can be computed and substituted back into Equation (9) to find the stresses at any field point inside the medium.

Solving the System of Equations
Equation (17) can be expressed in matrix form as, The left-hand matrix [K] is termed the kernel matrix, B is the solution Burgers vectors of the padded generally-prismatic dislocation loops and F is termed the forcing vector. applies. Hence, to calculate the stress at any field point in the box, one needs to sum stresses from all surface loops on all surfaces plus the stresses rising from internal stress sources such as real crystal dislocations (Equation (9)).
Solving Equation (18) is the most crucial part of this method. For each collocation point (like point i in Figure 3), the state of stress from its own loop (loop i in this case) will dominate over stresses from other loops. This will make the Kernel matrix [K] to be dominant on its diagonal. In other words, ii Iterative methods, like the Gauss-Seidel method [29], can be employed for cases when the matrix is diagonally-dominant. Also, successive relaxation can be uti- In this study, the relaxation parameter λ is set equal to 1.35 (over-relaxation).
This iterative method works faster than decomposition methods like Gaussian elimination for example (refer to [29]). In general, decomposition techniques in linear algebra tend to be expensive computationally especially for large matrices.
In the case of curved crystal dislocations, each of them can be divided into linear segments. The stress from each segment can be found using Equation (1).
The right-side of Equations (17) can add up the contribution of all these segments. Such dislocation curve segmentation or discretization is routine practice in dislocation dynamics simulation codes (see for example [24] [25] [26] [30] [31]). The dislocation line discretization is capable of capturing the problem physics correctly. First, the current brute-force method with rectangular and triangular surface loops is verified. Consider a situation of a horizontally-situated linear dislocation segment AB that is below and parallel to a free surface (refer to Figure 5). This configuration has an analytical solution although the solution is for an infinite surface [9]. In this work, each side of the finite-length square surface used here is   (c) (d) Figure 6. Stresses compared to the analytical solution found in [9] as the number of rectangular surface dislocation loops is increased (see Figure 3(a)). MC refers to Maurissen and Capella [9] and CM refers to the current method.  Figure 7. Stresses compared to the analytical solution found in [9] as the number of triangular surface dislocation loops is increased (see Figure 3(b)). MC refers to Maurissen and Capella [9] and CM refers to the current method. Second, the current brute-force method with rectangular and triangular surface loops is verified for a vertical segment AB beneath the free surface (see Figure 8). This case is addressed analytically for an infinite surface dimensions by   [10]. Each side of the finite-length square surface used here is 20,000b. The depth of point B is 1000b. Segment AB has a length 100

Solution Time
For points A and B, the following exact coordinates were chosen ( ) 0, 0, 1,100b − and ( ) 0, 0, 1, 000b − , respectively. Also, b  for segment AB is chosen as Poisson ratio ν is 0.38. Consider a line, for field point calculations, beneath the surface and parallel to it at a depth of 400b (along 0 x = or the y-axis). Figure 6 and Figure 7 compare the analytical solution by [9], for a horizontal sub-surface dislocation segment, and the solution from this numerical method as we increase the surface mesh density for the rectangular and triangular surface loops. And Figure 9 and Figure 10  (c) (d) Figure 9. Stresses compared to the analytical solution found in [10] as the number of rectangular surface dislocation loops is increased (see Figure 3(a)). MC refers to Maurissen and Capella [10] and CM refers to the current method. (a) # dislocation loops 400; (b) # dislocation loops 900; (c) # dislocation loops 1600; (d) # dislocation loops 2500. Figure 10. Stresses compared to the analytical solution found in [10] as the number of triangular surface dislocation loops is increased (see Figure 3(b)). MC refers to Maurissen and Capella [10] and CM refers to the current method. (a) # dislocation loops 400; (b) # dislocation loops 900; (c) # dislocation loops 1600; (d) # dislocation loops 2500.
for a vertical sub-surface segment, and the solution from this numerical method as we increase the surface mesh density for the rectangular and triangular surface loops.
In the case of the horizontal sub-surface segment, the numerical solution approaches the analytical solution as the number of dislocation loops, i.e. collocation points, increase. Same goes to the vertical sub-surface segement. Saint-Venant's principle can be used to explain the stabilization of the numerical solution with increased number of surface loops (approximately 1600). According to this principle, an approximate solution (but equipollent) at a boundary and an exact one will show agreement at field points that are positioned from the boundary a "characteristic" distance away. In this work, the characteristic distance is the spacing, on average, between the surface collocation points (where the solution is enforced or the problem solved). Therefore, an increase in the number of surface loops will result in a smaller characteristic distance and thus more matching between the numerical and analytical solutions. The solution is not accurate for field points nearer to the surface than this characteristic distance. It may even exhibit oscillatory behavior in space with a wavelength equaling this distance. The mean value around the oscillations (e.g. by applying a least-square fit) will be equivalent to the correct solution.
The above was verification of the method for static cases. For dynamic cases, we use the paper's method into a three-dimensional dislocation dynamics (or DD) code (see the works of [24] [25] [26] [27] [30] [31]). Dislocation theory dictates that stress calculations cannot be made closer to the dislocation line than a core radius. Meyers and Chawla [32], for example, define the core radius as 5b (where b is the magnitude of the Burgers vector). Therefore, when calculating the Peach-Koehler (PK) force on a crystal dislocation, it has to be below the surface by at least a core distance. This is a known problem in dislocation dynamics codes as the overlapping of a dislocation core (say for the surface loops) with a field point causes numerical problems. However, this is not a serious issue as the PK force is so high that it quickly pulls the segment into the free surface to vanish it. Below, a figure will be shown of such force. Also, since the plastic strain computation in DD codes is homogenized over the volume of the computational cell, this small area sweep close to the surface contributes negligibly to the total strain value, i.e. the error in not fully capturing the right image stresses (which are the correction terms emanating from the surface dislocation loops in the current case) is so minimal that it does not impact the overall computational results. This was shown in these preliminaries works referenced.
In what follows, constant strain-rate experiments are simulated for a threedimensional crystal or cubic simulation box. The box has sides of 10,000b each. Each of the six box faces is padded with either rectangular or triangular generally prismatic surface dislocation loops (see Figure 2). The mesh density on the six free surfaces is also varied. For these simulations, the density = 2710 kg/m 3 Figure   11 and Figure 12, the stress and strain are the yz-components.  To study the effect of surface loop number, or equivalently loop density, on the stress-strain curves, refer to Figure 11 and to Figure 12. The former uses rectangular surface loops for the surface correction terms, and the latter uses triangular surface loops. An averaging of the output stress-strain curve, i.e. smoothing, is done using a cubic spline fit between the two linear portions of the curve (after linear regression of the plastic deformation part) and presented in a solid line for each case. The cases of surface dislocation loops are compared to the base case of no surface loop, i.e. no treatment of the zero-traction boundary condition on the surface. The case of no surface loops shows higher flow stress than the cases of surface loops. This stress-strain diagram is typical in DD computations for an unobstructed Frank-Read source. The initial bowing of the FR source occurs in an elastic fashion (i.e. the source will bounce back to its original dislocation shape once load is removed). However, after a critical bow radius or size, the source then bows plastically in the material, i.e. irreversibly. That is when the stress-strain diagram starts bending from its elastic region. It takes multiple bows or loops after the initial bowing for the stress-strain diagram to flatten, i.e. reach a flow stress value. When surface loops are present, they tend to pull on the dislocation source which causes it to bow earlier than the case of no surface loops. Meaningless applied strain is needed in the case of surface loops to reach critical bowing. This is due to the physical force that wants to pull on the dislocation and vanish it in the free surface to reduce the system's overall entropy, i.e. reduce the disorder in the crystal. This surface pull helps reduce the stress amount needed to cause this plastic flow and hence the observation in Figure 11 A. B. Siddique, T. Khraishi and to Figure 12.
If one measures the difference in flow stress for the above two cases of surface loops versus no surface loops, one finds this value approximately equal to 22%. The indicated percentage difference shall vary for different dislocation problems. If a problem has lower free surface area to volume ratio in the modeling cell then this percent amount should be lower due to lowered image stress effects. Another factor that will change this percentage is the exact initial position of the FR source inside the modeling cube/cell. Previous work has shown that the image stress effects depend on if the dislocation source is initially close to a free surface of the cube or far from it (i.e. positioned more internally in the cell). Discussions of this nature can be found in [25] and [26]. Figure 13 shows the effect of the free surfaces on the critical bowing of a Frank-Read source and the corresponding points on the stress-strain curves. Figure 13(b) is a zoomed-in figure from Figure 11(b). As is known [1], at the critical bowing moment, the dislocation becomes unstable and starts bowing quickly and irreversibly, i.e. plastic deformation ensues. From the simulations, the bowing dimension L 2 observed without surface effects is about 14 percent more than the bowing dimension L 1 when incorporating free surface effects. This is because the surface pulls the dislocation stronger as it gets closer to it. This causes the plastic flow to occur earlier than the case without free-surface effects (Figure 13(b)). As is observed in Figure 13, plastic flow started with a smaller radius of curvature for the dislocation source in the case where the free-surface effect is incorporated.  (Figure 14(b)). For Figure 14(b), the segment was centered with respect to the z-axis. As the dislocation moves close to the surface, the Peach Koehler force becomes stronger and pulls the dislocation towards it. So it is very important to treat the boundary condition properly to model this problem accurately. Figure 14(b) shows the z-component of the Peach Koehler force (acting at the center of the dislocation) as it varies with the dislocation segment length where the stress reaches a limiting value as the dislocation segment length is sufficiently high.

Conclusion
The current study presented a generalized numerical technique to properly enforce zero-traction condition on a free surface of a three-dimensional crystal containing dislocations. A numerical method was generated adopting this generalized formulation and termed here as the "brute-force" method for implementation in dislocation dynamics (DD) simulations. The method uses dislocation loops, as self-equilibrated sources of stress, meshing the free surface. We have shown that such meshing could involve triangular or rectangular elements, i.e. dislocation loops. The method was illustrated for dynamic and for static dislocation problems. This "brute-force" method can be adapted to model non-planar free surfaces too. The stress correction term can also be computed using an unstructured random triangular surface mesh, which enables the possibility of modeling more complicated non-planar free surfaces. Lastly, we have shown using dislocation dynamics that the flow stress in the model is considerably different in the cases of image stress incorporation or not. The areal density of the surface mesh, or the mesh type, was not a riding factor for determining the flow stress. Hence, a lower mesh density can be used to save computational time. Due to the flexibility of the current method on the mesh type, future work can potentially utilize the developed collocation-point method to solve problems with curved free surfaces.

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.