1 Introduction and Background↩︎

The flow stress and thereby the deformation rate of crystalline materials such as metals, is strongly influenced by the glide speed of dislocations [1]–[12]. If dislocations reach a limiting velocity, more dislocations need to be generated to accommodate higher deformation rates. Even if transsonic or supersonic dislocation glide is possible, the limiting velocities constitute a hard-to-overcome barrier. The limiting velocities can be determined from the continuum-limit theory of dislocations which predicts divergences in the elastic energy at crystal-geometry dependent limiting velocities \(v_L\), see review article [13]. Dislocation-core regularization [14], [15] and lattice treatments [16] have shown that these limiting velocities can be overcome in principle and molecular dynamics simulations have predicted the existence of supersonic (single) dislocations [17]–[22] .

Recent measurements of diamond have also shown supersonic dislocations to exist [23], and with the recent advances in femto-second x-ray radiography [24], [25], this might be just the first of several future discoveries of transsonic dislocations in crystals. Aside from the questions of how dislocations can overcome (even soft) gliding velocity barriers, it is essential to know the values of these limiting velocities, as the flow stress will increase non-linearly close to them [26]. In a larger crystal plasticity finite element (CPFE) or discrete dislocation dynamics (DDD) simulation, limiting velocities will need to be calculated fast and on-the-fly as they depend on the material density and elastic constants, and are thereby temperature and pressure dependent.

Based on the so-called integral method [27], [28] of computing the dislocation displacement field, Barnett et al. [29] presented a computationally efficient method on how to determine dislocation limiting velocities in general anisotropic crystals for arbitrary slip system and dislocation character angle. In general, there are three (slip system and dislocation character-dependent) gliding velocities where the dislocation self-energy diverges. The lowest of these limiting velocities separates the ‘sub-sonic’ from the ‘transonic’ regime of dislocation glide. However, Barnett’s method does not account for special cases where the plane perpendicular to the dislocation line is a reflection plane. In this case, the differential equations for edge and screw dislocations decouple. Therefore the screw dislocation has only one limiting velocity and the edge dislocation has two limiting velocities. Analytic solutions exist for the limiting velocity of screw dislocation in this case, as well as for the edge dislocation if elastic constants \(c_{16}=0=c_{26}\) in the rotated coordinate system where the \(\hat{z}\) axis is aligned with the dislocation and the \(\hat{y}\) axis is aligned with the slip plane normal, see [30]–[34] as well as the review article [13].

The one case that is not well covered is an edge dislocation with reflection symmetry and non-vanishing \(c_{16}\) and/or \(c_{26}\). A python code to calculate the limiting velocity in this case was previously implemented in the present author’s code PyDislocDyn [35], [36] and the method was based on Teutonico’s 1961 paper [31] using a combination of symbolic (sympy) and numeric methods.

In this short letter, we derive a specialized version of Barnett’s 1973 algorithm [29] adapted to the decoupled two-dimensional case of an edge dislocation with reflection symmetry. This improved method is almost two orders of magnitude faster to compute than the previous one that was based on Teutonico’s work [31].

We start by reviewing the general 3-dimensional Barnett algorithm which we subsequently simplify for the 2-dimensional case of an edge dislocation with reflection symmetry. It must also be noted that in the case of reflection symmetry, there are subtle cancellations in the numerator of the dislocation field which depend on polar angle \(\phi\) [13], [34]. Therefore, the true limiting velocity can be higher than min\((v_\text{limit}(\phi))\) computed from the full 3-dimensional determinant. In other words, the Barnett method fails in these special cases, which showcases the necessity of a 2-dimensional version specialized to the case of an edge dislocation with reflection symmetry. Despite the simplicity of the resulting equations, this new method has not appeared in the literature to the author’s best knowledge.

2 Review of Barnett’s method↩︎

The dislocation displacement field \(u_i\) follows from the differential equations \[\begin{align} \partial_i \sigma_{ij} &= \rho \ddot{u}_j \,, & \sigma_{ij} &= C_{ijkl} \partial_l u_{k} \,,\label{eq:diffeqns1} \end{align}\tag{1}\] where \(C_{ijkl}\) is the tensor of second order elastic constants (SOEC). The established numerically efficient method of computing \(u_i\) in the steady state limit, where the dislocation depends on time only via the linear combination \(x-vt\), is the so-called integral method [27], [28]: \[\begin{align} u_j(r,\phi)&=-\frac{b_l}{2\pi}\left\{S_{jl}\ln r-S_{il}\int_0^\phi\left[(nn)^{-1}(nm)\right]_{ji}d\phi'-B_{il}\int_0^\phi(nn)^{-1}_{ji}d\phi'\right\} \,,\nonumber\\ \mat{S}&=-\frac{1}{2\pi}\int_0^{2\pi}(nn)^{-1}(nm)\,d\phi \,,\nn\\ \mat{B}&=-\frac{1}{2\pi}\int_0^{2\pi}\left[(mn)(nn)^{-1}(nm)-(mm)\right]d\phi \,, \label{eq:uj-sol} \end{align}\tag{2}\] with the short-hand notation \((ab):=\left(a_iC_{ijkl}b_j - \rho v_ia_i\delta_{jk}v_lb_l\right)\). We note that all expressions in 2 are \(\pi\)-periodic in \(\phi\). Vectors \(m_i\) and \(n_i\) are orthonormal to one another and lie in the plane perpendicular to the dislocation line; for details we refer to the excellent review article of Bacon, Barnett, and Scattergood [27].

Since above depends on the inverse of matrix \((nn)\), a general mixed-character dislocation exhibits a divergence in its elastic self-energy whenever the following determinant is zero: \[\begin{align} \det \left[n_i C_{ijkl}n_l - \rho \left(v\cos\phi\right)^2\delta_{jk}\right]=0 \,, \end{align}\] where \(\phi\) is the angle between velocity vector \(\vec{v}\) and \(\vec{n}\). The latter vector lies in the plane perpendicular to the dislocation line and \(\phi\) is the polar angle in that plane. The determinant above can be solved in terms of \(v(\phi)\) and has three \(\phi\)-dependent solutions in general, which are the roots of a cubic polynomial [29]: \[\begin{align} &x^3+bx^2+cx+d=0 \,,\nonumber\\ &b = -\tr\left(\vec{n}\cdot\mat{C}\cdot\vec{n}\right)=-\sum_{j=1}^3n_i C_{ijjl}n_l \,,\nonumber\\ &c = \frac{1}{2}\left[b^2-\tr\left(\left(\vec{n}\cdot\mat{C}\cdot\vec{n}\right)^2\right)\right] = \frac{1}{2}\left[b^2-n_i C_{ijjl}n_ln_k C_{kjjr}n_r\right] \,,\nonumber\\ &d = -\det\left(\vec{n}\cdot\mat{C}\cdot\vec{n}\right) \,. \end{align}\] The transformation \(x=y-b/3\) eliminates the quadratic term leading to \[\begin{align} &y^3 + 3py + 2q=0 \,,\nonumber\\ p &= \frac{1}{3}\left(c-b^2/3\right)\,,& q &= \frac{d}{2}-\frac{bc}{6}+\frac{b^3}{27} \,. \end{align}\] The solutions, which are real for the elasticity problem at hand (and can be written in different ways), are \[\begin{align} y &= 2\sqrt{-p}\cos\left(\frac{\psi+2k\pi}{3}\right) & k=0,1,2 \,,\nonumber\\ \cos\psi &= \frac{-q}{\sqrt{-p^3}} \,,\nonumber\\ \rho\left(v\cos\phi\right)^2&=x=y-b/3 \,. \end{align}\] Thus the three limiting velocities of a dislocation of mixed character angle is determined by minimizing the three branches \(v(\phi)\) over the interval \([0,\pi]\) (since \(\vec{n}\) is a function of \(\phi\) and several terms in \(\pi\)-periodic are integrated over \(\phi\)).

3 Dislocations with reflection symmetry↩︎

If the plane perpendicular to the dislocation line is a reflection plane, the differential equations for screw and edge dislocations decouple. Hence, one is left with a single differential equation for pure screw dislocations and two coupled differential equations for pure edge dislocations. The equation for the screw dislocation can be solved analytically for the single limiting velocity [13], [31], [34]. An analytic solution for the two limiting velocities of the edge dislocation exists only if additionally the elastic constants \(c_{16}=0=c_{26}\) vanish, after rotating into a coordinate system where the \(\hat{z}\) axis is aligned with the dislocation and the \(\hat{y}\) axis is aligned with the slip plane normal [13], [31]. The more general case must be solved numerically and a method following Teutonico’s 1961 paper [31] was put forward in [13] and implemented in python code PyDislocDyn [35], [36]. However that implementation was fairly slow due to its use of a combination of symbolic (sympy) and (scipy) minimization routines.

On the other hand, the 3-dimensional general method of Barnett reviewed above in Section 2 fails in this high-symmetry case because the terms in the numerator of the dislocation field cancel out the divergence due to the vanishing determinant for certain of polar angles \(\phi\). In other words, it becomes unclear over what (limited) range of polar angles \(\phi\) one must minimize in this case. A prime example are the {112} slip planes in body-centered cubic (bcc) crystals, where reducing the minimization-interval to \([\pi/2,\pi]\) leads to the correct lowest limiting velocity for edge dislocations — at least for the handful of cases checked by the author.

A purely numerical solution along the lines of Section 2, but reduced to the two-dimensional problem at hand is therefore in order. The first step is to rotate all quantities into a coordinate system aligned with the dislocation. The rotation matrix we need for this purpose is determined by stacking the following three row vectors into a matrix, i.e.: \[\begin{align} U&=\left(\begin{array}{c} \hat{m}_0^T \\ \hat{n}_0^T \\ \hat{t}^T \end{array}\right) \,, \label{eq:rotationmatrix} \end{align}\tag{3}\] where \(\hat{m}_0\) is the slip direction, \(\hat{n}_0\) is the slip plane normal, and \(\hat{t}\) the edge dislocation line direction. Hence, \(U\cdot \hat{m}_0=\hat{x}\), \(U\cdot \hat{n}_0=\hat{y}\), and \(U\cdot \hat{t}=\hat{z}\). The tensor of second order elastic constants, which is always measured in crystal coordinates, is rotated into coordinates aligned with the dislocation via \[\begin{align} C'_{ijkl} = U_{ii'}U_{jj'}U_{kk'}U_{ll'}C_{i'j'k'l'} \,. \end{align}\] The edge dislocation displacement field \(u_i=U_{ij}\tilde{u}_j\) in this coordinate system has only two non-vanishing components, \(u_1\) and \(u_2\) (i.e. displacements only in the \(x\) and \(y\) directions). Furthermore, \[\begin{align} \vec{n}'=U_{ij}n_j = \left(\begin{array}{c} \cos\phi\\\sin\phi\\0 \end{array}\right) \end{align}\] because of 3 . The determinant of the 2-dimensional matrix we must consider for the pure edge dislocation is hence \[\begin{align} \det \left[n'_i C'_{i\alpha\beta j}n'_j - \rho \left(v\cos\phi\right)^2\delta_{\alpha\beta}\right]=0 \,, \end{align}\] where \(\alpha\), \(\beta\in[1,2]\) and \(i,j\in[1,2,3]\). The resulting quadratic polynomial equation in \(X=\rho \left(v\cos\phi\right)^2\) is \[\begin{align} &X^2-2pX+q=0 \,,\\ p&= \frac{1}{2}\left(n'_i C'_{i11 j}n'_j +n'_i C'_{i22 j}n'_j\right) \,,\\ q &= n'_i C'_{i11 j}n'_j n'_k C'_{k22 l}n'_l - n'_i C'_{i12 j}n'_j n'_k C'_{k21 l}n'_l \,, \end{align}\] with solutions \(X = p\pm\sqrt{p^2-q}\) (which are both real for the elasticity problem at hand). Hence, the two limiting velocities for pure edge dislocations are \[\begin{align} (v_L)^2 &= \min\left[\frac{p}{\rho\cos^2\phi}\pm\sqrt{p^2-q}\right] \,, \end{align}\] which must be numerically minimized over \(\phi\in[0,\pi]\) (since all expressions involved are \(\pi\)-periodic).

The new method has been implemented in PyDislocdyn 1.3.5 [35], [36] and has been verified to yield the same result as the slower previous method of Refs. [13], [31]. As an additional consistency check, one can calculate the line tension [37] and/or drag coefficient from phonon wind [38], [39] close to the limiting velocity to see that indeed they quickly increase in value as the correct limiting velocity is approached.

4 Conclusion↩︎

We have derived a computationally efficient fast method of determining the limiting velocities of gliding edge dislocations with reflection symmetry, but non-vanishing elastic constants \(c'_{16}\) or \(c'_{26}\) — the one special case where an efficient method was previously lacking in the literature.

A reference implementation of how to compute limiting velocities of gliding dislocations for all cases is given within the open source code PyDislocDyn [35], [36] developed by the present author. The method presented here for pure edge dislocations in anisotropic crystals where the plane perpendicular to the dislocation is a reflection plane, constitutes the last “missing piece” needed to calculate all dislocation limiting velocities for all slip systems and dislocation characters numerically efficient and without the need for any symbolic (sympy) routines. It is therefore now possible to implement on-the-fly calculations of dislocation limiting velocities into larger simulations codes. An implementation in Fortran (in addition to Python) has therefore been included in the current development version of PyDislocDyn, i.e. a Fortran library (or frontend) can be compiled and linked to (called from) other codes where speed is essential.

Acknowledgements↩︎

The author is grateful for the support of the Materials project within the Advanced Simulation and Computing, Physics and Engineering Models Program of the U.S. Department of Energy under contract 89233218CNA000001.

References↩︎