Introduction

The purpose of this document is to give a brief description of DDLab. If advanced

users are interested in the algorithms of the code, they may simply look at the
source code. We have tried to write and organize the subroutines in a compact yet
readable manner.
Warning: The main executable file dd3d.m requires the compilation of
mex -O SegSegForces.c

You need to have a C compiler installed on your computer to make it work. This is
ok when you use gcc-3.3. You may have difficulty compiling with higher version of
gcc.

Input parameters of DDLab

There are various input parameters in DDLab. In this section, we don't give a
detailed expla- nation of each parameter but we explain some. The relation
between parameters will be explained in more detail in the following sections.
We can divide all input parameters into the following categories.

Dislocation structure : rn, links* Mobility : mobility

Integrator : integrator
Topological changes : lmax, lmin, areamin, areamax, rmax, rann, rntol,
doremesh, docollision, desepraration
Time controls : totalstpes, dt0
Display controls : plotfreq, plim, viewangle, printfreq, printnode
Materials constant : MU (shear modulus), NU (Poisson's ratio), Ec (core
energy), a (core radius)
Miscellaneous : appliedstress, maxconnections

Nodal representation of dislocation structure

Figure 1: Nodal representation of dislocation structure Fig. 1 shows a simple

approach that can represent an arbitrary dislocation network. The dislocations are
specified by a set of nodes that are connected with each other by straight
segments. Each segment has a nonzero Burgers vector. Because the Burgers vector
is defined only after a sense of direction is chosen for the dislocation line, we can
define bij

Mobility
In DDLab, the motion of dislocation is expressed by the motion of nodes. Thus,
mobility gives the velocity. DDLab contains several mobility laws; mobfcc1.m,
mobfcc0.m, and mobbcc0.m.

Integrator
An integrator scheme is used to calculate the position of each nodes for a given
time step. The nodal forces and mobility model comprise the equation motion for
the nodes
where g_i sums up both the nodal force and the
mobility model. This is a first ordinary differential equation and we solve it by
integration. DDLab contains two numerical integrator: eulerforward.m,
newtonkrylov.m. ode15s.m, trapezoid.m. In this work, we use trapezoid.m is
used.

Topological Changes
In order to accurately represent dislocation behavir and also for numerical
efficiency, we need a way to handle topological changes, i.e. changes in the
connectivity of nodes. For example, should a line change its length and/or
curvature during the simulation, it may become necessary to adjust the number of
nodes used to represent this dislocation line. The following parameters are used
for topological changes in DDLAB.

lmax: Maximum length of a dislocation segment

lmin: Minimum length of a dislocation segment

areamin: Maximum area criterion

areamax: Minimum area criterion
rmax: maximum distance a node may travel in one cycle
rann: annihilation distance used to calculate the collision of dislocation
lines
rntol: solution tolerance used to control the automatic timestepping
doremesh: Toggle remesh.m for remeshing nodes. doremesh is set as zero to
disable

the use of remesh rule. When it is set as zero, it is enabled.

Time controls
The following parameters are for time controls.

totalsteps : number of cycles that are run for completion of dd3d command
dt0 : largest timestep that can be taken during a cycle

Display controls
The following parameters are for display controls

plotfreq : number of cycles between monitored node write statements

plim : limits of plotting space
viewangle : angle of viewpoint for the 3D plot of the geometry
printfreq : number of cycles between monitored node write statements
printnode : nodeid of the monitored node

Materials constants
The following parameters are Materials constants used in DDLAB.

Mu: Shear modulus

NU: Poisson's ratio
Ec: Dislocation core energy

Dislocation core radius used for non-singular force calculation

Miscellaneous

maxconnections : maximum number of segments a node may have

Flow of the code in DDLab

After setting up the initial configuration and other parameters, both DDLab and
ParaDiS will follow the same algorithm to simulate dislocation dynamics, as is
listed below:

calculate the force of each node

calculate the velocity of each node
calculate the new position of each node
make necessary topological changes
repeat (1) to (4) until the maximum step is reached

DDLab just track all nodal positions and plot them. However, to better understand
the code structure, two flow charts are presented in this primer: subroutine map
(Fig. 2) and function map (Fig. 3). The subroutine map consists of the MATLAB
subroutine file name, but the function map consists of the function of each
subroutine file. Thus, if you want to know the code in detail, just compare the
subroutine map with the function map, then find the file you want and take a look
at. Thus, these two maps are very useful to study the code structure of DDLab. All
functions of subroutines are explained in the following section which includes
some physics for the detailed explanations.

Subroutines details

SegSegforce.c

This subroutine calculates the force on the two end nodes between two dislocation
seg- ments (so, four nodes) in Fig. 4. The sum of these four force vectors should be
zero by action-reaction principle. f1 + f2 + f3 + f4 = 0
Figure 4: The interaction between the two dislocation segments

cleanupnodes.m

This subroutine start by removing all of the bad nodes in the system and cleans up
the data structures to include good nodes only. A bad node is a point without
connection to any dislocation segment. The point is usually generated by the
shrinkage of a dislocation loop as shown in Fig. 5. At the final stage, the loop
consists of two nodes, so the line direction of two segments is opposite, but their
Burgers vector is the same. Thus, two segments should annihilate. Then, only two
nodes remains. These two nodes are only a mathematical points in elastic medium.
In real material, the loop disappears after shrinkage, so the two nodes in DDLab
has to be removed. Thus, cleanupnodes.m performs the removal of these kinds of
isolated nodes.
Note:
This is very similar with to another case (two nodes and two segments) of
collision.m. Sometimes segments are not permitted to annihilate when the sum
of Burgers vectors are not zero. However, in the case of loop shrinkage, two
segments always have the opposite direction of line sense with the same Burgers
vector. Thus, the annihilation occurs. See collision.m section.

genconnectivity.m

This subroutine generates connectivity and linksinconnect. In fact, rn and links

matrices are enough to describe the dislocation structure. These two matrices give
linkid when we know nodeid. Conversely, connectivity and linksinconnect give
nodeid when

Figure 5 The shrinkage of a dislocation loop. we know linkid Here, nodeid is the
node number which is the row number of rn matrix. linkid is the link number which
is the row number of links matrix.

Integrator
Integrator calls drndt.m, drndt.m computes velocity and force on each node, and
force on each segment. Then, by using these values, Integrator integrates the
velocity with respect to time step. Then, we can get the nodal positions at the

current time step. The nodal forces and the mobility model comprise the equations
of motion of nodes:

where function gi sums up both the nodal force and the mobility model. This is first
ordinary differential equation (ODE), for which the simplest numerical integrator is
the so- called Euler forward method:
ri(t + t) = ri + gi(rj(t))t:
A numerical integrator with a higher order of accuracy is the trapezoid method, ri(t
+ t) = gi(rj(t)) + gi(rj(t + t))
Another important feature of integrator is the time step adjustment. When the
time step is large, we can perform simulation which require long physical time.
However, the large time step always bring large error. Conversely, when the time
step is too short, we can perform the more precise simulation, but it take too long
CPU time to get the result. Thus, it is better to adjust the time step in every
cycle. The present version of DDLab has the time step adjustment as below. This is
the part of the integrator code. ...err = rnvec1 rnvec0 (vnvec + vnvec0) / 2 *
dt;...errmag = max(abs(err));...if(errmag < rntol)break;elsedt = dt / 2;...maxchange = 1.2;... If
the defined err is larger than rntol, the time step is adjusted as the half of current
time step. If the defined err is smaller than rntol, the time step is adjusted as 1.2
times of current time step. Thus, DDLab changes its time step dynamically,
enabling a much more efficient simulation.
A series of time integration schemes are also included with the distribution named
int eulerforward, int eulerbackward, int newtonkrylov. int eulerforward is
an explicit forward integration scheme. <tt>int eulerbackward is a simple
implicit time integration scheme, and int newtonkylov is a more sophisticated
inexact implicit time integration scheme based on the Newton-Krylov matrix-free
solution method. The core of the int newtonkrylov employs a method written by
C. T. Kelley. The time integration scheme is specified in the input deck with the
integrator input parameter.

4.0.14 drndt.m drndt.m calls segforcevec.m. segforcevec.m calculates the force

on each segment, then return it to drndt.m. Then, drndt.m calls mobility. mobility
calculates the force and velocity on each node and return them to drndt.m Finally,
drndt.m return the velocity on node to integrator, which integrate the velocity
with respect to time step. Then, we can get the nodal positions at the current
time step. 4.0.15 segforcevec.m segforcevec calculates the force on each
segment. Note that this is the force not on the node but on the segment. It is
because DDLAB use Peach-Koehler formula fPK(x) = ((x)b)(x) . When we compute PK
force, we need to know the line sense vector () for segment. Thus, firstly DDLAB
calculates the force on each segment which has two end nodes. Then, mobility
function gives the half of force to each node. Finally, we can get the nodal force.
segforcevec.m includes the four functions; constructsegmentlist, pkforcevec,
selfforcevec, remoteforcevec. constructsegmentlist give us the list of all
information about segments as the format n0, n1, bx, by, bz, x0, y0, z0, x1, y1, z1,
nx, ny, nz, where n0 and n1 are nodeids of two end nodes. bx, by, bz is the Burgers
vector of one segment which consists of two end node, n0 and n1. x0, y0, z0 is the
position of node n0 and x1, y1, z1 is the position of node n1. nx, ny, nz is the
normal vector of the glide plane on which the dislocation segment lies. Based on
this list, segforcevec calculate the forces on each segment. selfforcevec calculates
the force on each segment by itself. remoteforcevec calculates the force on the
segment on each segment by the other segments. Finally, the sum of these forces
is the force on segment as below. fseg=[fpk, fpk]*0.5+[fs0, fs1]+[fr0, fr1]; The first
component of fseg (fpk, fs0, fr0) are the force vectors on node n0 and the second
ones are the force vectors on the node n1. The node n0 and n1 are the end nodes
of each segment as shown in Fig. 5.
Figure 7: (a) The example dislocation structure. (b) DDLAB calcuate the force on
the segment by segforcevec.m (c) Then, the half of the force is assigned to two
end nodes of each segment (d) The force on the node has several arms is just the
summation of forces on corresponding arms. 4.0.16 mobfcc1.m, mobfcc0.m,
mobfcc1.m A dislocation's response to forces may also depends on the line
orientation and varies significantly from one material to another. Similar to the
core energy, how dislocation respond to forces is a function of the non-linear
atomistic interactions in the dislocation core and, as such, is beyond the scope of
linear elasticity theory. The response functions have to be imported into the line
DD model in the form of external materials input, either from atomistic
simulations or from experiments A series of mobility functions are included in the
distribution named mobbcc0 mobbcc0b and mobfcc1. mobbcc0 and mobbcc0b are
intended to simulate bcc behavior in which the the screw dislocations are able to

glide in any plane normal to their Burgers vector. The difference between them is
slight and appears in the mobility of dislocations of mixed character. mobfcc1 is
intended to simulate fcc behavior in which the screw dislocation are not able to
cross-slip. The mobility function is chosen in the input deck with the mobility input
parameter.

4.0.17 separation.m For numerical and physical reasons, line-DD simulations need
to handle topological changes, i.e. changes on the connectivity between nodes
since we may want to adjust the number of nodes that represent a dislocation line
if the line gets longer or shorter during the simulation, or when two dislocation
lines meet in space, they may either annihilate or zip together to form a junction,
which also results in a change of nodal topology. Thus many types of topological
changes can be encountered in a line-DD simulation. Fortunately, since we use a
nodal representation here, all topological changes can be implemented through
two basic operators: split (one node split into two ) and merge (two nodes merge
into one) (See next section) and . The implementation of these two operators is
straightforward - all one needs to do is to make sure that at the end of the
operation the Burgers vector sum rule at every node and segment is still satisfied,
moreover, two nodes are either disconnected or connected only once, and each
node is connected with at least two other nodes, if a node has no segment, it will
be deleted. For example, for a 4-arm node such as P in there are 3 different ways
to partition its arms: (12)(34), (13)(24) and (14)(23). It is reasonable to expect that
the way nature would choose should be the one that gives rise to the maximum
energy dissipation rate, which is defined below. Suppose an n-arm mode i stays
intact (not splitting) and it feels a force fi and will move at velocity vi . Then the
local energy dissipation rate is, Wi = fivi. Now suppose that node i splits into two
nodes P and Q, such that node P retains 1, ... ,s of the original neighbors, and
node Q retains the remaining neighbors. Let fP and fQ be the forces on the two
nodes and vP and vQ be their velocities given by the mobility function. Then the
local energy dissipation rate is, WPQ = fPvP + fQvQ.
If WPQ > Wi, then node i prefers to split into two nodes P and Q instead
of moving

as a single node. The energy dissipation rate can be computed for all possible
(topological distince) modes to split i. The mode with the highest energy

dissipation rate is preferred. If a node will split in next step, the two new nodes
actually stay at the same location as the parent node at the current step. Because
the velocities of the two nodes is different (otherwise the node should not split),
the two nodes will be move away from each other in the next time step.

collision.m

remesh.m

Self consistency in DDLab units

The unit system of DDLab input file should be self-consistent. We recommend to
input all variables in S.I. units.
For instance, stresses in Pascal, distances in meter, times in seconds and velocities
in m/s.
The mobility law is such that
with b in meter, in Pascal, and v in m/s. That makes the mobility parameter M in
Pa 1s 1.
We often use the drag coefficient B = M 1 in the mobility law. The unit of B is

As another option, we can change the unit for distance (including both nodal
positions and Burgers vectors) to nm, whereas the unit for stress is still Pa, the
unit of time is still second, and the unit for mobility is still Pa 1s 1.

Introduction
Frank-Read source is a type of dislocation multiplication mechanism. Consider a
segment whose ends are pinned (corresponding to nodes in a network,
precipitates, or sites where the dislocation leaves the glide plane). Under a certain
applied stress the segment bows out by glide. As bow-out proceeds, the radius of
curvature of the line decreases and the line-tension forces tending to restore the
line to a straight configuration increase. For stress less than a critical value, a
metastable equilibrium configuration is attained, in which the line-tension force

balances that caused by the applied stress. For the large bow-out case, following
equilibrium condition holds:

(Hirth-Lothe, p. 752)
where r is the radius of the loop. The radius of the curvature is a minimum when r
= L / 2. Hence the maximum stress for which local equilibrium is possible is given
by the equation above with r = L / 2. For the typical case that L = 103 and = 0.33,
the critical stress for a dislocation initially pure edge and pure screw, respectively,
is * = 0.5b / L and * = 1.5b / L.
When the net local resolved shear stress( the applied stress plus the internal
stresses) exceeds * , the loop has no stable equilibrium configuration but passes
through the successive positions. Provided that the expanding loop neither jogs out
of the original glide plane because of intersections with other dislocations nor is
obstructed from rotating about the pinning point, it will annihilate over a portion
of its length, creating a complete closed loop and restoring the original
configuration. A sequence of loops then continues to form from the source until
sufficient internal stresses are generated for the net resolved shear stress at the
source to drop below * .
To compare the analytic result with DDLab simulation, we need to use the relation
= rc / 2. According to the non-singular continuum theory of dislocation (JMPS 54,
561-587, 2006), the energy of a circular prismatic dislocation loop is

According to Hirth and Lothe, the energy of the same loop is,

Therefore

Reference: J. P. Hirth and J. Lothe, Theory of Dislocations", 2nd ed. (Wiley, New
York, 1982)

DDLab Input Parameter

DDLab has included a sample Frank-Read Source input file in its latest version
(2008 January) at ~/Inputs/input_frank_read.m. Following is a brief introduction
to some parameters important to Frank-Read source test.
1. rn--it gives the initial coordinations of nodes. The default setting is a straight
dislocation line pinned (tag=7) at (1200,1200,1200) and (-1200,-1200,-1200) with a
mobile (tag=0) node in the middle, (0,0,0). Of course, 1200 can be substituted into
arbitrary numbers when we investigate the influence of initial segmental length.
2. links--a data structure which gives the information of dislocation segment
connectivity, burgers vectors( [0.5 0.5 0.5] default ) and glide planes (-1,1,0)
default ).
3. totalsteps--number of cycles that are run for completion of dd3d command. It
should be large enough to secure dislocation configuration has arrived at its
equilibrium configuration, for example, 1000 or above.
4. appliedstress--external applied stress in the form of 3*3 symmetric tensor.
Peach-Koehler formula will be needed to give the exact force on each segment,
i.e.

We need to use this formula to give F-R * under certain appliedstress.

5. a--dislocation core radius used for non-singular force calculation, sometimes
written as rc.
6. plotfreq--number of cycles between monitered node write statements, higher
number makes observation simpler. When the dislocation line goes beyond plotting
limits, simply Increasing plim can provide a better view.
7. lmin and lmax -- minimum and maximum length of a dislocation segment for
remeshing. The shorter segments are, the smoother dislocation line becomes.

8. mobility --mobility law based on which dislocation responds to forces (mobbcc0

for BCC, which is the default setting; mobfcc1 for FCC). For FCC crystals, glide
plane is uniquely defined rather BCC, i.e. dislocation motion is confined to a
preferred plane.
We can manipulate these parameters to investigate Frank-Read Sources via DDLab.

F-R Test
1. Criteria to define *
In the case of Frank-Read source, It is difficult to give a quantitative definition on
critical applied stress. When totalsteps is large enough, the bowing process will
eventually stop unless that critical stress has been reached. Let's define a certain
geometric configuration beyond which the dislocation line becomes unstable and
keeps growing, even forming loops. When such configuration is formed, we denote
the current external applied stress as the critical one. It is convenient to use
maximum distance criteria, that we can extract the coordinates of all nodes from
rn and find out the maximum distance from the center of initial dislocation line. If
the maximum distance exceeds an empirical value, let's say, L, which is the length
of initial dislocation line, we can regard the configuration as "critical" and
determine the critical stress, * . Several trials have proven that L is a safe one. As
long as totalsteps is large enough, 2L, 4L etc will also work well.
2. Basic settings
The initial dislocation line is pinned at (-l,-l,-l) ,(l, l, l). (0, 0, 0) is a mobile node
at the center.
We can change any parameter before running dd3d.m. Simple test cases have
shown that when l ranges from 1,000 to 10,000, the external applied stress is at
the order of 10 4.

rn

= [ 1000
1000
1000
7;
-1000
-1000 -1000
7;
0
0
0
0];
links = [1 3 0.5 0.5 0.5 -1 1 0;

%pinned at two ends

%mobile node in the middle

3 2 0.5 0.5 0.5 -1 1 0];

glide plane (-1 1 0)
MU = 1;
NU = 0.305;

%burgers vector [0.5 0.5 0.5],

%elasticity constants for isotropic media

maxconnections=8;
lmax = 1000;
lmin = 200;
% fixed lmax and lmin for all cases
areamin=lmin*lmin*sin(60/180*pi)*0.5;
areamax=20*areamin;
a=lmin/sqrt(3)*0.5;
%defined in terms of lmin
Ec = MU/(4*pi)*log(a/0.1);
%Ec=0 when compared with analytical
expression
totalsteps=200;
%enough for rough estimations
dt0=1e7;
mobility='mobfcc1';
%FCC mobility law
%Drag (Mobility) parameters
Bscrew=1e0;
Bedge=1e0;
Beclimb=1e2;
Bline=1.0e-4*min(Bscrew,Bedge);
integrator='int_trapezoid';
rann = 0.5*a;
rntol = 0.5*rann;
doremesh=1;
docollision=1;
doseparation=1;
elasticinteraction=1;
%line tension model when this is zero
plotfreq=1;
plim=10000;
appliedstress =1e-3.*[2 0 1; 0 2 -1; 1 -1 0];
% a value times a stress
tensor
viewangle=[45 -45];
printfreq=1;
%set as 10 for better view
printnode=3;
rmax=100;

3. Algorithm
The critical stress is determined digit by digit. Assuming the critical externally
applied stress is 2.35e-4 times a matrix A, to determine the first digit, simulation

is run for 400 steps when the applied stress is 1e-4*A. The stress will be increased
by 1e-4*A unless the maximum criteria is not reached. We will see, the first digit is
3. As for the second digit, the initial applied stress starts from 2.1e-4*A and
simulation is run for 1000 steps. In the similar manner, 2.2e-4*A will be applied if
critical criteria is not met. Finally, the second digit is 4. To determine the third
digit, simulation starts at 2.31e-4*A...Of course, the matrix A and totalsteps are
subject to change. More simulation steps will be needed when the stress increase
is smaller, i.e. better accuracy. For detailed coding, please refer to crstrfr.m.
4. Investigation
4.1. L dependence
Based on the maximum distance criteria, we can easily determine the critical
applied stress for different initial straight line lengths. Here is a graph showing the
relationship between c and L at a=50.

According to the formula, there is a linear relation between the critical stress and
segment length L, which was also proven by numerical results.

The relative error decreases from 0.26 to 0.20 when the initial segment coordinate
parameter l goes from 1000 to 10,000; however, the error can't be negligible.
Better convergence should occur at a smaller a value. Numerical tests shows that,
when a=1, relative error decreases from 0.20 to 0.18 when l goes from 1000 to
6000. Lower error is expected for longer segments.

Longer simulation time can help eliminate the intersections at large l.

A interesting phenomenon is that numerical result for a/5 well fits analytical result
for a. Relative error goes to 0.01 when l goes to 10,000.

It shows that analytical result will be more accurate at smaller a/L ratio.
4.2. a dependence
When a varies from 1 to 50 and l is fixed at 5000, numerical results are given,
shown in following figures

Same as the formula, the critical stress has a linear dependence on ln(L/a)/L

We can see the decrease in relative error.

When the length of initial segment increases, relative error seems to decrease
when a approaches 0.

4.3. lmax influence

Above results are based on the setting that lmax=l, lmin=l/10 other than the
default lmax=1000, limin=200. Though the dislocation line is less smoother,
numerical results are the same for the case l=5000. For computation simplicity, we
can set lmax and lmin proportional to l and get reasonable results.

Conclusion
After conducting several test cases of Frank-Read Source using DDLab, we have
shown:
1, The critical stress is linearly dependent on ln(L/a)/L;
2. Numerical results are closer to theoretical ones at small a and large L.
3. Numerical critical stress for a/5 well fits theoretical results for a

Some of the plots above will be updated when more accurate data are available
for a=1 (or less) and L>6000.

Introduction
This tutorial will:

Study the formation and dissociation of Lomer-Cottrell (binary) junctions in

fcc crystals using DDLab
Introduce running a batch of simulations, and calculating junction length
and critical resolved shear stress to cause dissociation of the junction
Determine the dependence of junction length and critical stress at
dissociation on the initial orientation of the intersecting dislocations
Determine the dependence of junction length and critical stress at
dissociation on the initial length of the dislocations
Discuss any limitations. Physically the junction should be sessile. But in the
simulation, the inner nodes of the junction segment are able to move. The
end nodes of the junction are confined to the intersection line as they
should.

This tutorial assumes that the reader is familiar with the basics of DDLAB. New
users should read DDLab_Manual_01 before proceeding.
A copy of the input decks used in this tutorial can be downloaded from the link at
the bottom of this page.
A tutorial for ParaDiS will be created in the future.

Background
Plastic deformation of metals is usually governed by the motion of dislocations and
their reactions, such as nucleation, pinning, and multiplication. One particularly
important interaction occurs when two dislocations on different glide planes
approach each other. These dislocations could either repel or attract each other
depending on their Burgers vectors and orientations. If they are pulled together,
the dislocations will combine in a way that minimizes their energy. This often
results in the "zipping" of the dislocations along the line of intersection of the glide
planes to form a junction known as a Lomer-Cottrell lock. These junctions have
been shown to be a governing cause of strain hardening in fcc crystals.

Simulation Design
Initial Orientation

Figure 1: Formation of a Lomer-Cottrell lock

These simulations are performed using the material constants for gold in order to
conform with past computational and experimental results obtained by Madec et.
al<ref name="Madec">From Dislocation Junctions to Forest Hardening, Phys. Rev.
Lett. 89 255508 (2002), R. Madec, B. Devincre, and L.P. Kubin</ref>, although any
fcc metal could be used. In order to simplify the simulations, two dislocations with
identical lengths, but on differing slip planes, are initially placed such that they
intersect at their center. The first dislocation glides on the
plane henceforth referred to as dislocation 1. On the other hand, dislocation 2 glides on
the

plane. Both of these slip planes are part of the

in fcc crystals. The dislocations have Burgers vectors of

respectively, regardless of their line direction

slip system
and

. A Lomer-Cottrell junction could

form at the intersection of these two planes, i.e. along the

direction, if
the conditions are correct, as shown in Figure 1. Thus, it is desired to determine
the junction behavior at all combinations of line directions.
If we consider all possible orientations of a straight dislocation of constant length
and we take the origin to be the initial point of contact of the dislocations, we can
define an end point (i.e. fixed node),
, of the dislocation by a parametric
equation of :

Figure 2: Discretization of dislocation orientations

where

is half the length of the dislocation;

is a unit vector in the direction of

the line of intersection,

; and is the unit normal of the glide plane. The
vectors and
are then the end points of the dislocation. It is obvious from
the equation that varying from
generates a circle on the glide plane.
We have chosen to simulate 21 evenly distributed line orientations for both
dislocations in this study. This results in a discretization for both dislocations
shown in Figure 2, where the points represent one fixed node in each orientation
and the initial line direction is taken as the vector from the origin to that point.
LC_make_dis.m calculates the position of an end node of each dislocation in this
manner using the following code:
% Make dislocations
t=[-1:0.1:1]*pi;
%Creates end points of dislocations on [1 1 1]
x1=dis_length/2*(-cos(t)/sqrt(2)-sin(t)/sqrt(6));
y1=dis_length/2*(cos(t)/sqrt(2)-sin(t)/sqrt(6));
z1=dis_length/2*(2*sin(t)/sqrt(6));
%Creates end points of dislocations on [1 1 -1]
x2=dis_length/2*(-cos(t)/sqrt(2)+sin(t)/sqrt(6));
y2=dis_length/2*(cos(t)/sqrt(2)+sin(t)/sqrt(6));
z2=dis_length/2*(2*sin(t)/sqrt(6));

Junction Length Calculation

Figure 3: Typical zipped binary junction

The length of the dislocation junction along the line of intersection is thought to
be related to the strength of the Lomer-Cottrell lock. Therefore, we desire to
calculate the junction length after the dislocations have been allowed to relax. We
also want this calculation to be automated, because we hope to simulate 441
configurations. It is obvious from Figure 3 that the nodes at the ends of this
"zipped" binary junction are unique in that they possess more than two arms.
Calculating the length of the junction is as simple as finding the magnitude of the
vector between the two free nodes possessing three or more arms from the first
column of the linksconnect matrix. However, if the dislocations are oriented in
such a way that junction formation is unfavorable, they could remain crossed at a
single point, as illustrated in Figure 4. If there are less than two free nodes with
three or more arms, then the junction length is zero. A simple algorithm in MATLAB
is:
inter_nodes = find(-10^-2<rn(:,3) & rn(:,3)<10^-2 );
inter_nodes_seg = connectivity(inter_nodes,1);
junc_end_nodes= find(inter_nodes_seg>2);
if length(junc_end_nodes)<2
junc_length(dis1_no,dis2_no)=0;
else
junc_length(dis1_no,dis2_no)=norm(inter_coord(junc_end_nodes(1),1:3)...
-inter_coord(junc_end_nodes(2),1:3));
end

Figure 4: Unfavorable junction orientation

This algorithm works when the dislocations have only partially zipped a junction,
which encompasses the majority of situations. However, if one of the dislocations
is in the
or
direction, it is coincident with the intersection line
and, thus, the dislocations have the possibility to fully zip the junction, as shown
in Figure 5. In these situations a fixed end node is also the end node of the zipped
junction. The end nodes of the junction now only possess two arms, because the
fixed nodes were placed at the end of the dislocation. Therefore, the
aforementioned algorithm would find the junction length to be zero, rather than
the entire length of the original dislocation (4000). Since the fixed nodes possess a
flag in rn, a simple improvement to the previous algorithm is:
fixed_nodes = find(rn(:,4)~=0);
inter_fixed_nodes = find(inter_coord(:,4)~=0);
...
inter_nodes_seg(inter_fixed_nodes)=inter_nodes_seg(inter_fixed_nodes)
+1;

Figure 5: Fully zipped junction

This would add one to the connectivity of the end nodes such that the value of a
zipped end node is three, which is detected as a junction end node.

Figure 6: Repulsive coincident dislocations

Unfortunately, this algorithm still fails to address when both dislocations overlap,
i.e. they both lie along the intersection line. We know physically that the
dislocations would combine to form a single dislocation if they are attractive. This
is exactly what will happen in the simulations, because DDLab and ParaDiS both
handle topological changes. Therefore, this junction length should be considered
the entire length of the dislocation, but the previous algorithm would find it to be
zero, because there are no nodes with more arms than they originally possessed.
However, when the dislocations are repulsive the results in DDLab do not
correspond to reality. The dislocations are not able to simply move apart, because
their end nodes are pinned. The DDLab results for this situation can be seen in
Figure 6. In this case, only the end nodes will merge. An entirely different
approach must be used to calculate junction length for these special cases.
Initially, rn had four fixed nodes with a flag of 7; after the topological changes,
there are only two fixed nodes, which is unique to these scenarios. In the case of
repulsion the fixed nodes have two arms, while they only have one in attraction.
These additions can be seen in LC_junction.m

Batch Execution
Executing a batch of 441 jobs in DDLab is very simple due to the versatility of
MATLAB. A separate script named LC_junction.m was created to call the DDLab
input deck LC_relax.m for all iterations of the dislocations. Although it is possible
to place the loops inside a DDLab input deck, it is best to keep the DDLab portions
of the code separate from the other MATLAB scripts and retain the standard DDLab
file layout.* LC_junction.m then stores the resulting nodal positions, rn; segment
data, links; and junction length, junc_length, for all iterations in the MATLAB
binary file LC_juction_data.mat for later analysis, and also as a restart point for
the dissociation simulations explained in the next section.

*It should also be noted that the input decks for ParaDiS cannot be modified from
their usual format.

Resolved Shear Stress Calculation

Figure 7: Coordinate system for resolved shear stress

After the dislocations are allowed to relax and form junctions. A shear stress is
applied on the
plane in order to determine the resolved shear stress
necessary to break the junction. Thus far, we have been using the Cartesian
coordinate system with basis
,
, and
A stress transformation is needed in order to apply the shear directly on the

plane. The new orthogonal coordinate system was chosen as

,

, and

, as illustrated in Figure 7.

The stress transformation can be performed by the relation:

where

and

are the stress tensor defined in the original and new coordinate

systems, respectively.
is the transformation matrix, where its columns are the
components of the basis written in terms of the basis . If and are chosen
as orthonormal bases, then

is an orthogonal matrix consisting of directional

cosines, and

equals

The resolved shear stress,

in basis (Figure 7) is applied to the formed binary
junction. This stress is transformed using the following algorithm in LC_stress.m:
%stress in coordinate system 2
sigma = [ 0
0 tau_13
0
0 0
tau_13 0 0 ] * 1e6;

%in Pa

%coordinate system 1 (cubic coordinate system)

e1p = [1 0 0]; e2p = [0 1 0]; e3p = [0 0 1];
e1p=e1p/norm(e1p); e2p=e2p/norm(e2p); e3p=e3p/norm(e3p);
%coordinate system 2
e1 = [-1 -1 2]; e2 = [-1 1 0]; e3 = [-2 -2 -2];
e1=e1/norm(e1); e2=e2/norm(e2); e3=e3/norm(e3);
%rotation matrix
Q = [ dot(e1,e1p) dot(e2,e1p) dot(e3,e1p)
dot(e1,e2p) dot(e2,e2p) dot(e3,e2p)
dot(e1,e3p) dot(e2,e3p) dot(e3,e3p) ];
%Transform sigma into coordinate system 1
appliedstress = Q*sigma*Q';

Figure 8: Dissociation of a Lomer-Cottrell lock

This stress causes the junction to unzip by the bowing out of the dislocation arms,
until a force equilibrium is reached. If the stress happens to be greater than the
critical stress, then the binary junction will dissociate by fully unzipping the
dislocation junction, as illustrated in Figure 8. To avoid developing a complex
iterative scheme to find the critical stress, this tutorial asks the user to input the
magnitude of the applied stress until a satisfactory stress near the critical resolved
shear stress is found.
The strength of the junction has been shown to depend on the length of the
dislocation arm, , because the arms bow out like a Frank-Read Source, as shown
in the sequence of Figure 8. The critical resolved shear stress to break the junction
can be roughly obtained by the relationship:

where is the shear modulus and is the magnitude of the Burgers vector. From
the relation, it is obvious that dissociation can occur more easily with longer
dislocation arms. If the initial length of the two dislocations increases, the
resulting junction length increases. Also, the length of dislocation arms become
longer (not shown), resulting in lower critical resolved shear stresses.

Description of Files
This section provides a brief description of all the files used in this tutorial and
their hierarchy. The main bullets are the primary scripts, and the sub-bullets are
the functions or scripts they call.

LC_Junction.m- This script calculates the length of the junction formed

between the two dislocations in the absence of stress, as described in the

Junction Length Calculation section. It generates the surface plot of
normalized junction length as a function of orientation shown in Figure 9. It
stores the final dislocation data and junction length ( rn, links, and
junc_length) in LC_junction_data.mat.
o LC_make_dis.m- The script that generates the coordinates of one
fixed node for all 21 initial orientations of both dislocations from
to . The midpoint of all the dislocations is the origin
o LC_relax.m- Script with the input deck for DDLab with zero applied
stress. It requires D1 and D2, one fixed node position for dislocations
1 and 2 respectively; and dis_length, the initial length of the
dislocations, as inputs.
o LC_junc_length.m- Function that calculates the length of the
resulting binary junction. It takes the resulting rn and connectivity
matrices from LC_relax.m as inputs and returns only junction length.

LC_dissociate_orientation.m- This script determines the critical stress to

dissociate a binary junction based on the orientation of the dislocations.

The orientation angle of dislocation 1,

dislocation 2,

, equals

and those of

, are

,
, and
. The junction data is loaded from
LC_junction_data.mat. The user inputs a guess at the critical shear stress
in coordinate system 2 (see the Resolved Shear Stress Calculation section).
The resulting critical shear stresses are plotted with
and junction length.
It is best to start with a conservative guess to avoid having the dislocations
dissociate and begin acting as Frank-Read sources, which is costly
computationally.
o LC_stress.m- Script with the DDLab input deck for the user specified
in coordinate system 2. It performs the stress transformation
described in the Resolved Shear Stress Calculation section. It also
requires that dis_length, the initial length of the dislocations prior
to junction formation, is specified.

LC_dissociate_length.m- This script determines the critical stress to

dissociate a binary junction based on the initial length of the dislocations.

and
are both chosen as
. The user inputs a guess at the critical
shear stress in coordinate system 2 (see the Resolved Shear Stress
Calculation section). The resulting critical shear stresses are plotted with
junction length.
o LC_relax.m- See above.
o LC_junc_length.m- See above.
o LC_stress.m-See above.

Results and Conclusions

Junction Formation

Figure 9: Normalized junction length

Figure 9 shows the normalized resulting junction length between the dislocations
under no external stresses. The junction length is highly dependent on both

and

. These results compare well with those published by Madec et. al.<ref
name="Madec"/> Closer examination of the results reveals that reversing the line
direction, , of both dislocations creates junctions with the same length.
Although, the situations are not truly symmetric because the Burgers vectors do
not flip, the lack of external stress means that the only forces are due to the
interaction of the dislocations. Therefore in terms of junction formation, there is
no difference between the junction of two dislocations and that of their opposites.
For example, when

and

0.5862, which is the same as when

their normalized junction length is

and

. This means the

entire plot could have been generated by only varying

and

over the range

, thus cutting the number of trials in half.

Junction Dissociation
Note: These are only preliminary results. Feel free to improve upon them. For the
sake of time only a few junctions were manually dissociated.
The first comparisons in LC_dissociate_orientation.m are to see the effect of
orientation on the junction strength.

is fixed at

and

is varied from

to

. Figure 10 shows the critical resolved shear stress compared to

. The stress
appears to decrease linearly with increasing angle, but the data are insufficient.
The decrease in stress is expected, because the junction length decreases with
increasing angle which means the that dislocation arms, , are longer. Figure 11
compares the critical stress with the junction length. Again the relationship
appears linear, but no conclusion can be made.

Figure 10: Critical stress compared

with orientation angle

Figure 11: Critical stress compared with

junction length at different orientations

The effect of the initial dislocation length on the junction strength is studied in
LC_dissociate_length.m The dislocation arm length, should scale with junction
length; thus, the longer dislocations should have lower critical resolved shear
stresses. The junction length scales linearly with dislocation length, as illustrated
in Figure 12. As expected, the junction strength decreases with increased
increased junction length. Figure 13 shows the critical resolved shear stress
appears to scale with , which is predicted by the relationship in the Resolved
Shear Stress Calculation section.

Figure 12: Comparison of junction length

with initial dislocation length

Figure 13: Critical stress at various

dislocation lengths

We can discuss in more detail which dislocation length is the controlling parameter
for critical stress. Right now, the relaxed dislocation structure has five segments

and their lengths are all related in this particular series of simulations. More
simulations with different geometries are needed to determine which length is the
controlling parameter.

Junction Dislocations in DDLab

Figure 14: Bowing of the binary junction should be avoided

Physically, the dislocation formed by the junction between two intersecting
dislocations on differing slip planes is sessile, because it is composed of a partial
dislocation in both slip planes. If one of the partials were to glide in its respective
slip plane, it would force the other to climb. However in DDLab, the junction is
formed using its normal split and merge topological changes, which means the
nodes of the junction are not sessile. These operations enforce Burgers vector
conservation, but calculate the slip plane by taking the cross product of the
resulting Burgers vector and line direction. In the case of this tutorial,
and
for the junction dislocation. The nodes at the
end of the junction are constrained to only move along the line of intersection,
because they are connected to segments in both slip planes. On the other hand,
the nodes on the interior of the junction do not have such a restraint and are free
to move in the
plane. If an external shear stress is applied to the junction,
a situation such as that in Figure 14 can result. This bowing of the junction is
clearly not physical and should be avoided.

References
<references/> [1] From Dislocation Junctions to Forest Hardening, Phys. Rev. Lett.
89 255508 (2002), R. Madec, B. Devincre, and L.P. Kubin
Currently, the references only show up automatically to logged in users. We are
looking for a solution to this problem. For now we are just writing the references
by hand.

Downloads
LC_Junction_DDLab.tar.gz - DDLab tutorial files. (Note: best if extracted as a
subdirectory of the DDLab Inputs folder)
Note: there are currently some issues with the fcc mobility law in the latest
release of DDLab (12/18/2007) that we hope to work out in the future. There is
also a typo in mobfcc0.m involving BScrew, Bedge, etc. If you have any trouble
resolving this issue please contact the author for a corrected version of the code.

LC_lock.pdf. - Past version of this tutorial. (Note: There are some errors associated
with the junction length calculation and the codes have been significantly
modified)

* Please contact William Cash with any questions or concerns

DDLab has a function consistencycheck.m that checks the self-consistency of the
input dislocation structure before starting the DD simulation. The function
consistencycheck.m is called near the very beginning of dd3d.m
This function checks not only the input structures rn and links but also the
auxiliary data structures connectivity and linksinconnect, which are created by
function genconnectivity.m.
In the version ddlab-2007-12-18.tar.gz or newer, the consistencycheck.m does
the following checks. Items 3, 4, 5 are general checks that are in principle
necessary for any DD simulation program. Other items are specific to the DDLab
data structure and algorithm implementation.
1. the length of the connectivity matrix should be the same as that of
links
2. connectivity matrix: each row correspond to a node,
only the first 2*numNbrs+1 entries in each row can be non-zero
where numNbrs is the number of neighbors for this node

3. Burgers vector conservation at every node (if flag != 7)

4. a node cannot be connected to another node twice
5. a link cannot have zero Burgers vector
6. a node cannot be connected to the same link twice in the connectivity
matrix
7. consistency between "connectivity" matrix and "links" matrix
if node i has a link j, then link j must have node i as one of its
end nodes
8. check the consistency between the "connectivity" matrix and the
"linkinconnect" matrix