Methods to compute dislocation line tension energy and force in anisotropic elasticity

Some dislocation interactions, e.g. bow-out, cross-slip and junction formation, are well captured by computing only an approximation of the total force on the dislocation. This approximation, the orientation-dependent line tension model, represents dislocations as flexible strings. It assumes a simple energy per unit length of dislocation line and neglects the interactions between dislocations. It has been used to compute the critical stress at which a dislocation bows out [1] and to determine the shape of a dislocation glide loop at equilibrium [2] for general anisotropy. It has also proven effective in determining binary dislocation interactions as a function of their angle of incidence [3–7].

Line tension expressions for the force and the energy in anisotropic elasticity have been derived by Bacon et al in [8]. In sections 2 and 4 respectively, we revisit the formulation of energy and force in the line tension approximation and rewrite the expressions in a coordinate-independent form. The formulation proposed by Bacon et al [8] is extended to all Burgers' vectors and line directions instead of being defined only in the plane containing the line direction of the dislocation and its Burgers' vector.

The cost of computing energy in the line tension model is not negligible. A comparison is made between three methods of evaluating the energy and the force: a direct numerical integration method (section 3.1); a spherical harmonics expansion (section 3.2); and an interpolation table method (section 3.3). We study the cost and accuracy of each method as a function of the anisotropy ratio (defined as A = 2C₄₄/(C₁₁ − C₁₂) in cubic crystals where C₁₁, C₁₂ and C₄₄ are the elastic constants; A = 1 for an isotropic material) to determine which one is most suitable for implementation in large scale dislocation dynamics codes.

The strain energy per unit length of dislocation line in anisotropic elasticity is given in Bacon et al [8] by

$$\begin{equation} {\cal E} = E \ln \left( \frac{R}{a} \right), \end{equation} \tag{ 1 }$$

where R and a are the outer and inner cut-off radii, and E is the pre-logarithmic energy factor given by

$$\begin{equation} E({\bit b}, {\bit t}) = b_i B_{ij} ({\bit t}) b_j. \end{equation} \tag{ 2 }$$

E is uniquely determined by the Burgers' vector b of the dislocation line, and the 3 × 3 matrix B(t) given by

$$\begin{equation} B_{ij} ({\bit t})= \frac{1}{4 \pi^2} \int_0^{\pi} \left[ (mm)_{ij} - (mn)_{ir} (nn)_{rk}^{-1} (nm)_{kj} \right]\,{\rm d}\omega, \end{equation} \tag{ 3 }$$

where the unit vectors m and n are orthogonal to t, the dislocation line direction. The notation (mm) is such that (mm)_jk = m_iC_ijklm_l where C_ijkl is the fourth order elasticity tensor.

Similarly to Bacon et al [8], the two orthogonal unit vectors m and n can be defined by

$$\begin{eqnarray} {\bit m} & = & \cos\omega \,{\bit M} + \sin\omega \, {\bit N}\nonumber\\ {\bit n} & = & -\sin\omega \, {\bit M} + \cos\omega \, {\bit N}, \end{eqnarray} \tag{ 4 }$$

where ω is the integration variable. It varies in [0, π] (see figure 1).

Zoom In Zoom Out Reset image size

N is defined as the normalized vector

$$\begin{equation} N_i = \frac{D_{i \alpha}}{|D_{i\alpha}|}, \end{equation} \tag{ 5 }$$

where the matrix D is such that D ≡ I − t ⊗ t, where I is the identity matrix, the rows of D correspond to a set of vectors orthogonal to t. Define α as the index for which the component j of the vector v defined as

$$\begin{equation} v_j = \sqrt{\sum_i D_{ij}^2}, \end{equation} \tag{ 6 }$$

is maximum. D_α = (D_iα)_i=1,2,3 is one non-vanishing vector orthogonal to t that can always be defined.

Defining M = N × t, we have two unit vectors (M, N) orthogonal to t that are well defined for all orientations. By construction, (t, M, N) form an orthonormal basis independent of the Burgers' vector b.

In Bacon et al [8], equation (3) and its derivatives are given as angular derivatives of the line direction t with respect to an angle between the line direction t and some arbitrary fixed direction. In definitions (4)–(6), the fixed direction is D_α and it depends on t which allows M, N, m and n to always be well defined.

The complexity of calculating the energy ${\cal E}$, equation (1), resides in the definition of the matrix B, equation (3). Aside from the first term in equation (3), the matrix B cannot be computed analytically. Several numerical methods can be used to approximate its calculation. A direct approach using trapezoidal integration, spherical harmonic expansions, and interpolation methods are among the fastest techniques that can be used to evaluate B for a given set of elastic constants and dislocation line direction.

3.1. Direct numerical integration

The first term of the integral on the right-hand side of equation (3) can be determined analytically. It is given by

$$\begin{equation} \int_0^{\pi} ({\bit mm}) \, {\rm d}\omega = \frac{\pi}{2} \left[ {\bit MM} + {\bit NN} \right]. \end{equation} \tag{ 7 }$$

The other term cannot be determined analytically because it involves a quotient of polynomials of order 4 in the numerator and 6 in the denominator. It can be calculated using a numerical method such as trapezoid integration.

The accuracy of the numerical integration depends on the number of angles ω at which the second term in the integrand of B is evaluated and on the anisotropy ratio of the material. It does not depend on the length of the dislocation segment considered.

3.2. Spherical harmonics expansion methods

The matrix B defined in equation (3) may equivalently be described by an expansion in a spherical harmonics series as

$$\begin{equation} {\bit B}({\bit t}) = \sum\limits_{\ell=0}^{\infty} \sum\limits_{m=-\ell}^{\ell} {\bit B}^{\ell m} Y_\ell^m (\theta, \phi) \end{equation} \tag{ 8 }$$

which, by definition, uniformly converges on the unit sphere [9].

The spherical harmonics $Y_\ell^m$ are defined as the complex functions

$$\begin{equation} Y_\ell^m(\theta,\phi) = M_\ell^m \rme^{\rmi m\phi} P_\ell^m (\cos\theta), \end{equation} \tag{ 9 }$$

where $M_\ell^m = \sqrt{\frac{(2\ell+1)}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}}$, (θ, φ) are the spherical coordinates of t, and $P_\ell^m(\cos\theta)$ are the associated Legendre polynomials [10].

The expansion coefficients $B^{\ell m}_{ij}$ are independent of the line direction t(θ, φ) and are defined as

$$\begin{equation} B^{\ell m}_{ij} = \int\limits_0^{2\pi} \int\limits_0^\pi B_{ij}(\theta, \phi) Y_\ell^{m*}(\theta, \phi) \, \sin\theta\,{\rm d}\theta\,{\rm d}\phi. \end{equation} \tag{ 10 }$$

where $Y_\ell^{m*}$ is the complex conjugate of $Y_\ell^{m}$. The coefficients $B^{\ell m}_{ij}$ can thus be pre-calculated independently of t and stored.

An analysis of the properties of the expansion coefficients $B^{\ell m}_{ij}$ shows $B^{\ell m}_{ij} = 0$ when ℓ is odd.

3.2.1. Spherical harmonics series expansion using the recurrence relation of the associated Legendre polynomials (method I).

In equation (8), the definition of the matrix B involves associated Legendre polynomials $P_\ell^m(\cos\theta)$ (see equation (9)). A recurrence relation is often used to calculate the associated Legendre polynomials. Since only the even values of ℓ in equation (8) give non-vanishing values of B^ℓm, the following recurrence relation on even ℓ can be used

$$\begin{eqnarray} P_{\ell+2}^{\ell+2}(z) &=& (1-z^2)(2\ell+1)(2\ell+3) P_{\ell}^{\ell}(z) \nonumber\\ P_{\ell+2}^{\ell+1}(z) &=& -z\sqrt{1-z^2}(2\ell+1)(2\ell+3) P_{\ell}^{\ell}(z)\\ P_{\ell+2}^{m}(z) &=& \frac{-1}{(\ell+m+3)(\ell-m+2)} \left[\frac{2(m+1)z}{\sqrt{1-z^2}}P_{\ell+2}^{m+1}(z)+ P_{\ell+2}^{m+2}(z) \right],\nonumber \end{eqnarray} \tag{ 11 }$$

where z = cos θ and θ is one of the angles (θ, φ) of the spherical coordinates of the line direction t of the dislocation considered. The initial term of the recurrence is $P_0^0(z)=1$.

The recurrence relations for the associated Legendre polynomials, and the term e^imφ in the definition of $Y_\ell^m$, need to be computed for each dislocation segment, since they depend on the orientation of the line.

3.2.2. Spherical harmonics series expansion using pre-calculated functions (method II).

Another implementation using spherical harmonics involves pre-calculating part of the associated Legendre polynomials. This leads to another equality defining the B matrix of equation (3), see [11] for more details.

$$\begin{equation} B_{ij}({\bit t}) = \sum\limits_{\ell=0}^{\infty} \sum\limits_{m=0}^{\ell} \sum\limits_{k=0}^{[(\ell-m)/2]} \Re \left[ Q^{\ell m}_k B_{ij}^{\ell m} ({\bit t} \cdot {{\bit e}_{12}})^m ({\bit t} \cdot {{\bit e}_3})^{\ell-m-2k} \right], \end{equation} \tag{ 12 }$$

where the elastic constant matrix C_ijkl is defined in the basis (e₁, e₂, e₃) and e₁₂ is defined by e₁₂ = e₁ + ı e₂. In equation (12), ℜ[x] denotes the real part of x. The coefficients $Q^{\ell m}_k$ are defined as

$$\begin{equation} Q^{\ell 0}_k = \frac{(-1)^{k}}{4 \pi^2} \frac{1}{2^\ell} \sqrt{\frac{(2 \ell +1)}{4\pi}} \left(\begin{array}{@{}c@{}} \ell \\ k \end{array} \right) \left( \begin{array}{@{}c@{}} 2 \ell -2k\\ \ell \end{array}\right), \end{equation} \tag{ 13 }$$

when m = 0 and

$$\begin{equation} Q^{\ell m}_k = 2 \frac{(-1)^{m+k}}{4 \pi^2} \frac{m!}{2^\ell} \sqrt{\frac{(2 \ell +1)}{4\pi} \frac{(\ell-m)!}{(\ell +m)!}} \left( \begin{array}{@{}c@{}} \ell \\ k \end{array} \right) \left( \begin{array}{@{}c@{}} 2 \ell -2k\\ \ell \end{array} \right) \left( \begin{array}{@{}c@{}} \ell-2k\\ m \end{array}\right), \end{equation} \tag{ 14 }$$

when m > 0.

The definition of B in equation (12) involves a product of terms that depends only on the two variables m and ℓ − 2k and can be simplified further [11].

$$\begin{equation} B_{ij}({\bit t}) = \sum\limits_{q=0}^{\infty} \sum\limits_{m=0}^{2q} \Re \left[ S_{ij}^{q,m} ({\bit t} \cdot {{\bit e}_{12}})^m ({\bit t} \cdot {\bit e}_3)^{2q-m} \right], \end{equation} \tag{ 15 }$$

where $S^{q,m}_{ij}$ is given as

$$\begin{equation} S^{q,m}_{ij} = \sum_{2p=2q}^{\infty} Q_{p-q}^{2p,m} B^{2p,m}_{ij} \end{equation} \tag{ 16 }$$

and is independent of t, so it can be pre-calculated independently of t. Only the terms of the form (t · e₁₂)^m (t · e₃)^2q − m depend on t in this expression.

3.3. Interpolation using look-up tables

As noted earlier, the matrix B depends only on the line direction t of a given dislocation segment and not on its length. Another way of evaluating the matrix B is to pre-calculate its values on a series of grid points on the half-sphere defined by θ ∈ [0, π/2] and φ ∈ [0, 2π], and then use an interpolation scheme to approximate the value of B at the spherical coordinates of t = t(θ, φ). Grid points on the half-sphere are needed because of the symmetry property: B(t) = B(−t).

Let (θ_i, φ_j) be the uniformly distributed grid of n_θ × n_φ angles on the unit half-sphere. Let ${{\hat {\bit B}}}_{ij} = {{\bit B}}({\bit t}^{ij})$ be the evaluation of B at the grid points (θ_i, φ_j) with t^ij = (sin θ_i cos φ_j, sin θ_i sin φ_j, cos θ_i). The distance between two angles is denoted by δθ = θ₁ − θ₀ and δφ = φ₁ − φ₀.

The linear interpolation (interpolation of order d = 1) for B(θ, φ) is

$$\begin{eqnarray} &&\fl {\bit B}(\theta, \phi) = (1-p)(1-q){{\hat {\bit B}}}_{ij} + (1-p)q {{\hat {\bit B}}}_{i,,j+1}+ p(1-q) {{\hat {\bit B}}}_{i+1,,j} + p q{{\hat {\bit B}}}_{i+1,j+1}, \end{eqnarray} \tag{ 17 }$$

where $i=[\frac{\theta}{\delta\theta}]$ and $j=[\frac{\phi}{\delta\phi}]$, $p=\frac{\theta}{\delta\theta}-i$ and $q=\frac{\phi}{\delta\phi}-j$, with $[\frac{\theta}{\delta\theta}]$ the largest previous integer of $\frac{\theta}{\delta\theta}$.

For the same i, j, p and q, the interpolation of order d of B(θ, φ) is given by

$$\begin{equation} {\bit B}(\theta, \phi) = {\bit v}(p) {\bit A}^{-1} {{\hat {\bit B}}^{\rm g}} {\bit A}^{-1'} {\bit v}(q), \end{equation} \tag{ 18 }$$

where the vector v is defined by v_m(p) = p^d − m+1, and the notation A^−1' is the transpose of the inverse of A.

The matrix A is defined by

$$\begin{equation} A_{m, (d+2-p)} = m^{p-1}, \end{equation} \tag{ 19 }$$

where p and m vary in [1, d + 1] and in [(1 − d)/2, (1 + d)/2] respectively. The matrix ${{\hat {\bit B}}^g}$ is the (d + 1) × (d + 1) matrix defined as

$$\begin{equation} {{\hat {\bit B}}^g}_{mn} = {{\hat {\bit B}}}_{i+n, j+m} \end{equation} \tag{ 20 }$$

where m and n vary in [(1 − d)/2, (1 + d)/2].

Interpolation methods require only a few matrix-vector operations which make them fast. Their main drawback is the necessity of storing the grid of values ${{\hat {\bit B}}}_{ij}$, a set of n_θ × n_φ values.

In the dislocation dynamics code ParaDiS [12], dislocations are discretized into straight segments, and forces are evaluated at the end nodes of these segments. In a line tension calculation, the force coming from m = 1, ..., N dislocation segments (x^m, xⁿ), evaluated at the end node xⁿ is defined as

$$\begin{equation} {\bit F}^n= \sum_{m=1}^N {\bit F}^{nm} \ln \left( \frac{R}{a} \right), \end{equation} \tag{ 21 }$$

where F^nm is the pre-logarithmic force factor evaluated at xⁿ. F^nm is defined as the derivative of the energy with respect to this node position. For a dislocation segment length L^nm = |xⁿ − x^m|, it is given by

$$\begin{equation} {\bit F}^{nm} = \frac{{\rm d} E^{nm}}{{\rm d}{\bit t}^{nm}} - \left( \frac{{\rm d} E^{nm}}{{\rm d} {\bit t}^{nm}} \cdot {\bit t}^{nm} \right) {\bit t}^{nm} + E^{nm} {\bit t}^{nm}, \end{equation} \tag{ 22 }$$

where ${\bit t}^{nm} = \frac{({\bit x}^n - {\bit x}^m)}{L^{nm}}$ is the dislocation segment line direction.

In terms of the matrix B, the pre-logarithmic components of the force factor F^nm is

$$\begin{equation} F_k = b_i \left[ \frac{{\rm d} B_{ij}}{{\rm d} t_k} - \left( \frac{{\rm d} B_{ij}}{{\rm d} t_l} t_l \right) t_k + B_{ij} t_k \right] b_j, \end{equation} \tag{ 23 }$$

where superscripts (m, n) were dropped for clarity.

To define the line tension force, we need to compute the pre-logarithmic force factor F^nm defined in equations (22) and (23).

4.1. The derivative of energy with respect to the line direction

The line tension force involves the third order tensor d B/dt which needs to be determined in the same basis (t, M, N) described before (figure 1). It is given by

$$\begin{eqnarray} &&\fl \frac{{\rm d} B_{ij}}{{\rm d} t_l} = \frac{1}{4 \pi^2} \int_0^{\pi} \frac{{\rm d} (mm)_{ij}}{{\rm d} t_l} {\rm d}\omega - \frac{1}{4 \pi^2} \int_0^{\pi} \left[ \left\{\left( \frac{{\rm d} m}{{\rm d} t_l} n \right)_{ir} + \left( m \frac{{\rm d} n}{{\rm d} t_l} \right)_{ir} \right\} (nn)_{rk}^{-1} (nm)_{kj} \right.\\ &&-(mn)_{ir} \left( (nn)^{-1} \frac{{\rm d} (nn)}{{\rm d}t_l} (nn)^{-1} \right)_{rk} (nm)_{kj} \nonumber\\ &&\left. +(mn)_{ir} (nn)_{rk}^{-1} \left( \left\{\frac{{\rm d} n}{{\rm d} t_l} m \right\}_{kj} + \left\{m \frac{{\rm d} n}{{\rm d} t_l}\right\}_{kj} \right) \right]\,{\rm d}\omega.\nonumber \end{eqnarray} \tag{ 24 }$$

The derivatives of the vectors m and n are given in terms of the derivatives of M and N with respect to the line direction by

$$\begin{eqnarray} \frac{{\rm d}{\bit m}}{{\rm d}{\bit t}} & = & \cos\omega \,\frac{{\rm d}{\bit M}}{{\rm d}{\bit t}} + \sin\omega \, \frac{{\rm d}{\bit N}}{{\rm d}{\bit t}}\nonumber\\ \frac{{\rm d}{\bit n}}{{\rm d}{\bit t}} & = & -\sin\omega \, \frac{{\rm d}{\bit M}}{{\rm d}{\bit t}} + \cos\omega \, \frac{{\rm d}{\bit N}}{{\rm d}{\bit t}}, \end{eqnarray} \tag{ 25 }$$

and those of M and N by

$$\begin{eqnarray} \frac{{\rm d} N_i}{{\rm d} t_j} & = & \frac{N_i N_j t_\alpha}{\sqrt{1-t_{\alpha}^2}} -\frac{(\delta_{ij}t_{\alpha}+ \delta_{\alpha j} t_i)}{\sqrt{1-t_{\alpha}^2}} \\ \frac{{\rm d} M_i}{{\rm d} t_j} & = & \frac{\epsilon_{i\alpha j} - M_i N_j t_\alpha}{\sqrt{1-t_{\alpha}^2}}.\nonumber \end{eqnarray} \tag{ 26 }$$

Note that in the previous expressions (equations (26)), there is no sum on the index α as it is fixed by equation (6) and t_α is the α component of the line direction t.

Combining equations (25) and (26), the derivative of the matrix B given in equation (24) is fully determined.

Note that, like in the case of the energy, the first term in the integral in equation (24) can be computed analytically

$$\begin{equation} \fl \int_0^{\pi} \frac{{\rm d} ({\bit mm})}{{\rm d}{\bit t}}\,{\rm d}\omega = \frac{\pi}{2} \left[ \left(\frac{{\rm d}{\bit M}}{{\rm d}{\bit t}} {\bit M}\right) + \left({\bit M} \frac{{\rm d}{\bit M}}{{\rm d}{\bit t}}\right) +\left(\frac{{\rm d}{\bit N}}{{\rm d}{\bit t}} {\bit N}\right) + \left({\bit N} \frac{{\rm d}{\bit N}}{{\rm d}{\bit t}}\right) \right]. \end{equation} \tag{ 27 }$$

4.2. Spherical harmonics expansion of the line tension force

The line tension force, like the energy, can be decomposed in spherical harmonics. Two methods can be used to express it: one uses the recurrence relation for the associated Legendre polynomials, and the other pre-calculates as many terms in the expansion as possible to reduce the cost of its calculations.

4.2.1. Spherical harmonics series expansion using recurrence of the associated Legendre polynomials, (method I).

The chain rule for the derivative of the spherical harmonics expansion of B, equation (8) can be used to derive a recurrence relation for the line tension force, equation (23). It is given by

$$\begin{equation} \frac{{\rm d}{\bit B}}{{\rm d}{\bit t}} = \sum\limits_{\ell=0}^{\infty} \sum\limits_{m=-\ell}^{\ell} {\bit B}^{\ell m} \left[ im Y_\ell^m(t) \frac{\partial \phi}{\partial {\bit t}} + M_\ell^m \rme^{\rmi m\phi} P_\ell^{m'} (z) \frac{\partial z}{\partial {\bit t}} \right], \end{equation} \tag{ 28 }$$

where we have posed z = cos θ and where $P_\ell^{m'}(z)$ is the derivative of the Legendre polynomial $P_\ell^m(z)$ with respect to z.

Using the following identity for the associated Legendre polynomials

$$\begin{equation} (z^2-1)P_\ell^{m'}(z) = \sqrt{1-z^2} P_\ell^{m+1}(z) + mz P_\ell^{m}(z), \end{equation} \tag{ 29 }$$

and that the derivatives of the spherical coordinate φ and z = cos θ are

$$\begin{equation} \frac{\partial \phi}{\partial {\bit t}} = \frac{({\bit t} \cdot {\bit e}_1){\bit e}_2 -({\bit t} \cdot {\bit e}_2){\bit e}_1}{({\bit t} \cdot {\bit e}_1)^2 + ({\bit t} \cdot {\bit e}_2)^2} = \frac{\cos\phi \, {\bit e}_2 - \sin\phi \, {\bit e}_1}{\sin\theta} \end{equation} \tag{ 30 }$$

and

$$\begin{equation} \frac{\partial z}{\partial {\bit t}} = -(\sqrt{1-z^2}) {\bit e}_3, \end{equation} \tag{ 31 }$$

we obtain

$$\begin{equation} \fl \frac{{\rm d}{\bit B}}{{\rm d}{\bit t}} = \sum\limits_{\ell=0}^{\infty} \sum\limits_{m=-\ell}^{\ell} {\bit B}^{\ell m} M_\ell^m \left[ \left( \rmi \frac{\partial \phi}{\partial {\bit t}} - \frac{z}{\sqrt{1-z^2}} {\bit e}_3 \right) m P_\ell^{m}(z) - P_\ell^{m+1}(z) {\bit e}_3 \right] \rme^{\rmi m\phi} \end{equation} \tag{ 32 }$$

The same relation of recurrence as equations (12) can be used to compute the polynomials $P_\ell^m$ and the calculation of $\frac{{\rm d}{\bit B}}{{\rm d}{\bit t}}$ follows.

4.2.2. Spherical harmonics series expansion using pre-calculated functions (method II).

Differentiating the spherical harmonic expansion for B in equation (15) gives

$$\begin{eqnarray} &&\fl \frac{\rmd B_{ij}}{\rmd{\bit t}} = \sum\limits_{q=1}^{\infty} \sum\limits_{m=0}^{2q} \Re \left\{S^{q,m}_{ij} ({\bit t} \cdot {{\bit e}_{12}})^{m-1} ({\bit t} \cdot {\bit e}_3)^{2q-m-1} [m ({\bit t} \cdot {\bit e}_3) {{\bit e}_{12}} + (2q-m)({\bit t} \cdot {{\bit e}_{12}}) {\bit e}_3) ] \right\}. \nonumber\\ \end{eqnarray} \tag{ 33 }$$

Combining equations (15) and (33), the line tension force given in equation (23) is fully determined.

When the direct numerical integration method is used, a relative error of 10⁻¹⁶ in energy calculations is found for approximately 100 integration angles for any anisotropy ratio. The corresponding energy value is set as a reference value for estimating the relative error in energy calculations for the other methods. The corresponding time to compute this reference energy value is set to t_ref. In the next figures, 'relative cost' means the computational cost of the method considered t divided by the reference time t_ref, relative cost = t/t_ref.

The relative cost of computing the energy for the three methods described in section 3 for a fixed anisotropy ratio of A = 7.45 (corresponding to α-Fe at 900 °C) is shown in figure 2. Figure 2(a) shows the relative cost for the direct integration method (section 3.1) as well as a linear fit. The relative cost depends on the number of angles ω in [0, π] chosen to evaluate the integral. Figure 2(b) shows the relative cost as a function of the expansion order q_max for the two spherical harmonics methods (sections 3.2.1 and 3.2.2) as well as a quadratic fit.

Zoom In Zoom Out Reset image size

The relative cost of the direct method increases linearly as a function of the number of angles, and the relative cost increases quadratically as a function of the expansion order q_max when using the spherical harmonics methods. The spherical harmonics method I (section 3.2.1) is 1.85–2.7 times more expensive than the spherical harmonics method II (section 3.2.2).

The speed increase from the spherical harmonics method I to method II is a result of the larger number of precalculations done in the spherical harmonics expansion method II. In method I, only the coefficients $B^{lm}_{ij}$ are pre-computed. In method II, both the coefficients $B_{ij}^{\ell m}$ and $Q^{\ell m}_k$ are pre-calculated through the coefficients $S^{qm}_{ij}$. Also, the recurrence relations to calculate the Legendre polynomial are fast but the recurrence relation to obtain the product (t · e₁₂)^m (t · e₃)^2q−m is faster for lower order products.

The cost of computing the energy in the interpolation method (section 3.3) should depend on the number of grid points n_θ × n_φ and the order of interpolation. Figure 2(c) shows the relative cost of the interpolation method averaged over the number of angles as well as the standard deviation and a linear fit.

The relative cost of the interpolation method does not depend on the anisotropy ratio nor on the number of grid points chosen to interpolate B, since these values are pre-tabulated. Also, it does not appear to vary much with the order of interpolation. As the order of interpolation increases, larger matrices and vectors are involved in the matrix/vector products of equation (18) but this is negligible compared to the cost of retrieving the elements of the interpolated grid.

Figure 3(a) shows the error in energy calculations as a function of interpolation order and number of grid points for a fixed anisotropy ratio of 7.45. Figure 3(b) shows the base 10 log of the energy calculations error as a function of the relative cost for a fixed anisotropy ratio and a fixed interpolation order of 11. The interpolation method of order 11 needs about 2000 grid points to reach an error in energy calculations of 10⁻⁷ and about 25 000 points to reach an error of 10⁻¹². The interpolation method is fast compared to the two other methods but requires a large stored grid to converge. The relative cost of computing the energy using linear interpolation is about 20 times less expensive than the direct method.

Zoom In Zoom Out Reset image size

Figure 4 compares the error in energy calculations as a function of cost and anisotropy ratio for the direct integral calculations (figure 4(a)), the spherical harmonics method I (figure 4(b)) and II (figure 4(c)) described in the previous sections and for the interpolation method (figure 4(d)). The relative cost represented in figure 4 is deduced from figure 2.

Zoom In Zoom Out Reset image size

The interpolation method converges faster to a lower error in the energy calculation than in all the other methods. The spherical harmonics method II is the second fastest and requires much less storage. The error in energy calculation decreases exponentially for the direct integration and the spherical harmonics methods as a function of anisotropy ratio while the error in energy calculations for the interpolation method remains constant.

We have also computed the cost and accuracy of calculating the line tension force at a dislocation node for the three numerical methods: direct integration, spherical harmonic expansions and interpolation as a function of anisotropy ratio. The same conclusions seen for energy can be made for forces. Aside from the actual values, the trends are the same. The comparison results for the force error calculations are not shown in this paper.

A coordinate-independent line tension formulation for dislocation energy and force has been proposed for anisotropic elasticity. The expressions for the energy and its derivative are extended and are valid for dislocation line directions that can vary anywhere on the unit sphere. They are not limited to a plane containing the line direction and a fixed direction.

Three methods have been investigated to compute the energy in the line tension model in an anisotropic elastic medium. Choosing which method is the most cost efficient and accurate depends on the problem to be solved and the accuracy desired.

The direct integration method needs about 40–50 angles about the target line direction to reach a precision of 10⁻⁷–10⁻¹⁶ in energy, depending on the anisotropy ratio. It does not use any stored values, except perhaps the sines and cosines of the integration angles, which can be pre-computed and stored, and elastic constants.

The two spherical harmonics methods converge faster per unit computational cost than the direct method. As the anisotropy ratio moves further away from one, the spherical harmonics method becomes more expensive since the accuracy grows quadratically with q_max. It also needs to store more data than the direct integration method (such as the $S^{q,m}_{ij}$ coefficients described in the text). These coefficients represent 9(q_max + 1)(q_max + 2) real values. For q_max = 10, it is about 1200 values and for q_max = 20, this is about 4200 values.

The least computationally expensive method is the high order interpolation method. It requires approximately 1/20 of the relative cost to reach an accuracy of about 10⁻¹³ in energy error calculations but it requires the storage of angular grid points (24 500 values corresponding to 90 × 270 angle pairs), the storage of the interpolation matrix A (equation (19)) and evaluations of B to converge to this precision. The interpolation method can be used if storage of the grid points is not a problem, for instance the case of small simulations involving only a few dislocations on a single processor.

There are advantages and drawbacks in using any of these methods, and it might be useful to switch from one method to another for a given problem. The spherical harmonics method (method II) is a good compromise between low relative cost and storage.

Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the US Department of Energy, National Nuclear Security Administration under Contract DE-AC52-07NA27344.