This trial describes how to create edge and screw dislocations in iron BCC strating with one unitcell containing two atoms
The elastic solution for displacement field of dislocations is provided in the paper Dislocation Displacement Fields in Anisotropic Media.
The paper mentioned in backgroud subsection deals with only one dislocation. Here we describe how to extend the solution to periodic array of dislocations. Since we are dealing with linear elasticity we can superpose (sum up) the displacement field of all the individual dislocations. Looking at the Eqs. (2-8) of abovementioned reference this boils done to finding a closed form soloution for
$$\sum_{m=-\infty}^{\infty} \log\left(z-ma \right).$$Where $z= x+yi$ and $a$ is a real number, equivakent to $\mathbf{H}_{00}$ that defines the periodicity of dislocations on x direction.
Let us simplify the problem a bit further. Since this is the component displacement field we can add or subtract constant term so for each $\log\left(z-ma \right)$ we subtract a factor of $log\left(a \right)$, leading to
$$\sum_{m=-\infty}^{\infty} \log\left(\frac{z}{a}-m \right).$$Lets change $z/a$ to $z$ and when we arrive the solution we will change ot back
$$\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$Objective is to find a closed form solution for
$$f\left(z\right)=\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$First note that
$$ f'\left(z\right)=\frac{1}{z}+\sum_{m=1}^{\infty}\frac{1}{z-m}+\frac{1}{z+m}, $$and also $$ \frac{1}{z\mp m}=\mp \frac{1}{m}\sum_{n=0}^{\infty} \left(\pm \frac{z}{m}\right)^n. $$
This leads to $$ \frac{1}{z-m}+\frac{1}{z+m}=-\frac{2}{z}\sum_{n=1}^{\infty}\left(\frac{z}{m}\right)^{2n}, $$ and subsequently $$ f'\left(z\right)=\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\sum_{m=1}^{\infty}m^{-2n}, $$ $$ =\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\zeta\left(2n\right). $$ Where $\zeta$ is Riemann zeta function. Since $\zeta\left(0\right)=-1/2$, it simplifies to: $$ f'\left(z\right)=-\frac{2}{z}\sum_{n=0}^{\infty}\left(z\right)^{2n}\zeta\left(2n\right) $$ Note that $$ -\frac{\pi z\cot\left(\pi z\right)}{2}=\sum_{n=0}^{\infty}z^{2n} \zeta\left(2n\right) $$
I have no idea how I figured this out but it is true. Therefore,
$$ f'\left(z\right)=\pi\cot\left(\pi z\right). $$At this point one can naively assume that the problem is solved (like I did) and the answer is something like: $$ f\left(z\right)=\log\left[\sin\left(\pi z\right)\right]+C, $$ Where $C$ is a constant. However, after checking this against numerical vlaues you will see that this is completely wrong.
The issue here is that startegy was wrong at the very begining. The sum of the displacelment of infinte dislocations will not converge since we have infinite discountinuity in displacement field. In other words they do not cancel each other they feed each other.
But there is still a way to salvage this. Luckily, displacement is relative quantity and we are dealing with crystals. We can easily add a discontinuity in form an integer number burger vectors to a displacement field and nothing will be affected.
So here is the trick: We will focus only on the displacement field of one unit cell dislocation (number 0). At each iteration we add two dislocation to its left and right.
At $n$th iterations we add a discontinuity of the form
$$ -\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] \pi i $$and a constant of the form: $$ -2\log n. $$
In other words and we need to evaluate: $$ \lim_{m\to\infty}\sum_{n=-m}^{m} \biggl\{ \log\left(z-n\right) -\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] \pi i -2\log\left(n \right) \biggr\} + \pi, $$
which simplifies to $$ \lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-n\right) -\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i -2\log\left(\frac{m\!!}{\sqrt{\pi}} \right) $$
Note that we added an extra $\pi$ to displacement field for aesthetic reasons. After a lot of manipulations and tricks (meaning I dont't remember how I got here) we arrive at the following relation: $$ \lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-n\right) -\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i -2\log\left(\frac{m\!!}{\sqrt{\pi}} \right)=\log\left[\sin\left(\pi z\right)\right] $$ However, this is only valid when $$-1/2 \le\mathrm{Re}\left(z\right)\lt 1/2.$$
If one exceeds this domain the answer is:
$$ \boxed{ \log\left[\sin\left(\pi z\right)\right]-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right]\left \lceil{\mathrm{Re}\left(\frac{z}{2}\right)}-\frac{3}{4}\right \rceil 2 \pi i } $$Where $\lceil . \rceil$ is the cieling function. Of course there is probably a nicer form. Feel free to derive it
To account for peridicity of dislocations in $x$ direction, the expression $\log\left(z\right)$ in Eqs(2-7) of the paper, it should be replaced by:
$$\lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-na\right) -\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i -2\log\left(\frac{m\,\,\!!}{\sqrt{\pi}} \right),$$which has the closed form: $$ \boxed{ \log\left[\sin\left(\pi\frac{z}{a}\right)\right]-\mathrm{Sign}\left[\mathrm{Im}\left(\frac{z}{a}\right)\right]\left \lceil{\mathrm{Re}\left(\frac{z}{2a}\right)}-\frac{3}{4}\right \rceil 2 \pi i. } $$
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import mapp4py
from mapp4py import md
from lib.elasticity import rot, cubic, resize, displace, HirthEdge, HirthScrew
In [ ]:
from mapp4py import mpi
if mpi().rank!=0:
with open(os.devnull, 'w') as f:
sys.stdout = f;
md.export_cfg has a call method that we can use to create quick snapshots of our simulation box
In [ ]:
xprt = md.export_cfg("");
In [ ]:
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
nlyrs_fxd=2
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
Since mapp4py.md.atoms.cell_chenge() only accepts integer values start by creating a $\langle110\rangle\times\langle112\rangle\times\langle111\rangle$ cell
In [ ]:
sim.cell_change([[1,-1,0],[1,1,-2],[1,1,1]])
Now one needs to cut the cell in half in $[111]$ direction. We can achive this in three steps:
In [ ]:
H=np.array(sim.H);
def _(x):
if x[2] > 0.5*H[2, 2] - 1.0e-8:
return False;
else:
x[2]*=2.0;
sim.do(_);
_ = np.full((3,3), 0.0)
_[2, 2] = - 0.5
sim.strain(_)
In [ ]:
displace(sim,np.array([sim.H[0][0]/6.0,sim.H[1][1]/6.0,0.0]))
In [ ]:
max_natms=100000
H=np.array(sim.H);
n_per_area=sim.natms/(H[0,0] * H[1,1]);
_ =np.sqrt(max_natms/n_per_area);
N0 = np.array([
np.around(_ / sim.H[0][0]),
np.around(_ / sim.H[1][1]),
1], dtype=np.int32)
sim *= N0;
In [ ]:
H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[1][1] += 50.0
resize(sim, H_new, np.full((3),0.5) @ H)
In [ ]:
C_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
Q=np.array([np.cross(s,b)/np.linalg.norm(np.cross(s,b)),s/np.linalg.norm(s),b/np.linalg.norm(b)])
hirth = HirthScrew(rot(C_Fe,Q), rot(b*0.5*a,Q))
In [ ]:
ctr = np.full((3),0.5) @ H_new;
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
def _(x,x_d,x_dof):
sy=(x[1]-ctr[1])/H[1, 1];
x0=(x-ctr)/H[0, 0];
if sy>s_fxd or sy<=-s_fxd:
x_dof[1]=x_dof[2]=False;
x+=b_norm*hirth.ave_disp(x0)
else:
x+=b_norm*hirth.disp(x0)
sim.do(_)
In [ ]:
H = np.array(sim.H);
H_inv = np.array(sim.B);
H_new = np.array(sim.H);
H_new[0,0]=np.sqrt(H[0,0]**2+(0.5*b_norm)**2)
H_new[2,0]=H[2,2]*0.5*b_norm/H_new[0,0]
H_new[2,2]=np.sqrt(H[2,2]**2-H_new[2,0]**2)
F = np.transpose(H_inv @ H_new);
sim.strain(F - np.identity(3))
In [ ]:
xprt(sim, "dumps/screw.cfg")
In [ ]:
def make_scrw(nlyrs_fxd,nlyrs_vel,vel):
#this is for 0K
#c_Fe=cubic(1.5187249951755375,0.9053185628093443,0.7249256807942608);
#this is for 300K
c_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
#N0=np.array([80,46,5],dtype=np.int32)
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
Q=np.array([np.cross(s,b)/np.linalg.norm(np.cross(s,b)),s/np.linalg.norm(s),b/np.linalg.norm(b)])
c0=rot(c_Fe,Q)
hirth = HirthScrew(rot(c_Fe,Q),np.dot(Q,b)*0.5*a)
sim.cell_change([[1,-1,0],[1,1,-2],[1,1,1]])
displace(sim,np.array([sim.H[0][0]/6.0,sim.H[1][1]/6.0,0.0]))
max_natms=1000000
n_per_vol=sim.natms/sim.vol;
_=np.power(max_natms/n_per_vol,1.0/3.0);
N1=np.full((3),0,dtype=np.int32);
for i in range(0,3):
N1[i]=int(np.around(_/sim.H[i][i]));
N0=np.array([N1[0],N1[1],1],dtype=np.int32);
sim*=N0;
sim.kB=8.617330350e-5
sim.create_temp(300.0,8569643);
H=np.array(sim.H);
H_new=np.array(sim.H);
H_new[1][1]+=50.0
resize(sim, H_new, np.full((3),0.5) @ H)
ctr=np.dot(np.full((3),0.5),H_new);
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
s_vel=0.5-0.5*float(nlyrs_vel)/float(N0[1])
def _(x,x_d,x_dof):
sy=(x[1]-ctr[1])/H[1][1];
x0=(x-ctr)/H[0][0];
if sy>s_fxd or sy<=-s_fxd:
x_d[1]=0.0;
x_dof[1]=x_dof[2]=False;
x+=b_norm*hirth.ave_disp(x0)
else:
x+=b_norm*hirth.disp(x0)
if sy<=-s_vel or sy>s_vel:
x_d[2]=2.0*sy*vel;
sim.do(_)
H = np.array(sim.H);
H_inv = np.array(sim.B);
H_new = np.array(sim.H);
H_new[0,0]=np.sqrt(H[0,0]**2+(0.5*b_norm)**2)
H_new[2,0]=H[2,2]*0.5*b_norm/H_new[0,0]
H_new[2,2]=np.sqrt(H[2,2]**2-H_new[2,0]**2)
F = np.transpose(H_inv @ H_new);
sim.strain(F - np.identity(3))
return N1[2],sim;
In [ ]:
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
nlyrs_fxd=2
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
In [ ]:
sim.cell_change([[1,1,1],[1,-1,0],[1,1,-2]])
H=np.array(sim.H);
def _(x):
if x[0] > 0.5*H[0, 0] - 1.0e-8:
return False;
else:
x[0]*=2.0;
sim.do(_);
_ = np.full((3,3), 0.0)
_[0,0] = - 0.5
sim.strain(_)
displace(sim,np.array([0.0,sim.H[1][1]/4.0,0.0]))
In [ ]:
max_natms=100000
H=np.array(sim.H);
n_per_area=sim.natms/(H[0, 0] * H[1, 1]);
_ =np.sqrt(max_natms/n_per_area);
N0 = np.array([
np.around(_ / sim.H[0, 0]),
np.around(_ / sim.H[1, 1]),
1], dtype=np.int32)
sim *= N0;
In [ ]:
# remove one layer along ... direction
H=np.array(sim.H);
frac=H[0,0] /N0[0]
def _(x):
if x[0] < H[0, 0] /N0[0] and x[1] >0.5*H[1, 1]:
return False;
sim.do(_)
In [ ]:
H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[1][1] += 50.0
resize(sim, H_new, np.full((3),0.5) @ H)
In [ ]:
C_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
_ = np.cross(b,s)
Q = np.array([b/np.linalg.norm(b), s/np.linalg.norm(s), _/np.linalg.norm(_)])
hirth = HirthEdge(rot(C_Fe,Q), rot(b*0.5*a,Q))
In [ ]:
_ = (1.0+0.5*(N0[0]-1.0))/N0[0];
ctr = np.array([_,0.5,0.5]) @ H_new;
frac = H[0][0]/N0[0]
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
def _(x,x_d,x_dof):
sy=(x[1]-ctr[1])/H[1, 1];
x0=(x-ctr);
if(x0[1]>0.0):
x0/=(H[0, 0]-frac)
else:
x0/= H[0, 0]
if sy>s_fxd or sy<=-s_fxd:
x+=b_norm*hirth.ave_disp(x0);
x_dof[0]=x_dof[1]=False;
else:
x+=b_norm*hirth.disp(x0);
x[0]-=0.25*b_norm;
sim.do(_)
In [ ]:
H = np.array(sim.H)
H_new = np.array(sim.H);
H_new[0, 0] -= 0.5*b_norm;
resize(sim, H_new, np.full((3),0.5) @ H)
In [ ]:
xprt(sim, "dumps/edge.cfg")
In [ ]:
def make_edge(nlyrs_fxd,nlyrs_vel,vel):
#this is for 0K
#c_Fe=cubic(1.5187249951755375,0.9053185628093443,0.7249256807942608);
#this is for 300K
c_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
#N0=np.array([80,46,5],dtype=np.int32)
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
# create rotation matrix
_ = np.cross(b,s)
Q=np.array([b/np.linalg.norm(b), s/np.linalg.norm(s), _/np.linalg.norm(_)])
hirth = HirthEdge(rot(c_Fe,Q),np.dot(Q,b)*0.5*a)
# create a unit cell
sim.cell_change([[1,1,1],[1,-1,0],[1,1,-2]])
H=np.array(sim.H);
def f0(x):
if x[0]>0.5*H[0][0]-1.0e-8:
return False;
else:
x[0]*=2.0;
sim.do(f0);
_ = np.full((3,3), 0.0)
_[0,0] = - 0.5
sim.strain(_)
displace(sim,np.array([0.0,sim.H[1][1]/4.0,0.0]))
max_natms=1000000
n_per_vol=sim.natms/sim.vol;
_=np.power(max_natms/n_per_vol,1.0/3.0);
N1=np.full((3),0,dtype=np.int32);
for i in range(0,3):
N1[i]=int(np.around(_/sim.H[i][i]));
N0=np.array([N1[0],N1[1],1],dtype=np.int32);
N0[0]+=1;
sim*=N0;
# remove one layer along ... direction
H=np.array(sim.H);
frac=H[0][0]/N0[0]
def _(x):
if x[0] < H[0][0]/N0[0] and x[1]>0.5*H[1][1]:
return False;
sim.do(_)
sim.kB=8.617330350e-5
sim.create_temp(300.0,8569643);
H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[1][1] += 50.0
ctr=np.dot(np.full((3),0.5),H);
resize(sim,H_new, np.full((3),0.5) @ H)
l=(1.0+0.5*(N0[0]-1.0))/N0[0];
ctr=np.dot(np.array([l,0.5,0.5]),H_new);
frac=H[0][0]/N0[0]
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
s_vel=0.5-0.5*float(nlyrs_vel)/float(N0[1])
def f(x,x_d,x_dof):
sy=(x[1]-ctr[1])/H[1][1];
x0=(x-ctr);
if(x0[1]>0.0):
x0/=(H[0][0]-frac)
else:
x0/= H[0][0]
if sy>s_fxd or sy<=-s_fxd:
x_d[1]=0.0;
x_dof[0]=x_dof[1]=False;
x+=b_norm*hirth.ave_disp(x0);
else:
x+=b_norm*hirth.disp(x0);
if sy<=-s_vel or sy>s_vel:
x_d[0]=2.0*sy*vel;
x[0]-=0.25*b_norm;
sim.do(f)
H = np.array(sim.H)
H_new = np.array(sim.H);
H_new[0, 0] -= 0.5*b_norm;
resize(sim, H_new, np.full((3),0.5) @ H)
return N1[2], sim;
In [ ]:
nlyrs_fxd=2
nlyrs_vel=7;
vel=-0.004;
N,sim=make_edge(nlyrs_fxd,nlyrs_vel,vel)
In [ ]:
xprt(sim, "dumps/edge.cfg")
In [ ]:
_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;
In [ ]:
C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)
In [ ]:
B = np.linalg.inv(
np.array([
[C[0, 0, 0, 0], C[0, 0, 1, 1], C[0, 0, 0, 1]],
[C[0, 0, 1, 1], C[1, 1, 1, 1], C[1, 1, 0, 1]],
[C[0, 0, 0, 1], C[1, 1, 0, 1], C[0, 1, 0, 1]]
]
))
In [ ]:
_ = np.roots([B[0, 0], -2.0*B[0, 2],2.0*B[0, 1]+B[2, 2], -2.0*B[1, 2], B[1, 1]])
mu = np.array([_[0],0.0]);
if np.absolute(np.conjugate(mu[0]) - _[1]) > 1.0e-12:
mu[1] = _[1];
else:
mu[1] = _[2]
alpha = np.real(mu);
beta = np.imag(mu);
p = B[0,0] * mu**2 - B[0,2] * mu + B[0, 1]
q = B[0,1] * mu - B[0, 2] + B[1, 1]/ mu
K = np.stack([p, q]) * np.array(mu[1], mu[0]) /(mu[1] - mu[0])
K_r = np.real(K)
K_i = np.imag(K)
In [ ]:
Tr = np.stack([
np.array(np.array([[1.0, alpha[0]], [0.0, beta[0]]])),
np.array([[1.0, alpha[1]], [0.0, beta[1]]])
], axis=1)
def u_f0(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) + x[0])
def u_f1(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) - x[0]) * np.sign(x[1])
def disp(x):
_ = Tr @ x
return K_r @ u_f0(_) + K_i @ u_f1(_)
In [ ]:
_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;
C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)
disp = crack(C)
In [ ]:
n = 300;
r = 10;
disp_scale = 0.3;
n0 = int(np.round(n/ (1 +np.pi), ))
n1 = n - n0
xs = np.concatenate((
np.stack([np.linspace(0, -r , n0), np.full((n0,), -1.e-8)]),
r * np.stack([np.cos(np.linspace(-np.pi, np.pi , n1)),np.sin(np.linspace(-np.pi, np.pi , n1))]),
np.stack([np.linspace(-r, 0 , n0), np.full((n0,), 1.e-8)]),
), axis =1)
xs_def = xs + disp_scale * disp(xs)
fig, ax = plt.subplots(figsize=(10.5,5), ncols = 2)
ax[0].plot(xs[0], xs[1], "b-", label="non-deformed");
ax[1].plot(xs_def[0], xs_def[1], "r-.", label="deformed");