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

Abu Bakar Siddique, Tariq Khraishi

Mechanical Engineering Department, University of New Mexico, New Mexico, USA.

DOI: 10.4236/mnsms.2021.111001                

Abstract

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. Several 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.

Keywords

Free Surfaces, Dislocation Theory, Dislocation Dynamics, Peach-Koehler Force

1. 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 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.

2. 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,

${\sigma }_{\alpha \beta }^{AB}={\left({\sigma }_{\alpha \beta }^{AB}\right)}_{{\stackrel{\to }{r}}^{\prime }=OB}-{\left({\sigma }_{\alpha \beta }^{AB}\right)}_{{\stackrel{\to }{r}}^{\prime }=OA}$ (1)

and ${\sigma }_{\alpha \beta }^{AB}\left(\stackrel{\to }{r}\right)$ is described in [14].

$\begin{array}{l}{\sigma }_{\alpha \beta }^{AB}\left(\stackrel{\to }{r}\right)=\frac{\mu }{\pi Y^2}[{\left[\stackrel{\to }{b}\stackrel{\to }{Y}\stackrel{\to }{t}\right]}_{\alpha \beta }^s-\frac{1}{\left(1-\nu \right)}{\left[\stackrel{\to }{b}\stackrel{\to }{t}\stackrel{\to }{Y}\right]}_{\alpha \beta }^s\\ -\frac{\left(\stackrel{\to }{b},\stackrel{\to }{Y},\stackrel{\to }{t}\right)}{2\left(1-\nu \right)}\left[{\delta }_{\alpha \beta }+t_{\alpha }{t}_{\beta }+\frac{2}{Y^2}\left[{\rho }_{\alpha }{Y}_{\beta }+Y_{\alpha }{\rho }_{\beta }+\frac{L}{R}{Y}_{\alpha }{Y}_{\beta }\right]\right]]\end{array}$ (2)

where the material is isotropic with a modulus for shear μ and a Poisson ratio indicated by ν. In addition,

$\begin{array}{l}\stackrel{\to }{R}=\stackrel{\to }{r}-{\stackrel{\to }{r}}^{\text{'}}\\ L=\stackrel{\to }{R}\cdot \stackrel{\to }{t}\\ \stackrel{\to }{\rho }=\stackrel{\to }{R}-L\stackrel{\to }{t}\\ \stackrel{\to }{Y}=\stackrel{\to }{R}+R\stackrel{\to }{t}\end{array}$ (3)

Here, $\left(\stackrel{\to }{a},\stackrel{\to }{b},\stackrel{\to }{c}\right)$ is called a scalar triple product and is given by:

$\left(\stackrel{\to }{a},\stackrel{\to }{b},\stackrel{\to }{c}\right)=\left(\stackrel{\to }{a}\times \stackrel{\to }{b}\right)\cdot \stackrel{\to }{c}$ (4)

and ${\left[\stackrel{\to }{a}\stackrel{\to }{b}\stackrel{\to }{c}\right]}_{\alpha \beta }^s$ is called symmetric tensor operator and is given by:

${\left[\stackrel{\to }{a}\stackrel{\to }{b}\stackrel{\to }{c}\right]}_{\alpha \beta }^s=\frac{1}{2}\left({\left(\stackrel{\to }{a}\times \stackrel{\to }{b}\right)}_{\alpha }{c}_{\beta }+{\left(\stackrel{\to }{a}\times \stackrel{\to }{b}\right)}_{\beta }{c}_{\alpha }\right)$ (5)

The total stress ${\sigma }^P$, in linear elasticity, atP (any field point) due to a number of dislocations $\gamma$ can be written as,

${\sigma }^P={\displaystyle {\sum }_l^{\gamma }{\sigma }^{l\to P}}$ (6)

Here, the stress of dislocation l is given by ${\sigma }^{l\to 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 $\stackrel{\to }{T}=\sigma \stackrel{\to }{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 $\stackrel{\to }{n}$, and the stress state at a surface point is given by $\sigma$. To match this physical condition at any point M on the free surface, a correction term can be added by re-writing Equation (6):

${\sigma }^M={\displaystyle {\sum }_l^{\gamma }{\sigma }^{l\to M}}+{\sigma }_{corr}^M$ (7)

where, ${\sigma }_{corr}^M$ 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:

${\stackrel{\to }{F}}^q=\left(\left({\displaystyle {\sum }_{l\ne q}^{\gamma }{\sigma }^{l\to qc}}+{\sigma }_{corr}^{qc}\right){\stackrel{\to }{b}}^q\right)\times {\stackrel{\to }{t}}^q$ (8)

where, ${\sigma }^{l\to qc}$ is the stress tensor at the center of dislocation q due to dislocation l. ${\stackrel{\to }{b}}^q$, ${\stackrel{\to }{t}}^q$ and ${\sigma }_{corr}^{qc}$ 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

${\sigma }^P={\displaystyle {\sum }_l^{\gamma }{\sigma }^{l\to P}}+{\sigma }_{corr}^P$ (9)

A collocation point method was provided in [25] and [26] to capture ${\sigma }_{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 ${\sigma }_{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 ${\sigma }_{corr}$, as a global system is used for all computations (such as Equation (7) or Equation (9)).

Consider Figure 2. It shows a sub-surface dislocation (or dislocation segment more specifically) below a traction-free surface. The dislocation is in a 3D domain of finite dimensions. The surface is meshed or padded by elements. Each surface element stands for a dislocation loop. Each dislocation loop is called a generally-prismatic loop since its Burgers vector is angled with respect to the normal unit vector for the loop. Consider N such loops. These are not real dislocations, i.e. not crystal dislocations, but rather fictitious or mathematical dislocations used to solve the boundary-value problem at hand. These elements or

Figure 1. The figure defines the quantities used in calculating the σ_αβ stress component due to a linear dislocation segment labelled AB. The vector $\stackrel{\to }{r}$ represents the position vector of an arbitrary field point P. The point O is the origin of a generally-situated global Cartesian coordinate system (xyz-system). The line sense vector for segment AB is denoted by the vector $\stackrel{\to }{t}$.

(a) (b)

Figure 2. A cubic crystal used for the computation (showing surface loops) with segment AB shown inside of it. (a) Padded rectangular loops on the surfaces; (b) Padded triangular loops on the surfaces.

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 $\stackrel{\to }{T}=\sigma \stackrel{\to }{n}=\stackrel{\to }{0}$ for the surface traction translates to ${\sigma }_{yx}={\sigma }_{yy}={\sigma }_{yz}=0$ 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,

$\left({\displaystyle {\sum }_l^{\gamma }{\sigma }^{l\to i}}+{\displaystyle {\sum }_j^N{\sigma }^{j\to i}}\right)\stackrel{\to }{n}=\stackrel{\to }{0}$ (10)

where, ${\sigma }^{l\to i}$ and ${\sigma }^{j\to 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:

$\begin{array}{l}{\displaystyle {\sum }_{j=1}^N{\sigma }_{yx}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yx}^{l\to i}},i=1,2,\cdots ,N\\ {\displaystyle {\sum }_{j=1}^N{\sigma }_{yy}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yy}^{l\to i}},i=1,2,\cdots ,N\\ {\displaystyle {\sum }_{j=1}^N{\sigma }_{yz}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yz}^{l\to i}},i=1,2,\cdots ,N\end{array}$ (11)

where ${\sigma }_{\alpha \beta }^{j\to i}$ is the stress computed at point i due to loop j. Also, ${\sigma }_{\alpha \beta }^{l\to 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] ):

${\sigma }^{j\to i}={\sigma }^{j_{\left(1\right)}\to i}+{\sigma }^{j_{\left(2\right)}\to i}+{\sigma }^{j_{\left(3\right)}\to i}+{\sigma }^{j_{\left(4\right)}\to i}$ (12)

where ${\sigma }^{j_{\left(1\right)}\to i}$ is the stress of dislocation segment (1) of loop j (shown as DEFG in Figure 3(a)) evaluated at collocation point i, and so forth. Similarly, for a

(a) (b)

Figure 3. The figure shows a linear dislocation segment AB lying below a free surface. (a) This figure shows a meshing of the free surface with rectangular elements. In this method, each element is a dislocation loop with three Burgers vector components (b_x, b_y and b_z; see the inset). The area S covered by the mesh is where we are seeking to annul the traction in this problem. The inset also illustrates one of the dislocation loops, i.e. DEFG. (b) This figure is similar to (a) but the elements/dislocation loops here are triangular in nature. One such loop is shown in the inset along with its Burgers vector components, i.e. loop DEF. The center of each element is a collocation point. There is N collocation point. For example, point i in the sketch is a collocation point. (a) Rectangular surface loops; (b) Triangular surface loops.

generally-prismatic triangular dislocation loop (see Figure 3(b)) in the global coordinate system one can write,

${\sigma }^{j\to i}={\sigma }^{j_{\left(1\right)}\to i}+{\sigma }^{j_{\left(2\right)}\to i}+{\sigma }^{j_{\left(3\right)}\to i}$ (13)

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 ( $b_x$, $b_y$ and $b_z$ ). Therefore, in Equation (11), one can then insert:

${\sigma }_{\alpha \beta }^{j\to i}\left(x,y,z\right)=K_{x\alpha \beta }^{j\to i}\left(x,y,z\right)b_x^j+K_{y\alpha \beta }^{j\to i}\left(x,y,z\right)b_y^j+K_{z\alpha \beta }^{j\to i}\left(x,y,z\right)b_z^j$ (14)

Here, the stress ${\sigma }_{\alpha \beta }^{j\to i}\left(x,y,z\right)$ is due to Loop j acting on collocation point i. Also, the x-component of the Burgers vector is labelled here $b_x^j$ for loop j. Lastly, the kernel term associated with this Burgers vector component is labelled $K_{x\alpha \beta }^{j\to i}\left(x,y,z\right)$. 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, $K_{x\alpha \beta }$ can be calculated by putting $b_y=0$, $b_z=0$ and $b_x=1$ in Equation (14). In general,

$\begin{array}{l}{K}_{x\alpha \beta }^{j\to i}\left(x,y,z\right)={{\sigma }_{\alpha \beta }^{j\to i}\left(x,y,z\right)|}_{b_x^j=1,b_y^j=0,b_z^j=0}\\ K_{y\alpha \beta }^{j\to i}\left(x,y,z\right)={{\sigma }_{\alpha \beta }^{j\to i}\left(x,y,z\right)|}_{b_x^j=0,b_y^j=1,b_z^j=0}\\ K_{z\alpha \beta }^{j\to i}\left(x,y,z\right)={{\sigma }_{\alpha \beta }^{j\to i}\left(x,y,z\right)|}_{b_x^j=0,b_y^j=0,b_z^j=1}\end{array}$ (15)

Take Equation (14) and insert it into Equations (11). The result is a system of linear algebraic Equations (replacing Equations (11)):

$\begin{array}{l}{\displaystyle {\sum }_{j=1}^N{b}_x^j{K}_{xyx}^{j\to i}+b_y^j{K}_{yyx}^{j\to i}+b_z^j{K}_{zyx}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yx}^{l\to i}},i=1,2,\cdots ,N\\ {\displaystyle {\sum }_{j=1}^N{b}_x^j{K}_{xyy}^{j\to i}+b_y^j{K}_{yyy}^{j\to i}+b_z^j{K}_{zyy}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yy}^{l\to i}},i=1,2,\cdots ,N\\ {\displaystyle {\sum }_{j=1}^N{b}_x^j{K}_{xyz}^{j\to i}+b_y^j{K}_{yyz}^{j\to i}+b_z^j{K}_{zyz}^{j\to i}}=-{\displaystyle {\sum }_l^{\gamma }{\sigma }_{yz}^{l\to i}},i=1,2,\cdots ,N\end{array}$ (16)

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:

$\left[\begin{array}{ccc}\begin{array}{ccc}{K}_{xyx}^{1\to 1}& K_{yyx}^{1\to 1}& K_{zyx}^{1\to 1}\\ K_{xyy}^{1\to 1}& K_{yyy}^{1\to 1}& K_{zyy}^{1\to 1}\\ K_{xyz}^{1\to 1}& K_{yyz}^{1\to 1}& K_{zyz}^{1\to 1}\end{array}& \cdots & \begin{array}{ccc}{K}_{xyx}^{N\to 1}& K_{yyx}^{N\to 1}& K_{zyx}^{N\to 1}\\ K_{xyy}^{N\to 1}& K_{yyy}^{N\to 1}& K_{zyy}^{N\to 1}\\ K_{xyz}^{N\to 1}& K_{yyz}^{N\to 1}& K_{zyz}^{N\to 1}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{ccc}{K}_{xyx}^{1\to N}& K_{yyx}^{1\to N}& K_{zyx}^{1\to N}\\ K_{xyy}^{1\to N}& K_{yyy}^{1\to N}& K_{zyy}^{1\to N}\\ K_{xyz}^{1\to N}& K_{yyz}^{1\to N}& K_{zyz}^{1\to N}\end{array}& \cdots & \begin{array}{ccc}{K}_{xyx}^{N\to N}& K_{yyx}^{N\to N}& K_{zyx}^{N\to N}\\ K_{xyy}^{N\to N}& K_{yyy}^{N\to N}& K_{zyy}^{N\to N}\\ K_{xyz}^{N\to N}& K_{yyz}^{N\to N}& K_{zyz}^{N\to N}\end{array}\end{array}\right]\left\{\begin{array}{c}{b}_x^1\\ b_y^1\\ b_z^1\\ ⋮\\ b_x^N\\ b_y^N\\ b_z^N\end{array}\right\}=-\left\{\begin{array}{c}{\displaystyle \sum {\sigma }_{yx}^{l\to 1}}\\ {\displaystyle \sum {\sigma }_{yy}^{l\to 1}}\\ {\displaystyle \sum {\sigma }_{yz}^{l\to 1}}\\ ⋮\\ {\displaystyle \sum {\sigma }_{yx}^{l\to N}}\\ {\displaystyle \sum {\sigma }_{yy}^{l\to N}}\\ {\displaystyle \sum {\sigma }_{yz}^{l\to N}}\end{array}\right\}$ (17)

The matrix containing the kernel terms is labelled here the kernel matrix. Its size is 3N × 3N (=9N²). 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 ${\sigma }_{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.

3. Numerical Considerations

3.1. Solving the System of Equations

Equation (17) can be expressed in matrix form as,

$\left[K\right]B=F$ (18)

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. Notice that [K] and its elements couple B and F. The elements of [K] come from the surface loops. The elements of F come from the crystal dislocations. The surface loop density determines the dimension of [K]. The elements of [K] are computed from the surface dislocation loops with no contribution from any crystal dislocations inside the computational domain/cell. For F, however, its elements are computed (just like [K] at the collocation points), from the effect of the real dislocations in the computational domain/cell. This is important to observe since if more dislocations are present inside the crystal, this will affect only the calculation of F and won’t affect the matrix [K]. This is very helpful in dislocation dynamics simulations as the kernel terms only need to be evaluated at the beginning of the simulation, whereas the elements of F need to be updated as time evolves, i.e. with every time stepping. Hence, as the dislocations glide in the crystal and change their configuration, so do the values for the B vector elements. Consider Figure 2 which shows a computational box with six free surfaces. On each of these surfaces a zero stress-traction condition ( $\stackrel{\to }{T}=\stackrel{\to }{0}$ ) 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, $\left|k_{ii}\right|\ge {\displaystyle {\sum }_{j\ne i}{k}_{ij}}$. Iterative methods, like the Gauss-Seidel method [29], can be employed for cases when the matrix is diagonally-dominant. Also, successive relaxation can be utilized to solve the system of linear Equations (see Equations (18)). For faster computing, we have chosen a relaxation parameter λ such that $0<\lambda <2$. If $b_f$ is a solution of Equation (18) after iteration f then the new value of $b_f$ is modified by a weighted average of the results from the previous and present iterations as:

$b_f^{new}=\lambda b_f^{new}+\left(1-\lambda \right)b_f^{old}$ (19)

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.

3.2. Solution Time

Figure 4 shows a time comparison between two different shapes of surface loops (see Figure 2 and Figure 3) used to compute the correction term. The normalized time is recorded while computing the correction term using 2500 surface loops of both shapes. The figure shows about 30% time improvement for using a triangular surface mesh. This will make a remarkable time improvement in dynamic simulations by saving a significant amount of time over many time steps.

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 20,000b. The magnitude of the Burgers vector $\stackrel{\to }{b}$ is b. The depth of the linear segment below the surface is 1000b. The length of the linear segment is 100b. For points A and B, the following exact coordinates were chosen $\left(100b,0,0\right)$ and $\left(0,0,0\right)$, respectively. Also, $\stackrel{\to }{b}$ for segment AB is chosen for simplicity as

$\left(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\right)$. For the material, take Poisson ratio ν as 0.383. Consider a line,

for field point calculations, beneath the surface and parallel to it at a depth of 400b (along $y=0$ or the x-axis) (See Figure 6 and Figure 7).

Figure 4. Computational time of any field stresses for two surface loop geometries used in computing the correction term.

Figure 5. A Horizontal subsurface segment.

(a) (b) (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. (a) # dislocation loops 400; (b) # dislocation loops 900; (c) # dislocation loops 1600; (d) # dislocation loops 2500.

(a) (b) (c) (d)

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. (a) # dislocation loops 400; (b) # dislocation loops 900; (c) # dislocation loops 1600; (d) # dislocation loops 2500.

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 Maurissen and Capella (1974) [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 $L=100b$. For points A and B, the following exact coordinates were chosen $\left(0,0,-1,100b\right)$ and $\left(0,0,-1,000b\right)$, respectively. Also, $\stackrel{\to }{b}$ for segment AB is

chosen as $\left(\frac{1}{\sqrt{3}},0,\frac{1}{\sqrt{3}}\right)$. 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 $x=0$ 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 compare the analytical solution by [10],

Figure 8. A Vertical subsurface segment.

(a) (b) (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.

(a) (b) (c) (d)

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 three- dimensional 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³. Also,

Poisson ratio $\nu =0.3$, along with a dislocation mobility equaling $1000\frac{1}{\text{Pa}\cdot \text{s}}$. A

dislocation source (of the Frank-Read type) is placed along the x-axis. The coordinates for its end points are $\left(-2000b,0,0\right)$ and $\left(2000b,0,0\right)$. The Burgers vector was $\left(0,1,0\right)$. Lastly, the applied strain-rate amount was 10 s⁻¹. For Figure 11 and Figure 12, the stress and strain are the yz-components.

Figure 11. The effect of varying the number of surface dislocation loops on the stress- strain curve (using constant strain-rate loading). The correction term is calculated from rectangular surface loops. (a) # dislocation loops 64; (b) # dislocation loops 144; (c) # dislocation loops 256; (d) Effect of surface loop density.

Figure 12. The effect of varying the number of surface dislocation loops on the stress-strain curve (using constant strain-rate loading). The correction term is calculated from triangular surface loops. (a) # dislocation loops 64; (b) # dislocation loops 144; (c) # dislocation loops 256; (d) Effect of density of the surface loops.

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 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₂ observed without surface effects is about 14 percent more than the bowing dimension L₁ 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(a) shows the z-component of the Peach-Koehler force acting at the center of a linear dislocation segment varying the dislocation depth from the free surface. The parameters for the calculations here are the same as those used earlier for a horizontal sub-surface segment with the exception of varying the depth (Figure 14(a)) or the length (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.

(a) (b)

Figure 13. Simulation snapshots at the critical bowing of a Frank-Read source of length L. (a) Illustration of the critical bowing of simulated Frank-Read sources. The view is perpendicular to the slip plane; (b) Markers showing the critical bowing stress and strain values, i.e. start of plasticity.

(a) (b)

Figure 14. The z-component of the Peach-Koehler force on a dislocation segment as the z-depth and length are varied. (a) Varying the dislocation depth; (b) Varying the dislocation length.

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.

4. 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 authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hull, D. and Bacon, D.J. (2011) Introduction to Dislocations. Butterworth-Heinemann, Oxford. https://doi.org/10.1016/B978-0-08-096672-4.00002-5
[2]	Hirth, J.P. and Lothe, J. (1982) Theory of Dislocations. John Wiley & Sons, Hoboken.
[3]	Eshelby, J.D. and Stroh, A.N. (1951) CXL. Dislocations in Thin Plates. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42, 1401-1405. https://doi.org/10.1080/14786445108560958
[4]	Jeans, J. (1925) Mathematical Theory of Electricity and Magnetism. CUP Archive.
[5]	Yoffe, E.H. (1961) A Dislocation at a Free Surface. Philosophical Magazine, 6, 1147-1155. https://doi.org/10.1080/14786436108239675
[6]	Bastecká, J. (1964) Interaction of Dislocation Loop with Free Surface. Czechoslovak Journal of Physics, 14, 430-442. https://doi.org/10.1007/BF01689476
[7]	Groves, P.P. and Bacon, D.J. (1970) The Dislocation in a Semi-Infinite Isotropic Medium. Fundamental Aspects of Dislocation Theory, 1, 35-45.
[8]	Groves, P.P. and Bacon, D.J. (1970) The Dislocation Loop near a Free Surface. Philosophical Magazine, 22, 83-91. https://doi.org/10.1080/14786437008228153
[9]	Maurissen, Y. and Capella, L. (1974) Stress Field of a Dislocation Segment Parallel to a Free Surface. Philosophical Magazine, 29, 1227-1229. https://doi.org/10.1080/14786437408226608
[10]	Maurissen, Y. and Capella, L. (1974) Stress Field of a Dislocation Segment Perpendicular to a Free Surface. Philosophical Magazine, 30, 679-683. https://doi.org/10.1080/14786439808206591
[11]	Comninou, M. and Dundurs, J. (1975) The Angular Dislocation in a Half Space. Journal of Elasticity, 5, 203-216. https://doi.org/10.1007/BF00126985
[12]	Lothe, J., Indenbom, V.L. and Chamrov, V.A. (1982) Elastic Field and Self-Force of Dislocations Emerging at the Free Surfaces of an Anisotropic Halfspace. Physica Status Solidi (B), 111, 671-677. https://doi.org/10.1002/pssb.2221110231
[13]	Gosling, T.J. and Willis, J.R. (1994) A Line-Integral Representation for the Stresses Due to an Arbitrary Dislocation in an Isotropic Half-Space. Journal of the Mechanics and Physics of Solids, 42, 1199-1221. https://doi.org/10.1016/0022-5096(94)90032-9
[14]	Devincre, B. (1995) Three Dimensional Stress Field Expressions for Straight Dislocation Segments. Solid State Communictions, 93, 875-878. https://doi.org/10.1016/0038-1098(94)00894-9
[15]	Kubin, L., Canova, G., Condat, M., Devincre, B., Pontikis, V. and Brechet, Y. (1992) Dislocation Microstructures and Plastic Flow: A 3D Simulation. Solid State Phenomena, 23-24, 455-472. https://doi.org/10.4028/www.scientific.net/SSP.23-24.455
[16]	Canova, G.R. and Fivel, M.C. (1999) Developing Rigorous Boundary Conditions to Simulations of Discrete Dislocation Dynamics. Modelling and Simulation in Materials Science and Engineering, 7, 753-768. https://doi.org/10.1088/0965-0393/7/5/308
[17]	Hartmaier, A., Fivel, M.C., Canova, G.R. and Gumbsch, P. (1999) Image Stresses in a Free-Standing Thin Film. Modelling and Simulation in Materials Science and Engineering, 7, 781-793. https://doi.org/10.1088/0965-0393/7/5/310
[18]	El-Azab, A. (2000) The Boundary Value Problem of Dislocation Dynamics. Modelling and Simulation in Materials Science and Engineering, 8, 37-54. https://doi.org/10.1088/0965-0393/8/1/304
[19]	Deng, J., El-Azab, A. and Larson, B.C. (2008) On the Elastic Boundary Value Problem of Dislocations in Bounded Crystals. Philosophical Magazine, 88, 3527-3548. https://doi.org/10.1080/14786430802558544
[20]	Zhou, C., Biner, S.B. and LeSar, R. (2010) Discrete Dislocation Dynamics Simulations of Plasticity at Small Scales. Acta Materialia, 58, 1565-1577. https://doi.org/10.1016/j.actamat.2009.11.001
[21]	Crone, J.C., Munday, L.B. and Knap, J. (2015) Capturing the Effects of Free Surfaces on Void Strengthening with Dislocation Dynamics. Acta Materialia, 101, 40-47. https://doi.org/10.1016/j.actamat.2015.08.067
[22]	Jing, P., Khraishi, T., Young, J.A. and Wirth, B.D. (2005) Multi-Scale Simulations of the Effects of Irradiation-induced Voids and Helium Bubbles on the Mechanical Properties of Aluminum. Philosophical Magazine, 85, 757-767. https://doi.org/10.1080/14786430412331319958
[23]	Jing, P., Khraishi, T., Zepeda-Ruiz, L.A. and Wirth, B.D. (2009) The Elastic Fields of Sub-Surface Dislocation Loops: A Comparison between Analytical Continuum-Theory Solutions and Atomistic Calculations. International Journal of Theoretical and Applied Multiscale Mechanics, 1, 71-85. https://doi.org/10.1504/IJTAMM.2009.022472
[24]	Khraishi, T.A. and Zbib, H.M. (2001) Dislocation Dynamics Simulations of the Interaction between a Short Rigid Fiber and a Glide Circular Dislocation Pile-Up. Computational Materials Science, 24, 310-322. https://doi.org/10.1016/S0927-0256(01)00253-1
[25]	Khraishi, T.A. and Zbib, H.M. (2002) Free-Surface Effects in 3D Dislocation Dynamics: Formulation and Modeling. Journal of Engineering Materials and Technology, 124, 342-351. https://doi.org/10.1115/1.1479694
[26]	Yan, L., Khraishi, T.A., Shen, Y.-L. and Horstemeyer, M.F. (2004) A Distributed-Dislocation Method for Treating Free-Surface Image Stresses in Three-Dim-ensional Dislocation Dynamics Simulations. Modelling and Simulation in Materials Science and Engineering, 12, S289-S301. https://doi.org/10.1088/0965-0393/12/4/S01
[27]	Siddique, A.B. and Khraishi, T. (2020) Numerical Methodology for Treating Static and Dynamic Dislocation Problems near a Free Surface. Journal of Physics Communications, 4, Article ID: 055005. https://doi.org/10.1088/2399-6528/ab8ff9
[28]	Siddique, A.B. and Khraishi, T.A. (2020) A Dislocation Near a Cylindrical Hole: A Numerical Treatment. Proceedings of the 2020 ASEE Gulf-Southwest Annual Conference, 23-29 April 2020, 3. https://peer.asee.org/35939
[29]	Chapra, S.C. and Canale, R.P. (1998) Numerical Methods for Engineers. McGraw-Hill, New York.
[30]	Khraishi, T.A., Yan, L. and Shen, Y.-L. (2004) Dynamic Simulations of the Interaction between Dislocations and Dilute Particle Concentrations in Metal-Matrix Composites (MMCs). International Journal of Plasticity, 20, 1039-1057. https://doi.org/10.1016/j.ijplas.2003.10.003
[31]	Zbib, H.M., Rhee, M. and Hirth, J. (1998) On Plastic Deformation and the Dynamics of 3D Dislocations. International Journal of Mechanical Sciences, 40, 113-127. https://doi.org/10.1016/S0020-7403(97)00043-X
[32]	Meyers, M.A. and Chawla, K.K. (1998) Mechanical Behavior of Materials. Cambridge University Press, Cambridge.