# coding: utf-8
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
# Distributed under the terms of "New BSD License", see the LICENSE file.
from functools import cached_property
from typing import Annotated
import numpy as np
from semantikon.converter import units
__author__ = "Sam Waseda"
__copyright__ = (
"Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH "
"- Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Sam Waseda"
__email__ = "waseda@mpie.de"
__status__ = "development"
__date__ = "Aug 21, 2021"
[docs]
class Dislocation:
"""
Anisotropic elasticity theory for dislocations described by
[Eshelby](https://doi.org/10.1016/0001-6160(53)90099-6).
All notations follow the original paper.
"""
def __init__(self, elastic_tensor: np.ndarray, burgers_vector: np.ndarray) -> None:
"""
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
burgers_vector ((3,)-array): Burgers vector
"""
assert np.shape(elastic_tensor) == (3, 3, 3, 3)
assert np.shape(burgers_vector) == (3,)
self.elastic_tensor = elastic_tensor
self.burgers_vector = burgers_vector
self.fit_range = np.linspace(0, 1, 10)
def _get_pmat(self, x: np.ndarray) -> np.ndarray:
return (
self.elastic_tensor[:, 0, :, 0]
+ np.einsum(
"...,ij->...ij",
x,
self.elastic_tensor[:, 0, :, 1] + self.elastic_tensor[:, 1, :, 0],
)
+ np.einsum("...,ij->...ij", x**2, self.elastic_tensor[:, 1, :, 1])
)
[docs]
@cached_property
def p(self) -> np.ndarray:
coeff = np.polyfit(
self.fit_range, np.linalg.det(self._get_pmat(self.fit_range)), 6
)
p = np.roots(coeff)
p = p[np.imag(p) > 0]
return p
[docs]
@cached_property
def Ak(self) -> np.ndarray:
Ak = []
for mat in self._get_pmat(self.p):
values, vectors = np.linalg.eig(mat.T)
Ak.append(vectors.T[np.absolute(values).argmin()])
return np.array(Ak)
[docs]
@cached_property
def D(self) -> np.ndarray:
F = np.einsum("n,ij->nij", self.p, self.elastic_tensor[:, 1, :, 1])
F += self.elastic_tensor[:, 1, :, 0]
F = np.einsum("nik,nk->ni", F, self.Ak)
F = np.concatenate((F.T, self.Ak.T), axis=0)
F = np.concatenate((np.real(F), -np.imag(F)), axis=-1)
D = np.linalg.solve(F, np.concatenate((np.zeros(3), self.burgers_vector)))
return D[:3] + 1j * D[3:]
@property
def dzdx(self) -> np.ndarray:
return np.stack((np.ones_like(self.p), self.p, np.zeros_like(self.p)), axis=-1)
def _get_z(self, positions: np.ndarray) -> np.ndarray:
z = np.stack((np.ones_like(self.p), self.p), axis=-1)
return np.einsum("nk,...k->...n", z, np.asarray(positions)[..., :2])
[docs]
def get_displacement(self, positions: np.ndarray) -> np.ndarray:
"""
Displacement vectors
Args:
positions ((n,3)-array): Positions for which the displacements are
to be calculated
Returns:
((n,3)-array): Displacement vectors
"""
return np.imag(
np.einsum(
"nk,n,...n->...k", self.Ak, self.D, np.log(self._get_z(positions))
)
) / (2 * np.pi)
[docs]
def get_strain(self, positions: np.ndarray) -> np.ndarray:
"""
Strain tensors
Args:
positions ((n,3)-array): Positions for which the strains are to be
calculated
Returns:
((n,3,3)-array): Strain tensors
"""
strain = np.imag(
np.einsum(
"ni,n,...n,nj->...ij",
self.Ak,
self.D,
1 / self._get_z(positions),
self.dzdx,
)
)
strain = strain + np.einsum("...ij->...ji", strain)
return strain / 4 / np.pi
[docs]
@units
def get_dislocation_displacement(
elastic_tensor: Annotated[np.ndarray, {"units": "=e"}],
positions: np.ndarray,
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
) -> Annotated[np.ndarray, {"units": "=b"}]:
"""
Displacement field around a dislocation according to anisotropic elasticity theory
described by [Eshelby](https://doi.org/10.1016/0001-6160(53)90099-6).
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
positions ((n,2) or (n,3)-array): Position around a dislocation. The third axis
coincides with the dislocation line.
burgers_vector ((3,)-array): Burgers vector
Returns:
((n, 3)-array): Displacement field (z-axis coincides with the dislocation line)
"""
return Dislocation(elastic_tensor, burgers_vector).get_displacement(positions)
[docs]
@units
def get_dislocation_strain(
elastic_tensor: Annotated[np.ndarray, {"units": "=e"}],
positions: Annotated[np.ndarray, {"units": "=p"}],
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
) -> Annotated[np.ndarray, {"units": "=b/p"}]:
"""
Strain field around a dislocation according to anisotropic elasticity theory
described by [Eshelby](https://doi.org/10.1016/0001-6160(53)90099-6).
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
positions ((n,2) or (n,3)-array): Position around a dislocation. The third axis
coincides with the dislocation line.
burgers_vector ((3,)-array): Burgers vector
Returns:
((n, 3, 3)-array): Strain field (z-axis coincides with the dislocation line)
"""
return Dislocation(elastic_tensor, burgers_vector).get_strain(positions)
[docs]
@units
def get_dislocation_stress(
elastic_tensor: Annotated[np.ndarray, {"units": "=e"}],
positions: Annotated[np.ndarray, {"units": "=p"}],
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
) -> Annotated[np.ndarray, {"units": "=e*b/p"}]:
"""
Stress field around a dislocation according to anisotropic elasticity theory
described by [Eshelby](https://doi.org/10.1016/0001-6160(53)90099-6).
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
positions ((n,2) or (n,3)-array): Position around a dislocation. The third axis
coincides with the dislocation line.
burgers_vector ((3,)-array): Burgers vector
Returns:
((n, 3, 3)-array): Stress field (z-axis coincides with the dislocation line)
"""
strain = get_dislocation_strain(elastic_tensor, positions, burgers_vector)
return np.einsum("ijkl,...kl->...ij", elastic_tensor, strain)
[docs]
@units
def get_dislocation_energy_density(
elastic_tensor: Annotated[np.ndarray, {"units": "=e"}],
positions: Annotated[np.ndarray, {"units": "=p"}],
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
) -> Annotated[np.ndarray, {"units": "=e*b**2/p**2"}]:
"""
Energy density field around a dislocation (product of stress and strain, cf. corresponding
methods)
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
positions ((n,2) or (n,3)-array): Position around a dislocation. The third axis
coincides with the dislocation line.
burgers_vector ((3,)-array): Burgers vector
Returns:
((n,)-array): Energy density field
"""
strain = get_dislocation_strain(elastic_tensor, positions, burgers_vector)
return np.einsum("ijkl,...kl,...ij->...", elastic_tensor, strain, strain)
[docs]
@units
def get_dislocation_energy(
elastic_tensor: Annotated[np.ndarray, {"units": "=e"}],
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
r_min: Annotated[float, {"units": "=r"}],
r_max: Annotated[float, {"units": "=r"}],
mesh: int = 100,
) -> Annotated[float, {"units": "=e*b**2"}]:
"""
Energy per unit length along the dislocation line.
Args:
elastic_tensor ((3,3,3,3)-array): Elastic tensor
burgers_vector ((3,)-array): Burgers vector
r_min (float): Minimum distance from the dislocation core
r_max (float): Maximum distance from the dislocation core
mesh (int): Number of grid points for the numerical integration along the angle
Returns:
(float): Energy of dislocation per unit length
The energy is defined by the product of the stress and strain (i.e. energy density),
which is integrated over the plane vertical to the dislocation line. The energy density
:math:`w` according to the linear elasticity is given by:
.. math:
w(r, \\theta) = A(\\theta)/r^2
Therefore, the energy per unit length :math:`U` is given by:
.. math:
U = \\log(r_max/r_min)\\int A(\\theta)\\mathrm d\\theta
This implies :math:`r_min` cannot be 0 as well as :math:`r_max` cannot be infinity. This
is the consequence of the fact that the linear elasticity cannot describe the core
structure properly, and a real medium is not infinitely large. While :math:`r_max` can
be defined based on the real dislocation density, the choice of :math:`r_min` should be
done carefully.
"""
if r_min <= 0:
raise ValueError("r_min must be a positive float")
theta_range = np.linspace(0, 2 * np.pi, 100, endpoint=False)
r = np.stack((np.cos(theta_range), np.sin(theta_range)), axis=-1) * r_min
strain = get_dislocation_strain(elastic_tensor, r, burgers_vector=burgers_vector)
return (
np.einsum("ijkl,nkl,nij->", elastic_tensor, strain, strain)
* np.diff(theta_range)[0]
* r_min**2
* np.log(r_max / r_min)
)
[docs]
@units
def get_dislocation_force(
stress: Annotated[np.ndarray, {"units": "=s"}],
glide_plane: np.ndarray,
burgers_vector: Annotated[np.ndarray, {"units": "=b"}],
) -> Annotated[np.ndarray, {"units": "=s*b"}]:
"""
Force per unit length along the dislocation line.
Args:
stress ((n, 3, 3)-array): External stress field at the dislocation line
glide_plane ((3,)-array): Glide plane
burgers_vector ((3,)-array): Burgers vector
Returns:
((3,)-array): Force per unit length acting on the dislocation.
"""
g = np.asarray(glide_plane) / np.linalg.norm(glide_plane)
return np.einsum(
"i,...ij,j,k->...k", g, stress, burgers_vector, np.cross(g, [0, 0, 1])
)