Nonsingular Stress Distribution of Edge Dislocations near Zero-Traction Boundary

Abstract

Among many types of defects present in crystalline materials, dislocations are the most influential in determining the deformation process and various physical properties of the materials. However, the mathematical description of the elastic field generated around dislocations is challenging because of various theoretical difficulties, such as physically irrelevant singularities near the dislocation-core and nontrivial modulation in the spatial distribution near the material interface. As a theoretical solution to this problem, in the present study, we develop an explicit formulation for the nonsingular stress field generated by an edge dislocation near the zero-traction surface of an elastic medium. The obtained stress field is free from nonphysical divergence near the dislocation-core, as compared to classical solutions. Because of the nonsingular property, our results allow the accurate estimation of the effect of the zero-traction surface on the near-surface stress distribution, as well as its dependence on the orientation of the Burgers vector. Finally, the degree of surface-induced modulation in the stress field is evaluated using the concept of the $L_2$-norm for function spaces and the comparison with the stress field in an infinitely large system without any surface.

1. Introduction

Continuum mechanics theory enables the calculation of the spatial distribution of stress (and strain) inside a material and its deformation under loading. If the interior of the material is ideally uniform and free from defects, then the elastic field is smooth and continuous over the system. However, in actual materials, a certain number of defects are often embedded in a discrete and inhomogeneous manner. Thus, contributions from those defects to the mechanical properties of the system in the realm of continuum theory should be considered [1,2].

Of many types of defects, dislocations are ubiquitous in crystalline metals and alloys and are typically found at a high density of ${10}^8$ or more per 1 cm². Each dislocation is a line defect that disturbs the lattice arrangement of the perfect crystal, locally breaking the translational symmetry and acting as the source of stress and strain that are distributed in the materials [3]. The broken symmetry near the dislocation-core and its movement due to external forces or thermal excitations play primary roles in the physical phenomena. For instance, the plastic deformation of a crystal is caused by a large flow of dislocations [4,5], and work hardening is caused by the dislocation arrangement formed in the crystal [6]. Otherwise, the nucleation and mobility of dislocations near the tip of a crack have a strong impact on the fracture mechanics of the material [7]. Furthermore, an in-depth understanding of dislocation dynamics and dislocation-induced elastic fields near traction-free boundaries [8,9,10] and grain boundaries [9,11,12,13,14], as well as the assessment of the boundary conditions used in the existing formulations [15] are critical to predicting the performance of advanced materials: they include gold nanowires [16,17], barium titanium nanomaterials [18], functionally graded materials [19], and other metallic and semiconductor materials [20,21,22,23].

The mechanics of dislocations and their long-range correlated behaviors [24] have been addressed successfully using conventional (traditional) dislocation theory, as well as coarse-grained approaches [25], at least from a macroscopic viewpoint [26,27,28]. However, when considering nanometric-scale mechanics, conventional dislocation theory may not work adequately. The most significant drawback of the conventional theory is that stress and strain fields diverge infinitely near the dislocation-core. Because of this divergence, the theory does not fully describe the elastic field near the dislocation-core. This loss of accuracy hinders a detailed analysis of the near-core elastic field. Still, this analysis is crucial in the cases where a large number of dislocations are densely packed in a nanoscopic region or close to the interface (free surfaces, grain boundaries, etc.), as experimentally observed in the middle and late stages of the plastic deformation of metals under cyclic loading [4,29]. Furthermore, the conventional dislocation theory is intrinsically size-independent, and thus, no characteristic length appears. This complicates the theoretical interpretation of the characteristic length of the self-organized structures of dislocation ensembles that occur spontaneously during the cyclic deformation of metals [30,31].

To overcome these difficulties, various nonclassical continuum theories of dislocations have been developed in recent decades. These include theories of nonlocal elasticity [32,33], gradient elasticity [34,35], the gauge-field approach [36,37], and other sophisticated numerical approaches [38,39]. One notable feature of these nonstandard approaches is that stress singularities at the dislocation-core are eliminated in the obtained solutions so that the stress distribution becomes smooth and continuous without a divergence [33,40]. In addition, all these approaches engender characteristic length scales defined by theory-dependent parameters (e.g., nonlocality parameter, gradient coefficient, or related material parameters), which correspond to the spatial dimension of dislocation-cores. These nonsingular solutions turned out to be consistent with the numerical simulations of elastic fields realized in actual materials [41,42] and those caused by uniformly moving dislocations in an infinite body [43]. The practical significance of the nonsingular solutions has also been confirmed in the study of stress fields near crack tips in isotropic elastic materials [44]. Yet, most theoretical studies reported on the nonsingular problem so far have considered the elastic field of dislocations existing in an infinitely large system, and only a few studies have handled a finite or semi-infinite system endowed with a surface boundary. For microsized and nanosized materials [4,30], the ratio of the surface area to sample volume is so high that the effect of a surface boundary should be more prominent than that observed in macroscale samples.

Against this backdrop, we were motivated to derive a theoretical representation of the nonsingular stress field generated by an edge dislocation near the free surface of a semi-infinitely large elastic medium. The Burgers vector of the edge dislocation was assumed to be oblique to the free surface, considering that, in the actual plastic deformation process, the slip plane of edge dislocations in metals is often oriented diagonally to the load direction. Our analytical expression allows us to quantify the spatial modulation of the nonsingular stress field of an edge dislocation in the presence of a zero-traction surface boundary. The stress-field formula derived in this study provides a theoretical basis for analyzing the behavior of dislocations accumulated near the free surface and the interaction between adjacent dislocations on micrometric and nanometric scales [4,30].

2. Nonsingular Elastic Field in an Infinite Medium

We assume an infinitely long, straight edge dislocation that extends along the z-axis in the Cartesian coordinate system. The dislocation-core is located at the origin of the x-y coordinate plane. For convenience, we first consider an infinite elastic system without boundaries, and then, a semi-infinite system, including a planar free surface with no traction. Because the system spreads infinitely in the z direction, the edge dislocation produces a two-dimensional state of plane stress, as well as plane strain defined by $u_z=0$ and $d{u}_i/\mathit{dz}=0$ for $i=x,y,z$, where $u_i$ is the i-directional displacement of an infinitesimal element of the system.

It is known that any in-plane stress problem can be reduced to exploring an appropriate biharmonic function [45], $\varphi (x,y)$, the so-called Airy’s stress function, that satisfies ${\nabla }^4\varphi =0$. Once the solution of $\varphi (x,y)$ in a particular domain of interest is obtained under the given boundary conditions, the stress components within the domain can be derived through partial differentiation of $\varphi (x,y)$ as

$${\sigma }_{\mathit{xx}}(x,y)=\frac{{\partial }^2\varphi }{\partial y^2},{\sigma }_{\mathit{yy}}(x,y)=\frac{{\partial }^2\varphi }{\partial x^2},{\tau }_{\mathit{xy}}(x,y)=\frac{-{\partial }^2\varphi }{\partial x\partial y}.$$

(1)

The existing work based on the gradient elasticity theory [34,35] and the gauge field theory [36,37] has unveiled that Airy’s stress function $\varphi (x,y)$ associated with a nonsingular stress field of an edge dislocation with the Burgers vector parallel to the x-axis obeys the following inhomogeneous Helmholtz equation:

$$({\nabla }^2-{\kappa }^2)\varphi (x,y)={\kappa }^2\mathcal{G}ylog\left(\sqrt{x^2+y^2}\right).$$

(2)

Here, $\mathcal{G}$ is a material-dependent constant defined by

$$\mathcal{G}=\frac{\mu b}{2\pi (1-{\nu }^{*})},$$

(3)

with the elastic shear modulus $\mu$, the magnitude of Burgers vector b, and Poisson’s ratio ${\nu }^{*}$. The solution of Equation (2) is represented by (see Appendix A)

$$\varphi (x,y)=-2\mathcal{G}y\left[\frac{log\left(\sqrt{x^2+y^2}\right)}{2}+\frac{1-\kappa \sqrt{x^2+y^2}·K_1\left(\kappa \sqrt{x^2+y^2}\right)}{{\kappa }^2\left(x^2+y^2\right)}\right].$$

(4)

In Equation (4), $K_n$ indicates the nth-order modified Bessel function of the second kind (see Appendix B). By taking the partial derivative of $\varphi (x,y)$ twice according to Equation (1), we obtain the nonsingular solution of the stress field around the edge dislocation located at the origin of an infinitely large elastic medium.

The factor $\kappa$ in Equation (2), having dimension of inverse length, is decisive for the elastic field obtained from $\varphi (x,y)$ to be free of singularities at the dislocation-core (i.e., the origin of the x–y coordinate plane). In the framework of the gauge field theory, $\kappa$ is defined by the coefficient of the constitutive relation between translational gauge-invariant physical state quantities; $\kappa$ determines the position and magnitude of the extrema of the nonsingular stress field. In contrast, in the gradient elasticity theory, ${\kappa }^{-2}$ serves as the gradient coefficient appearing in the modified constitutive relation. It should be emphasized that, in both theories, ${\kappa }^{-1}$ is regarded as a characteristic length scale of the system. In the limit of ${\kappa }^{-1}\to 0$ (or equivalently, $\kappa \to \infty$), the problem under consideration falls into a conventional one suffering from the dislocation-core singularity.

It is noteworthy that the nonsingular solution of $\varphi (x,y)$ given by Equation (4) covers the classical counterpart suffered from a dislocation-core singularity, as confirmed by considering the limiting behaviors of the modified Bessel function $K_n\left(u\right)$. In general, $K_n\left(u\right)$ is a monotonic decreasing function with u. In particular, when $n=1$, its asymptotic decays for $u\gg 1$ and $u\ll 1$ are approximated by [46]

$$K_1\left(u\right)\simeq \sqrt{\frac{\pi }{2u}}{e}^{-u}(u\gg 1)\mathrm{and}{K}_1\left(u\right)\simeq \frac{1}{u}(u\ll 1).$$

(5)

At $\sqrt{x^2+y^2}\gg {\kappa }^{-1}$, therefore, the second term in the square brackets in Equation (4) vanishes so that $\varphi (x,y)$ reduces to the classical stress function for an edge dislocation:

$$\varphi (x,y)=-\mathcal{G}ylog\left(\sqrt{x^2+y^2}\right).$$

3. Semi-Infinite Problem Solution

3.1. Outline of Strategy

The stress field in a semi-infinite plane with a free surface can be derived by considering the superposition of three Airy’s stress functions such that their sum satisfies the traction-free boundary condition. To proceed with this argument, let us suppose that a semi-infinitely large elastic medium spans in the region of $x<0$, and an edge dislocation is positioned at $(x,y)=(-d,0)$ with $d>0$. If there were no free surfaces and, thus, the elastic medium was infinitely large, the stress field caused by this edge dislocation could have been derived from Airy’s stress function ${\varphi }^{\left(\mathrm{re}\right)}(x,y)$ given by

$${\varphi }^{\left(\mathrm{re}\right)}(x,y)=\varphi (x+d,y),$$

(6)

where the superscript “re” indicates that the edge dislocation at $x=-d$ actually exists in the real elastic material under consideration. In Equation (6), $\varphi (x+d,y)$ is a function obtained by rewriting the argument x of $\varphi (x,y)$ expressed by Equation (4) to $x+d$.

As a matter of course, the stress field derived from ${\varphi }^{\left(\mathrm{re}\right)}(x,y)$ does not satisfy the traction-free boundary conditions expressed by ${\sigma }_{\mathit{xx}}=0$ and ${\tau }_{\mathit{xy}}=0$ along $x=0$. Of these two, the former requirement is accomplished by superposing another Airy’s stress function defined by

$${\varphi }^{\left(\mathrm{im}\right)}(x,y)=\varphi (x-d,y),$$

(7)

onto ${\varphi }^{\left(\mathrm{re}\right)}(x,y)$. The superscript “im” indicates a contribution from a virtual image dislocation, which is assumed to be located at $(x,y)=(+d,0)$ (outside the free surface), but it is not present in the real material. We also assume that the Burgers vector of the image dislocation has the same magnitude as that of the real dislocation, whereas their orientations can be different from each other. By tuning the relative orientations of the two Burgers vectors, we can make the stress field obtained from the sum of ${\varphi }^{\left(\mathrm{re}\right)}(x,y)+{\varphi }^{\left(\mathrm{im}\right)}(x,y)$ satisfy the condition of ${\sigma }_{\mathit{xx}}(x,y)=0$ at $x=0$. Here, the contribution ${\sigma }_{\mathit{xx}}^{\left(\mathrm{im}\right)}(x,y)$ from the image dislocation cancels out the contribution ${\sigma }_{\mathit{xx}}^{\left(\mathrm{re}\right)}(x,y)$ from the real dislocation at every point along $x=0$. This technique of canceling the stress component ${\sigma }_{\mathit{xx}}$ at a given boundary is called the image force method [47,48], which is a mechanical analog of the image charge construction employed in electromagnetism [49].

The remaining requirement of ${\tau }_{\mathit{xy}}=0$ along $x=0$ is met by introducing a third “excess” Airy’s stress function, designated by ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)$, such that the two components of the stress field derived from it, ${\sigma }_{xx}^{\left(\mathrm{ex}\right)}$ and ${\tau }_{\mathit{xx}}^{\left(\mathrm{ex}\right)}$, satisfy the following conditions at $x=0$:

$${\sigma }_{\mathit{xx}}^{\left(\mathrm{ex}\right)}(0,y)=0\mathrm{and}{\tau }_{\mathit{xy}}^{\left(\mathrm{ex}\right)}(0,y)=-\left\{{\tau }_{\mathit{xy}}^{\left(\mathrm{re}\right)}(0,y)+{\tau }_{\mathit{xy}}^{\left(\mathrm{im}\right)}(0,y)\right\}.$$

(8)

The latter in Equation (8) indicates that the contribution from the excess Airy’s function cancels out those from the two dislocations (real and image). If such ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)$ is found, the sum of three Airy’s functions given by

$${\varphi }^{\left(\mathrm{all}\right)}(x,y)={\varphi }^{\left(\mathrm{re}\right)}(x,y)+{\varphi }^{\left(\mathrm{im}\right)}(x,y)+{\varphi }^{\left(\mathrm{ex}\right)}(x,y)$$

(9)

satisfies both ${\sigma }_{\mathit{xx}}^{\left(\mathrm{all}\right)}(x,y)=0$ and ${\tau }_{\mathit{xy}}^{\left(\mathrm{all}\right)}=0$ along $x=0$. Thus, it can be used to obtain the three stress components distributed over the semi-infinite plane at $x<0$.

3.2. Stress Due to the Real Dislocation at $x=-d$

In metals, edge dislocations that move in a slanted direction with respect to the free surface play a vital role in plastic deformation [4,6,30]. Therefore, we developed an explicit and tractable formula for the nonsingular stress field of a semi-infinite system to infer the surface effect on the fundamental properties of the edge dislocations near the free surface.

To simplify the mathematical expressions, we henceforth use the following notations:

$$x_{\mathrm{p}}\equiv x+d,c_{\mathrm{p}}(x,y)\equiv \kappa \sqrt{{(x+d)}^2+y^2},$$

(10)

where the subscript “$\mathrm{p}$” indicates the plus sign “+” contained in the term $x+d$. Next, we calculate the three components of the stress, ${\sigma }_{\mathit{xx}}^{\left(\mathrm{re}\right)}(x,y)$, ${\sigma }_{\mathit{yy}}^{\left(\mathrm{re}\right)}(x,y)$, and ${\tau }_{\mathit{xy}}^{\left(\mathrm{re}\right)}(x,y)$, generated by the edge dislocation located at $(x,y)=(-d,0)$. Herein, the Burgers vector is oriented in the slanted direction by $\alpha$ with respect to the x-axis. The coordinate rotation of the stress components is obtained by partial differentiation of ${\varphi }^{\left(\mathrm{re}\right)}(x,y)$ as

$${\sigma }_{\mathit{xx}}^{\left(\mathrm{re}\right)}(x,y)=\mathcal{G}\frac{4{\chi }_{\mathit{xx}}^{\left[1\right]}+c_{\mathrm{p}}^2{\chi }_{\mathit{xx}}^{\left[2\right]}+2c_{\mathrm{p}}^2{K}_0\left(c_{\mathrm{p}}\right){\chi }_{\mathit{xx}}^{\left[3\right]}+2c_{\mathrm{p}}{K}_1\left(c_{\mathrm{p}}\right){\chi }_{\mathit{xx}}^{\left[4\right]}}{c_{\mathrm{p}}^2{\left(x_{\mathrm{p}}^2+y^2\right)}^2},$$

(11)

with

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xx}}^{\left[1\right]}& =& x_{\mathrm{p}}\left(x_{\mathrm{p}}^2-3y^2\right)sin\alpha +y\left(3x_{\mathrm{p}}^2-y^2\right)cos\alpha ,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xx}}^{\left[2\right]}& =& -x_{\mathrm{p}}\left(x_{\mathrm{p}}^2-y^2\right)sin\alpha -y\left(3x_{\mathrm{p}}^2+y^2\right)cos\alpha ,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xx}}^{\left[3\right]}& =& -{\chi }_{xx}^{\left[1\right]},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xx}}^{\left[4\right]}& =& -2{\chi }_{xx}^{\left[1\right]}+c_{\mathrm{p}}^2{y}^2\left(x_{\mathrm{p}}sin\alpha +ycos\alpha \right),\hfill \end{array}$$

(15)

and

$${\sigma }_{yy}^{\left(\mathrm{re}\right)}(x,y)=\mathcal{G}\frac{4{\chi }_{yy}^{\left[1\right]}+c_{\mathrm{p}}^2{\chi }_{yy}^{\left[2\right]}+2c_{\mathrm{p}}^2{K}_0\left(c_{\mathrm{p}}\right){\chi }_{yy}^{\left[3\right]}+2c_{\mathrm{p}}{K}_1\left(c_{\mathrm{p}}\right){\chi }_{yy}^{\left[4\right]}}{c_{\mathrm{p}}^2{\left(x_{\mathrm{p}}^2+y^2\right)}^2},$$

(16)

with

$$\begin{array}{ccc}\hfill {\chi }_{yy}^{\left[1\right]}& =& -x_{\mathrm{p}}\left(x_{\mathrm{p}}^2-3y^2\right)sin\alpha -y\left(3x_{\mathrm{p}}^2-y^2\right)cos\alpha \left(=-{\chi }_{xx}^{\left[1\right]}\right),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\chi }_{yy}^{\left[2\right]}& =& -x_{\mathrm{p}}\left(x_{\mathrm{p}}^2+3y^2\right)sin\alpha +y\left(x_{\mathrm{p}}^2-y^2\right)cos\alpha ,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\chi }_{yy}^{\left[3\right]}& =& -{\chi }_{yy}^{\left[1\right]},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\chi }_{yy}^{\left[4\right]}& =& -2{\chi }_{yy}^{\left[1\right]}+c_{\mathrm{p}}^2{x}_{\mathrm{p}}^2\left(x_{\mathrm{p}}sin\alpha +ycos\alpha \right),\hfill \end{array}$$

(20)

and

$${\tau }_{xy}^{\left(\mathrm{re}\right)}(x,y)=\mathcal{G}\frac{4{\chi }_{\mathit{xy}}^{\left[1\right]}+c_{\mathrm{p}}^2{\chi }_{\mathit{xy}}^{\left[2\right]}+2c_{\mathrm{p}}^2{K}_0\left(c_{\mathrm{p}}\right){\chi }_{\mathit{xy}}^{\left[3\right]}+2c_{\mathrm{p}}{K}_1\left(c_{\mathrm{p}}\right){\chi }_{\mathit{xy}}^{\left[4\right]}}{c_{\mathrm{p}}^2{\left(x_{\mathrm{p}}^2+y^2\right)}^2},$$

(21)

with

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xy}}^{\left[1\right]}& =& y\left(3x_{\mathrm{p}}^2-y^2\right)sin\alpha -x_{\mathrm{p}}\left(x_{\mathrm{p}}^2-3y^2\right)cos\alpha ,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xy}}^{\left[2\right]}& =& \left(x_{\mathrm{p}}^2-y^2\right)\left(-ysin\alpha +x_{\mathrm{p}}cos\alpha \right),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xy}}^{\left[3\right]}& =& -{\chi }_{xy}^{\left[1\right]},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\chi }_{\mathit{xy}}^{\left[4\right]}& =& -2{\chi }_{xy}^{\left[1\right]}+c_{\mathrm{p}}^2\left(-x_{\mathrm{p}}y\right)\left(x_{\mathrm{p}}sin\alpha +ycos\alpha \right).\hfill \end{array}$$

(25)

Note that the three expressions in Equations (11), (16) and (21), have the same functional form, whereas the definitions of ${\chi }_{ij}^{\left[\mathrm{k}\right]}$$(i,j=x\mathrm{or}y,1\le k\le 4)$ are different depending on the components. It also should be noted that, at sufficiently far distances from the dislocation-core (i.e., at ${(x-d)}^2+y^2\gg {\kappa }^{-2}$), only the term proportional to ${\chi }_{ij}^{\left[2\right]}$ remains non-negligible. Hence, the three components are reduced to the classical counterpart. These two noteworthy facts hold true for the stress components associated with the image dislocation, which will be described in the next subsection.

3.3. Stress Due to the Image Dislocation at $x=+d$

The three components of the stress generated by the image edge dislocation at $(x,y)=(+d,0)$, denoted by ${\sigma }_{xx}^{\left(\mathrm{im}\right)}(x,y)$, ${\sigma }_{yy}^{\left(\mathrm{im}\right)}(x,y)$, and ${\tau }_{xy}^{\left(\mathrm{im}\right)}(x,y)$, can be obtained by the variable transformation for all the mathematical expressions regarding the real dislocation counterparts given in the previous subsection. Once they are rewritten as $\mathcal{G}\to -\mathcal{G}$, $\alpha \to \beta$, and $x+d\to x-d$, we obtain the expressions of the stress component of the image dislocation drawn in Figure 1.

By setting the slanted angles to $\alpha =\theta$ and $\beta =-\theta$, it is easily proven that, at $x=0$, the obtained expressions of ${\sigma }_{xx}^{\left(\mathrm{re}\right)}$ and ${\sigma }_{xx}^{\left(\mathrm{im}\right)}$ cancel each other out. This cancellation means that the virtually introduced image dislocation makes it possible to eliminate the ${\sigma }_{xx}$ component at the surface, on which no traction force must be present.

3.4. Sum of the Shear Stress Components at a Free Surface

For later use in the derivation of the excess Airy function ${\varphi }^{\left(\mathrm{ex}\right)}$, we now compute the sum of the two shear components of the two above-mentioned dislocations along the line of $x=0$, denoted by ${\tau }_{xy}^{\left(\mathrm{re}\right)}(0,y)+{\tau }_{xy}^{\left(\mathrm{im}\right)}(0,y)$. Using the notation of

$$c_0\left(y\right)=\kappa \sqrt{y^2+d^2},$$

(26)

the sum of the shear components is rewritten as

$$\begin{array}{ccc}\hfill {\tau }_{xy}^{\left(\mathrm{re}\right)}(0,y)+{\tau }_{xy}^{\left(\mathrm{im}\right)}(0,y)& =& 2\mathcal{G}\frac{4{\omega }^{\left[1\right]}+c_0^2{\omega }^{\left[2\right]}+2c_0^2{K}_0\left(c_0\right){\omega }^{\left[3\right]}+2c_0{K}_1\left(c_0\right){\omega }^{\left[4\right]}}{c_0^2{\left(y^2+d^2\right)}^2},\hfill \end{array}$$

(27)

with

$$\begin{array}{ccc}\hfill {\omega }^{\left[1\right]}& =& -y\left(y^2-3d^2\right)sin\theta +d\left(3y^2-d^2\right)cos\theta ,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\omega }^{\left[2\right]}& =& \left(y^2-d^2\right)\left(ysin\theta -dcos\theta \right),\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\omega }^{\left[3\right]}& =& -{\omega }^{\left[1\right]},\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill {\omega }^{\left[4\right]}& =& -2{\omega }^{\left[1\right]}+c_0^2(-dy)(dsin\theta -ycos\theta ).\hfill \end{array}$$

(31)

In the limit of $c_0\to \infty$, the right-hand side of Equation (27) is reduced to

$$2\mathcal{G}\frac{\left(y^2-d^2\right)\left(ysin\theta -dcos\theta \right)}{{\left(y^2+d^2\right)}^2},$$

(32)

which agrees with the counterpart obtained by the classical dislocation theory.

3.5. Stress from Excess Airy’s Function

For now, we have successfully derived the nonsingular stress components associated with the real and image dislocations that satisfy one of the two free surface conditions, i.e., ${\sigma }_{xx}(0,y)=0$. In order to make another remaining component ${\tau }_{xy}(0,y)$ be zero keeping ${\sigma }_{xx}(0,y)=0$, we now consider the definition of excess Airy’s function ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)$.

The two stress components, ${\sigma }_{xx}^{\left(\mathrm{ex}\right)}(x,y)$ and ${\tau }_{xy}^{\left(\mathrm{ex}\right)}(x,y)$, associated with ${\varphi }^{\left(\mathrm{ex}\right)}$ must satisfy the following relations:

$${\sigma }_{xx}^{\left(\mathrm{ex}\right)}(0,y)=0,$$

(33)

and

$${\tau }_{xy}^{\left(\mathrm{ex}\right)}(0,y)=-2\mathcal{G}\frac{4{\omega }^{\left[1\right]}+c_0^2{\omega }^{\left[2\right]}+2c_0^2{K}_0\left(c_0\right){\omega }^{\left[3\right]}+2c_0{K}_1\left(c_0\right){\omega }^{\left[4\right]}}{c_0^2{\left(y^2+d^2\right)}^2}.$$

(34)

The solution satisfying the former condition can be derived by the separation of variables method, signified by ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)=X\left(x\right)Y\left(y\right)$ (see Appendix C). A straightforward calculation leads to the conclusion that the general solution of ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)$ that satisfies ${\sigma }_{xx}^{\left(\mathrm{ex}\right)}(0,y)=0$ is

$${\varphi }^{\left(\mathrm{ex}\right)}(x,y)={\int }_0^{\infty }{a}_1\left(k\right)x{e}^{kx}sinkydk+{\int }_0^{\infty }{a}_2\left(k\right)x{e}^{kx}coskydk,$$

(35)

where $a_1\left(k\right)$ and $a_2\left(k\right)$ are k-dependent coefficients. Its partial differentiations according to Equation (1) yield the stress components associated with ${\varphi }^{\left(\mathrm{ex}\right)}(x,y)$:

$$\begin{array}{ccc}\hfill {\sigma }_{xx}^{\left(\mathrm{ex}\right)}(x,y)& =& I_{s,1}(-kx;x,y)+I_{c,2}(-kx;x,y),\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill {\sigma }_{yy}^{\left(\mathrm{ex}\right)}(x,y)& =& I_{s,1}(2+kx;x,y)+I_{c,2}(2+kx;x,y),\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill {\tau }_{xy}^{\left(\mathrm{ex}\right)}(x,y)& =& I_{s,2}(1+kx;x,y)+I_{c,1}(-1-kx;x,y),\hfill \end{array}$$

(38)

where

$$\begin{array}{ccc}\hfill I_{s,i}(u;x,y)& =& {\int }_0^{\infty }uk{a}_i\left(k\right)e^{kx}sinkydk,\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill I_{c,i}(u;x,y)& =& {\int }_0^{\infty }uk{a}_i\left(k\right)e^{kx}coskydk,\hfill \end{array}$$

(40)

with $i=1$ or 2.

The remaining task is to find the k-dependences of $a_1\left(k\right)$ and $a_2\left(k\right)$ that suffice for ${\tau }_{xy}^{\left(\mathrm{ex}\right)}(0,y)$ to cancel the shear stress contribution of Equation (27) from the real and image dislocations. From Equations (34) and (38), this requirement regarding ${\tau }_{xy}^{\left(\mathrm{ex}\right)}(0,y)$ is expressed as

$$I_{s,2}(1;0,y)+I_{c,1}(-1;0,y)=-2\mathcal{G}\frac{4{\omega }^{\left[1\right]}+c_0^2{\omega }^{\left[2\right]}+2c_0^2{K}_0\left(c_0\right){\omega }^{\left[3\right]}+2c_0{K}_1\left(c_0\right){\omega }^{\left[4\right]}}{c_0^2{\left(y^2+d^2\right)}^2},$$

(41)

with the definitions of $c_0\left(y\right)$ and ${\omega }^{\left[\mathrm{j}\right]}(y,\theta )$$(j=1,2,3,4)$ given by Equations (26) and (29)–(31), respectively. We then apply the inverse cosine and sine Fourier transforms to both sides of Equation (41), in which the integrations are performed using Cauchy’s residue theorem, which is a powerful tool for evaluating real improper integrals (from −∞ to ∞) of analytic functions. As a result of the calculation, we obtain the following solution:

$$\begin{array}{ccc}\hfill k{a}_1\left(k\right)& =& 2\mathcal{G}\left(kd-\frac{2k^2}{{\kappa }^2}\right)e^{-kd}cos\theta +\frac{4\mathcal{G}}{\pi }{\int }_{-\infty }^{\infty }\mathcal{K}(y,\theta )coskydy,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill k{a}_2\left(k\right)& =& 2\mathcal{G}\left(kd-1-\frac{2k^2}{{\kappa }^2}\right)e^{-kd}sin\theta -\frac{4\mathcal{G}}{\pi }{\int }_{-\infty }^{\infty }\mathcal{K}(y,\theta )sinkydy.\hfill \end{array}$$

(43)

with

$$\mathcal{K}(y,\theta )=\frac{c_0^2{K}_0\left(c_0\right){\omega }^{\left[3\right]}+c_0{K}_1\left(c_0\right){\omega }^{\left[4\right]}}{c_0^2{\left(y^2+d^2\right)}^2},c_0=\kappa \sqrt{y^2+d^2}.$$

(44)

The integrals in Equations (42) and (43) cannot be solved analytically owing to the presence of $K_0$ and $K_1$. Therefore, we used numerical integration to derive the k-dependences of the coefficients $k{a}_1\left(k\right)$ and $k{a}_2\left(k\right)$. Once $k{a}_1\left(k\right)$ and $k{a}_2\left(k\right)$ are numerically obtained, the spatial distribution of the stress components, ${\sigma }_{xx}^{\left(\mathrm{ex}\right)}(x,y)$, ${\sigma }_{yy}^{\left(\mathrm{ex}\right)}(x,y)$, and ${\tau }_{xy}^{\left(\mathrm{ex}\right)}(x,y)$, can be computed using Equations (36)–(40). Eventually, we obtain the nonsingular stress field of the edge dislocation in the infinitely large system with a traction-free planar surface.

4. Numerical Results

4.1. Nonsingular Stress Field of an Edge Dislocation

Figure 2 shows the contour plots of the nonsingular stress distribution in a semi-infinite elastic medium depicted in Figure 1. The panels (a)–(c) show the component of ${\sigma }_{xx}$; (d)–(f) show ${\sigma }_{yy}$; (g)–(i) show ${\tau }_{xy}$. The Burgers vector orientation is fixed to be parallel to the x-axis. The distances from the free surface to the dislocation-core measured by $\kappa d$ (i.e., the ratio of d to the characteristic length scale ${\kappa }^{-1}$) are shifted as $\kappa d=8.0$, $4.0$, and $2.0$, from left to right in the figure.

When $\kappa d=8.0$, the profiles of the stress distribution are fairly close to those realized in an infinitely large system with no boundary, in particular for the components of ${\sigma }_{xx}$ and ${\tau }_{xy}$. The profile of ${\sigma }_{xx}$ in Figure 2a is close to a figure-eight-shaped one, endowed with the inversion symmetry with respect to the x-axis. The profile of ${\tau }_{xy}(x,y)$ in Figure 2g also exhibits a multipolar symmetry, analogous to the infinite system. An exception is the profile of ${\sigma }_{xy}$, where strong amplitude stresses occur near the free surface (all but $y=0$), even though the dislocation is far from the surface.

As the dislocation approaches the surface, the stress field distribution gradually deviates from that in the infinite system. The area with a high stress magnitude gradually decreases with decreasing $\kappa d$, indicating that the presence of the free surface suppresses the near-surface field. However, regardless of the proximity of the dislocation to the surface, the divergence near the dislocation-core that occurs in the classical theory does not occur here. The stress field derived in this work is always smooth and continuous, differing from the stress field derived from the classical theory, wherein the upward and downward peak magnitudes diverge infinitely, resulting in the field discontinuity at $y=0$. The disappearance of the singularity in our results can be visually observed in the 3D plot presented in Figure 3, where the component ${\sigma }_{xx}$ under the settings of $\kappa d=4.0$ and $\theta =0$ (the same as in Figure 2b) is plotted as an example.

Figure 4 depicts contour maps, which decompose the three constituent elements of the shear component, ${\tau }_{xy}^{\left(\mathrm{re}\right)}$, ${\tau }_{xy}^{\left(\mathrm{im}\right)}$, and ${\tau }_{xy}^{\left(\mathrm{ex}\right)}$, under the settings of $\kappa d=4.0$ and $\theta =0$. The superposition of the three contour plots reproduces Figure 2h. The stress distribution drawn by ${\tau }_{xy}^{\left(\mathrm{re}\right)}$ would have been realized by the edge dislocation at $x=-d$ if there was no surface boundary at $x=0$; thus, the system was infinitely large. Similarly, the contour diagram shown by ${\tau }_{xy}^{\left(\mathrm{im}\right)}$ is the stress field created by the (image) dislocation at $x=d$ if the surface boundary does not exist. The third element, ${\tau }_{xy}^{\left(\mathrm{ex}\right)}$, is a compensating element so that the sum of the shear stress contributions from the first two elements becomes zero at the surface boundary. In fact, for the case shown in Figure 4, the large negative shear stress of ${\tau }_{xy}^{\left(\mathrm{ex}\right)}$ within the region near the origin compensates for the large positive shear stress derived from the sum of the other two elements, ${\tau }_{xy}^{\left(\mathrm{re}\right)}+{\tau }_{xy}^{\left(\mathrm{im}\right)}$ in the same region. It should be emphasized that all three elements are smooth and continuous, free from singularities near the dislocation-core.

The above-mentioned facts are valid even if we set the Burgers vector to a different orientation from the x-axis. Figure 5 shows the case of $\theta =\pi /4$, in which the Burgers vector is oblique diagonally to the free surface. Again, the obtained stress fields are free from singularities, for any choice of $\kappa d$ and $\theta$.

As previously mentioned, only component ${\sigma }_{yy}$ can have a finite magnitude at the traction-free surface at $x=0$. Figure 6 shows the spatial distribution of the stress component ${\sigma }_{yy}(0,y)$ along the free surface of $x=0$. The closer the dislocation is to the surface, the sharper are the peaks near $y=0$. As is clear from Figure 6, the regions with strong ${\sigma }_{yy}$ components do not disappear and remain, even if the dislocations are far away from the free surface. This result can be said to be a manifestation of the long-range interaction nature between dislocations and free surfaces.

4.2. Spatial Field Modulation Induced by the Free Surface

In order to estimate the magnitude of spatial modulation in the stress field caused by the presence of the zero-traction surface, we use the concept of the $L_2$-norm of the function [50]. In general, the $L_2$-norm of a function $f(x,y)$, denoted by ${\parallel f\parallel }_2$, is defined by

$${\parallel f\parallel }_2\equiv {\left\{{\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }dy{\left|f(x,y)\right|}^2\right\}}^{\frac{1}{2}}.$$

(45)

This quantity can be regarded as a generalization of the length of a vector $\mathit{x}=(x_1,x_2,\dots ,x_n)$, in an n-dimensional space:

$$\left|\mathit{x}\right|={\left(\sum _{j=1}^n{x}_j^2\right)}^{\frac{1}{2}}.$$

(46)

Using this analogy, the distance between two functions $f(x,y)$ and $g(x,y)$ is measured by the $L_2$-norm of ${\parallel f-g\parallel }_2$, which quantifies the deviation of the spatial distribution of $f(x,y)$ from that of $g(x,y)$. Based on the discussion, we calculated the following three quantities:

$${\mathcal{S}}_{\mathit{xx}}=\frac{\parallel {\sigma }_{\mathit{xx}}-{\sigma }_{\mathit{xx}}^{\left(\mathrm{re}\right)}{\parallel }_2}{\parallel {\sigma }_{\mathit{xx}}^{\left(\mathrm{re}\right)}{\parallel }_2},{\mathcal{S}}_{\mathit{yy}}=\frac{\parallel {\sigma }_{\mathit{xy}}-{\sigma }_{\mathit{xy}}^{\left(\mathrm{re}\right)}{\parallel }_2}{\parallel {\sigma }_{\mathit{xy}}^{\left(\mathrm{re}\right)}{\parallel }_2},{\mathcal{T}}_{xy}=\frac{\parallel {\tau }_{\mathit{xy}}-{\tau }_{\mathit{xy}}^{\left(\mathrm{re}\right)}{\parallel }_2}{\parallel {\tau }_{\mathit{xy}}^{\left(\mathrm{re}\right)}{\parallel }_2},$$

(47)

in all of which the denominator takes the role of “unit length” for normalization.

Figure 7 shows the results of ${\mathcal{S}}_{\mathit{xx}}$, ${\mathcal{S}}_{\mathit{yy}}$, and ${\mathcal{T}}_{\mathit{xy}}$ for the cases of $\theta =0$ and $\pi /4$. As can be seen from the figure, the $L_2$-norm curves exhibit crossover behavior across $\kappa d=1$. When the dislocation is close to the surface ($d\ll {\kappa }^{-1}$), the $L_2$-norms are nearly constant. This indicates that the norm value does not change significantly because the surface effect is already sufficiently strong even if the distance from the dislocation to the surface changes slightly. In contrast, when the dislocation is far from the surface ($d\gg {\kappa }^{-1}$), the norm decreases in a logarithmic manner. This logarithmic heavy-tail decay in the norm represents again the long-range nature of the dislocation-to-surface interaction.

5. Discussion

In general, the modified Bessel functions of the second order, $K_n\left(u\right)$, are difficult to manipulate in an analytic manner because they cannot be expressed in a closed form. However, if the value of the variable u is not less than 1, the following approximation formula holds for the cases of $n=0$ and $n=1$:

$$\begin{array}{ccc}\hfill K_0\left(u\right)& \simeq & \sqrt{\frac{\pi }{2u}}{e}^{-u}\left(1-\frac{1}{8u}\right).\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill u{K}_1\left(u\right)& \simeq & \sqrt{\frac{\pi u}{2}}{e}^{-u}\left(1+\frac{3}{8u}\right).\hfill \end{array}$$

(49)

Therefore, when considering the case where the dislocation-core is located farther than the characteristic length ${\kappa }^{-1}$ from the free surface at $x=0$, the functions $K_0$ and $K_1$ involved in the present study can be approximated by Equations (48) and (49), respectively. This approximation will facilitate the analytical treatment of Airy’s stress functions and the resulting stress fields in the systems under study.

Figure 8 visually shows the high accuracy of this approximate expression. It shows the y-dependences of the functions:

$$K_0\left(\kappa \sqrt{y^2+d^2}\right)\mathrm{and}\kappa \sqrt{y^2+d^2}{K}_1\left(\kappa \sqrt{y^2+d^2}\right),$$

(50)

both of which play an important role in the formulation of the present nonsingular stress field. Clearly, in the settings of $\kappa d=2.0$ and $4.0$, these functions agree quite well with the approximation for all y values, as drawn by dotted curves. The accuracy of the approximation increases as the value of $\kappa d$ increases. In addition, these functions exhibit large values only in a limited area near the origin (i.e., $y=0$), showing a single peak at $y=0$ and decaying sharply as they move away from $y=0$. Such a localization property, as well as the high accuracy of the approximation formula of the modified Bessel functions for $\kappa d>1$ are useful in situations where simple mathematical expressions can facilitate understanding of the process under study.

6. Concluding Remark

We developed an explicit formula of the nonsingular stress distribution around an edge dislocation near the traction-free surface. Our formula is applicable to any value of the distance from the free surface to the dislocation center and for any orientation of the Burgers vector of the edge dislocation. The formulas allow us to evaluate the impact of the free surface on the interior stress field with high accuracy without the singularity problem that occurs at the dislocation-core in existing theories.

In this study, we only considered the stress field of single edge dislocations that exist alone, with a particular emphasis on the explicit and tractable mathematical expressions. We believe that an extension of our formulation from the single-dislocation system to mutually interacting dislocations systems will provide a theoretical basis for accomplishing singularity-free analysis of energetically stable configurations of multiple dislocations near the surface boundary of nanometals. Another possible extension of this work is related to the selection of boundary conditions. The main assumption used in this study was that the boundary of the elastic medium was a free surface with a planar geometry. The zero-traction conditions were expressed by ${\sigma }_{xx}=0$ and ${\tau }_{xy}=0$ over the boundary surface of $x=0$. On the other hand, actual metal specimens are often surface-treated; to deal with the latter case, appropriate boundary conditions different from those above need to be applied in theoretical analyses. Given the generality of the formulation used in this study, it can be extended to apply to the latter systems by replacing the imposed boundary conditions with appropriate ones; this issue will be considered in the near future.