Numerical methodology for treating static and dynamic dislocation problems near a free surface

Simulation of dislocation dynamics enables researchers and scientists to explore the plastic behavior of crystalline materials under loading. Analytic solutions for the stress field due to a linear dislocation segment near a free surface are case-specific, e.g. dealing with either a horizontal segment or a vertical segment, and therefore hard to implement in time-dependent dislocation dynamics simulations as different dislocation segments could be differently oriented. This article presents a generalized numerical framework to find the stress field beneath a free surface due to the presence of a dislocation segment. The framework can be expanded to non-flat surfaces. Also, three-dimensional discrete dislocation dynamics simulations are performed here, which clearly show the effect of free surfaces on the flow stress of a material.


Introduction
Simulation of dislocation dynamics enables researches and scientists to explore the phenomena that happen during plastic deformation, which may not be observed in experimental work. Simulation of dislocation dynamics satisfying the physical boundary conditions in a 3D representative volume element is always a challenge.
A dislocation segment causes stress in a crystalline material. Field stresses derived for an infinite medium are considered as classical problems in dislocation theory [1,2]. The use of stress formulation derived for infinite medium needs an additional term to satisfy the physical boundary conditions in the cases of the semi-infinite and finite media. To derive/evaluate this additional correction term, researchers and scientists implemented different techniques in their work.
Yoffe [3] derived the elastic fields of a dislocation meeting a surface in an angle for an arbitrary choice of Burgers vector. Baštecká [4] formulated the field stress due to a pure edge circular dislocation loop near a free surface with Burgers vector normal to that free surface. Bacon and Groves [5,6] showed the field displacements due to an infinitesimal dislocation loop of arbitrary shape, orientation, and Burgers vector in a semi-infinite isotropic medium. Maurissen and Capella [7,8] derived the field stress correction terms of a dislocation segment parallel and perpendicular to a free surface in a semi-infinite elastic medium. Comninou and Dundurs [9] presented the formulations of the elastic field of an angular dislocation segment in isotropic half-space. For an anisotropic medium Lothe et al [10] derived an integral form of field stress for the case of a dislocation terminating at the free surface of an anisotropic half-space. Gosling and Willis [11] expressed the stresses due to an arbitrary dislocation in a semi-infinite medium as a line integral along the dislocation. Devincre [12] developed the expression of three-dimensional stress field for a straight dislocation segment in a global coordinate system. Also prior to that, Hirth and Lothe [2] developed the stress field of a linear dislocation segment in a segment-attached coordinate system.
Kubin et al [13] introduced the framework of 3D dislocation dynamics simulation in a finite medium. Later Canova et al [14] and Hartmaier et al [15] implemented different field correction methods in discrete dislocation dynamics simulations to incorporate the free surface effect. In a single crystal, El-Azab et al [16,17] formulated traction and dislocation flux boundary conditions. Khraishi et al developed several methods for the treatment of traction-free boundary conditions in 3D dislocation dynamics [18][19][20]. In these referenced works, the dislocation field stress formulation developed in an infinite medium was augmented by image stress and surface correction terms, i.e. additional stress terms to satisfy the physical boundary conditions. In a recent article, Po et al [21] reviewed the nodal discrete dislocation methods for the investigation of plastic behavior in crystals with free surfaces.
In the works mentioned earlier [18][19][20], correction terms to mitigate the surface traction was developed locally and then transformed into a global coordinate system to compute the stress field due to the presence of a dislocation segment. Moreover, for each dislocation segment, the works just-mentioned required the placing of image dislocations in all six directions (across six free surfaces), see figure 1. In discrete dislocation dynamics simulations and during the evolution of dislocation segments, this reflection of image segments requires a substantial computational effort. For example, if the computational box includes 100 segments then the reflected image segments that have to be accounted for will be 600! In this article, the authors present a generalized numerical framework to tackle traction-free boundary conditions using generally-prismatic rectangular dislocation loops padding the free surface. All the stress and correction terms are computed in the global coordinate system to save the effort of first-rank, and second rank, tensor local-to-global transformations and vice versa (as required in the case of figure 1). In the Theory section, the theoretical formulation for the problem of nullifying tractions on free surfaces is established. In the Numerical Considerations section, different numerical implementations particular to this problem are discussed. In the Results section, both static and dynamic computations are shown. The static for verification purposes against a known analytical solution involving an infinite plane, and the dynamic using dislocation dynamics simulations capable of modeling constant strain-rate loading.

Theory
A dislocation is a stress source in a crystalline material. Stresses in a field point P due to M dislocations in a medium, subject to linear elasticity, is written as Where s P represents the total stress tensor at field point P due to stresses caused by the M dislocations and s i is the stress tensor at point P caused by dislocation i. Each σ is a 3×3 matrix of stress components. The field stress for a straight dislocation segment AB in an infinite medium is derived in [12] as, where ( )  s ab r AB is described in [12] and not reproduced here for brevity. At the free surfaces, the traction vector   s = T non each surface point must be zero which cannot be ensured from equation (2). In fact, it is violated since equation (2) was developed for an infinite medium.  n is a unit vector perpendicular to the free surface and σ is the total stress tensor at the point. To meet this physical boundary condition, we need to add a correction term and rewrite equation (1) as Where, σ corr is the required stress correction term to ensure zero traction on each field point of the free surface. The use of a correction term as in equation (3), as an approach to satisfy a surface boundary condition, was used in the previously mentioned references.
Khraishi et al [18,19] introduced a collocation point method to derive σ corr . In the collocation point method, it is sufficient to satisfy the boundary conditions at the distributed collocation points only as it is implied that an infinite number of collocation points will produce the exact analytic solution for the stress field if one existed.
In the collocation point method described in [18,19] σ corr is formulated using a local coordinate system i.e. the free surface is always considered as the local xy plane. Computed σ corr in this method needs to be secondranked tensor transformed to sum with the stress from the crystal dislocation segments computed in the global coordinate system (using [12]). Also, first-ranked tensor transformation is needed for coordinates of points (like points A and B in figure 2). In this article, we show a method to compute σ corr in a global coordinate system so no further transformation is needed to compute equation (3).
For a dislocation segment inside a finite three-dimensional medium and beneath a free surface as shown in figure 3, we pad N generally-prismatic fictitious (i.e. non-crystalline) rectangular dislocation loops on the free surfaces where the collocation points are located at the center of these fictitious dislocation loops (see figure 4).
As we compute the traction   s = T nfor the surface points, we find that the condition s s s = = = 0 xz yz zz must exist to ensure zero traction. Here, these three stress components are annulled on a set of N collocation points representing the generally-prismatic loops padding the surface. Mathematically, the problem is formulated as follows:   An examination of equation (2) shows that it is linear in the Burgers vector components b , x b y and b . z This allows the following to be inserted into equation (4): and ab K z have been derived from the work of Devincre [12] and provided all of them for the first time in the Appendix.
When inserting the last equation in the equation just before it, the following system of equations ensues: The above system represents a system of 3N equations, which is more time consuming to solve than the previous system of N equations when an image segment was used in the last method [18]. The unknowns in this system, equation (7), are the Burgers vector components of the N loops. The last system of equations can be written in expanded form to show more clearly the ensuing kernel matrix:  Here the size of the kernel matrix is 3N×3N which equals N 9 . 2 However and from dislocation theory [18,19], the kernels for s xz and s yz associated with the b z component of the Burgers vector are zero in the plane of the loop, i.e. the S area in figure 4 where = z 0. Moreover, from the same theory, the kernels for s zz associated with the b x and b y components of the Burgers vector are also zero. Also, any non-zero kernels on the surface where z=0 can be written in simplified forms to speed up computations (see [18,19]). All of this means that the last system of equations can be written alternatively as: The above represents two sets of decoupled systems of equations; one is N 2 equations and the other is N equations. The size of the kernel matrices in these two systems isŃ N 2 2 andŃ N, respectively, which equals N 5 . 2 These two systems can be solved much faster than the system of N 3 linear equations (e.g. (8)). This can only be achieved by recognizing the de-coupling between the prismatic and the shear components of the Burgers vectors in the plane of the dislocation loop i.e. the prismatic component of the Burgers vector b z only has a non-zero normal stress s zz in the loop plane, and the shear components of the Burgers vector b x and b y only have a non-zero shearing stress s xz and s yz in the loop plane. Based on the above, equation (9) can be written in expanded form as: And Where, s ab i is the stress due to real dislocation segments in the crystal and s ab  j P is the stress contribution by the generally-prismatic dislocation loops j at point P.

Numerical considerations
Equations (10) and (11)  F . The Kernel matrix is associated with the fictitious surface loops only, whereas the forcing vector is associated with real crystal dislocations or dislocation segments. The dimension of the kernel matrix depends on the areal loop density on a surface, and the elements of the kernel matrix are computed only taking into account the surface loops and not the dislocations inside the crystal. Elements of the forcing vector are computed at the collocation points, by only taking into account the dislocations or the dislocation segments inside the crystal. So, for the cases where more dislocations or dislocation segments are present, the Kernel matrices do not change, but the forcing vector does. In a simulation of dislocation dynamics, the Kernel matrix needs to be computed only once at the start of the simulation, but the forcing vector needs to be updated with every time step. The Burgers vectors of the fictitious surface loops thus change with every iteration depending on the positions of the dislocations or the dislocation segments relative to the surface. In a three-dimensional crystal volume, it is a must to ensure zero traction on all the free surfaces. For example, in a regular cube, traction-free boundary conditions need to be established on all the six cube surfaces. The stress at any material point is computed by summing the stresses from the dislocations in the crystal and all the fictitious loops from the six surfaces.
For the uniform surface mesh in figure 4, the stress from loop i will dominate the state-of-stress computed at collocation point i, which makes the Kernel matrix [ ] K diagonally dominant (i.e. | | | | å ¹  k k ii j i ij ). Since matrix [ ] K is diagonally dominant, it is possible to utilize an iterative method to solve the ensuing systems of linear equations (equations (10) and (11)). Here a Gauss-Seidel iterative method [22] was adopted in conjunction with a relaxation parameter λ (to hasten the computations) equal to 1.35. This method is faster than the regular Gaussian elimination or any decomposition techniques [22]. For large matrices, any decomposition technique is complex and computationally expensive.
Moreover, for a curved dislocation, we divide it into linear segments and calculate the forcing vector using equation (2) for each segment. For multiple segments, their contribution is added together into the right-hand side of equations (10) and (11). This discretization of the dislocation is still able to capture the correct physics of the problem. Again, this method can be readily implemented into the dislocation dynamics code by Khraishi et al [18,19] and by Zbib et al [23]. Figure 5 shows the time savings in solving the linear system of equation (8) versus the two systems of equations (10) and (11). The figure shows a significant time improvement of approximately 14% for one surface mesh of 2500 (50 50) loops in one iteration. This time improvement is significant because, in a dislocation dynamics simulation, a hefty amount of time is saved over many time steps.
As a summary, to implement the above method, the following steps need to be executed 1. Create the surface mesh, e.g. pad the rectangular loops on the free surfaces and determine the collocation points.
2. Calculate the [ ] K matrix using the equations presented in the appendix.
3. Find the forcing vector [ ] F using equation (2), adding contribution from as many existing segments. (10) and (11) to find the unknown Burgers vectors of the generally-prismatic loops on the free surfaces. 5. Compute stresses at any field point using equation (12).  (12) and (11). The N 9 2 bar represents the time to solve equation (8) and the N 5 2 bar represents the time to solve both equations (10) and (11).

Results
First, the current method is verified for the case of a horizontal subsurface linear segment AB parallel to the surface, see figure 6, as this case has an analytical solution [7]. The dimensions of the free surface are 20,000b × 20,000b, where b is the magnitude of the Burgers vector  b . The segment is located at a z-depth of b 1, 000 . Segment AB has a length = L b 100 . The coordinates for A and B are exactly ( b 100 , 0, 0) and (0, 0, 0), respectively. The dislocation segment is placed at the middle of the finite volume so that the numerical solution is comparable to the analytical solution. The Burgers vector of the segment is set to ( , , 3 ) for simplicity and Poisson's ratio n = 0.383. The stress field along a line parallel to the x-axis ( = y 0), with a z-depth of b 400 from the free surface, was computed. Figure 7 shows a comparison between the analytical solution by Maurissen and Capella [7] and the numerical solution from this work for an increasing number of surface loops. We see a gradual improvement in the numerical solution from, figures 7(a) to (f), as the number of loops or collocation points increases. We also observe that the numerical solution stabilizes after a certain number of surface loops (about40 40), which can be explained via the Saint-Venant's principle. Saint-Venant's principle states that agreement between an exact solution and an approximate one, but equipollent, at a boundary will be achieved for field points that lie more than a 'characteristic' distance away from the boundary. The characteristic distance, in this case, is the average separation distance between collocation points where the boundary condition is enforced (i.e. where the problem is solved). As we increase the surface loops, the characteristic distance becomes smaller, and the numerical solution has more matching to the analytical solution. Therefore, for field points closer to the surface than the characteristic length, the solution is not accurate and may exhibit oscillatory behavior with a wavelength equal to the characteristic length. It is worth mentioning here that although the solution might be oscillatory in this case, it still oscillates about a mean (obtained using least-squares curve fitting, for example) equal to the correct solution.
In the case of 3D discrete dynamics simulations, the dislocation segments are not always below the center of the free surfaces as their position evolves over time. In this case, the analytical solution in [7] is not expected to match the results of the numerical method presented here. To illustrate this point, let's consider a sub-surface segment near the edge of a free surface as shown in figure 8(a). In this case, the analytical solution ( figure 8(b)) does not match the numerical solution which is the proper solution in this case.
In addition to above static computations, we incorporated the current numerical method into the dislocation dynamics (DD) code by Zbib et al [23] and Khraishi et al [18,19]. From classical dislocation theory, the stress field of a dislocation can be calculated accurately beyond its core distance of approximately 5b [24]. With this limit from elasticity theory, the Peach-Koehler force on a sub-surface segment cannot be determined within a core distance, or depth, from the surface. This is in harmony with other DD calculations as the overlapping of a subsurface dislocation segment with the core of a surface dislocation loop could trigger numerical problems. This is not a serious limitation by any means because the force acting to pull the segment towards the surface at these small depths is tremendously high and causes rapid vanishing of such segments into the surface. Moreover, since the plastic strain calculation from the computational cell is homogenized over the cell's volume, any inaccuracies in quantifying the motion of dislocation segments (e.g. occurring from not capturing image stresses, i.e. stress correction terms, accurately in the immediate vicinity of the free surface) will have practically no appreciable impact on such calculation.
Here, we simulate a constant strain-rate experiment using a 3D crystal made of an isotropic material with dimensions of´b b b 10, 000 10, 000 10, 000 . Here, density = -2, 710 kg m , 3 Poisson's ratio n = 0.33, and  surface loops, the loops actually cause the FR source to bow critically earlier than before, i.e. at an earlier applied strain amount, as the stress needed to cause this flowing is reduced thanks to a pulling force onto the dislocation coming from the surface loops. This pulling force accelerates the movement of the bowed source as well as the vanishing of the dislocation into the surface.  For this particular problem setup, there is a difference of about 40% between the two cases (loops versus no loops). This percent difference will change for different problems. For example, the image stress effect is expected to be lower for problems involving a lower external surface area-to-volume ratio for the computational domain. This percentage will also change with the exact initial position of the dislocation source with respect to the interior of the computational domain (i.e. whether the source is initially close or far away from the free surfaces). See the discussions in [18,19].
In figure 10 we have superimposed all the stress-strain diagrams to accentuate the effect of surface loop density over the plastic behavior of the material. We observe no significant difference in the average stress-strain pattern for the different surface loop density, i.e. they all exhibit the same flow stress and flat plastic region onaverage. This result allows the use of lower surface loop density to save significant computational time.

Conclusion
The current study has presented an efficient mathematical formulation for treating the traction-free boundary condition of a finite three-dimensional crystal. A numerical framework was developed to implement the theory in dislocation dynamics simulation. The method was implemented for a mesh of rectangular/square surface dislocation loops for both the static and dynamic cases. The DD simulations showed that incorporating the image stresses results into a different flow stress for the crystal (a significant difference in this case). The simulations also showed that they are not sensitive to the areal mesh density. This present method can be extended to other regular polygonal surface loops. Moreover, it could be extended to adaptive surface meshing on planar surfaces or even be generalized to treat non-planar free surfaces. Future works are planned to address the last two statements.
xyy y z z y y z z y y y y y 2 Where, m is the isotropic material shear modulus, n is the Poisson's ratio and