Lagrangian Formulation of Nodal Dynamics

– Abacus Model^{∗}

Ju Li

April 15, 1999

In my problem – dislocation motion impeded by a random distribution of point-like obstacles, it may be still better to discretize the x-direction by regular, equally h-spaced mesh. However I am no longer assuming ∂y/∂x to be small, thus the canonical equation of motion (assuming no lattice barrier) for “nearly straight” dislocation^{1},

m\frac{∂^{2}y}{∂t^{2}} = ϵ\frac{∂^{2}y}{∂x^{2}} − B \frac{∂y}{∂t} + σb − F (x, y|c, a, K)  (1)

is no longer good. I shall briefly explain the physical parameters appearing in (1), whose meanings hold true in general even when (1) does not.

• m is the dislocation “mass” density ([kg/m]), which is related to the kinetic energy by dK = (ds)mv_{n}^{2}/2, where ds is the differential line length ([m]) and v_{n} is its “apparent”, normal velocity ([m/s]). Of course, dislocation is not a physical object and does not have mass in the literal sense; it does, however, describes a peculiar arrangement of atoms near it which carries quanta of materials deformation and has far-field elastic stress/strain. Thus, if the dislocation moves, the stress/strain field is going to be time-dependent, thus the atoms are going to throb up and down just like a boat in wavy water. So dK is actually the kinetic energy of atoms sorrounding the moving dislocation.

∗This is a term project for MIT 18.336, Numerical Methods of Applied Mathematics II, taught by Prof.

M. Brenner in fall, 1998. Links to complete source codes, results, and postscript file of this writeup can be found at http://www.mit.edu/∼liju99/paper archive.html.

1J. P. Hirth and J. Lothe, Theory of Dislocations, Wiley/Interscience, 1982; U. F. Kocks, A. S. Argon, M. F. Ashby, Thermodynamics and Kinetics of Slip, Pergamon Press, 1975.

• ϵ is the dislocation line tension or energy ([J/m]), which is the sum of core energy and elastic energy outside the core. In a not-too-bad idealization called constant line energy model, ϵ = μb^{2}/2 is a constant (does not depend on orientation), where μ is the shear modulus ([J/m^{3}]) of the isotropic continuum and b is the Burger’s vector ([m]) of the dislocation. Interestingly, one can derive for screw dislocation: ϵ = mv_{s}^{2}, in close analogy with the special relativity formula E = mc^{2}, where c is replaced by the shear wave speed v_{s} = \sqrt{μ/ρ} ([m/s]), and ρ is mass density of the material ([kg/m^{3}]). It is thus very convenient to use a reduced unit system in which the length unit is b, the velocity unit is v_{s}, and the energy unit is ϵb; from which all other units (time [b/v_{s}], stress [ϵ/b^{2}], ..) can be defined. And this unit system m = 1, b = 1, ϵ = 1 in (1).

• B is the phonon drag coefficient ([J·s/m^{3}] or [ϵ/v_{s}b]), which gives dissipation to the dislocation motion. It is similar in mechanism to the resistance a man feels when walking in the rain, who receives more raindrops at front than back. B can be shown to be on the order of k_{B}T/v_{s}b^{2}, where T is the material temperature.

• σ is the external stress ([J/m^{3}] or [ϵ/b^{2}]). It makes the dislocation move normal to itself to increase the area swept by since the material can be deformed that way, and σ likes the material to be deformed because the total free-energy is reduced.

• F (x, y|c, a, K) is the resistance force field due to impurities in the material. Here c is the volume ratio of the impurity, a is its average size (or force range) ([m] or [b]), and K is the maximum resistance the impurity offers ([J/m] or [ϵ]). There exist two opposite limits to F: one is c ≫ K/2ϵ, called “Mott statistics”, corresponding to weak but populous obstacles (solute atom strengthening). The theoretically predicted critical stress in this limit^{2} is

σ_{c} = \frac{2}{3} (\frac{K}{2ϵ} )^{4/3} c^{2/3}μ,  (2)

above which the dislocation can move steadily in spite of the obstacles and the material can deform plastically. The other limit, c ≪ K/2ϵ < 1, is called “Friedel statistics”^{3} (precipitate-particle strengthening), corresponding to fewer but stronger obstacles, which has a predicted flow stress of

2The theory was first proposed by Mott and Nabarro (1948), and later was analyzed by Schwarz and Labush, J. Appl. Phys. 49, 5174 (1978).

σ_{c} = \frac{2ϵ}{bL} (\frac{K}{2ϵ} )^{3/2} = \frac{2ϵ}{ba} (\frac{3c}{2π} )^{1/2} (\frac{K}{2ϵ} )^{3/2}  (3)

where L^{2} is the average area an obstacle “owns” on the slip plane ([m^{2}] or [b^{2}]), related to c and a. (3) is derived based on a geometrical argument by Friedel, who also predicts that at the critical point,

λ_{c} = L (\frac{K}{2ϵ} )^{−1/2} ,  (4)

where λ is the average separation between two nearest obstacles impeding the same dislocation.

I will focus on c ≪ K/2ϵ limit first, in testing my own numerical scheme, since (3) was well-confirmed by experiments and simulations and is also conceptually convincing; while the theoretical derivation of (2), to me, looks suspicious. The main significance of (3) is the scaling law

σ_{c} ∝ c^{1/2},  (5)

which says that the plastic strength of an alloy does not scale linearly with the concentration of strengthening agents, thus by adding a small amount of agents one can greatly enhance its strength (making dislocation harder to move and so the flow stress is increased); also, from (3) one can show that peak strength can be obtained by controlling precipitate-particle size a (as K also depend on a), through a thermal process called “aging”.

Now returning to (1), which in the reduced unit system is

\frac{∂^{2}y}{∂t^{2}} = \frac{∂^{2}y}{∂x^{2}} − B \frac{∂y}{∂t} + σ − F (x, y|c, a, K)  (6)

where actually everything is dimensionless in (6) but I am too lazy to add the˜ . Since the obstacles are now point-like, it makes sense to use its size a as the x-mesh size h. That is, I am going to neglect all details happening below a, effectively treating obstacle-dislocation interaction as a simple “arrest/free” rule, where a single mesh point (ih, y_{i}) is going to be

3A theory neglecting mass and dissipation was first proposed by Friedel (1956). A review was given by Ardell, Metal. Trans. A16, 2131 (1985).

Figure 1: Dislocation impeded by the obstacle, forming a cusp.

detained at the obstacle position (x_{k}^{o}, y_{k}^{o}), if it is found trying to cross y_{k}^{o} but the flanking angles on two sides are not enough to break through the obstacle, i.e.,

θ_{1} + θ_{2} < \frac{K}{ϵ}  (7)

in Fig. 1.

So, a cusp is going to form at the obstacle position (x_{k}^{o}, y_{k}^{o}), and may be released later if σ is large enough to make the flanks bow out beyond K/2ϵ, or if the inertia is large enough to induce that. The accuracy and stability of the numerical scheme is of foremost concern here. Although I am going to adopt periodic boundary condition, using the spectral method is out of the question because it loses all its magic when cusps exist in the solution. How about finite-difference schemes on a regular mesh? Surprisingly, the first thing I tried, an explicit, central-difference scheme, works well, for reasons I will explain later. The scheme is just

\frac{y_{i}^{n+1} − 2y_{i}^{n} + y_{i}^{n−1}}{Δt^{2}} = \frac{y_{i+1}^{n} − 2y_{i}^{n} + y_{i−1}^{n}}{Δx^{2}} − B \frac{y_{i}^{n+1} − y_{i}^{n−1}}{2Δt} + σ + O(Δx^{2}) + O(Δt^{2}), (8)

Figure 2: Numerical solution of cusp evolution using finite-difference scheme (9) with r = 0.9.

which leads to

y_{i}^{n+1} ≈ \frac{(y_{i+1}^{n} − 2y_{i}^{n} + y_{i−1}^{n})r^{2} + σΔt^{2} − (1 − BΔt/2)y_{i}^{n−1} + 2y_{i}^{n}}{1 + BΔt/2}  (9)

if node (ih, y_{i}^{n}) is free, and

y_{i}^{n+1} ≡ y_{k}^{o}  (10)

if otherwise. That is, the arrest and release of the dislocation node (segment) by the obstacle is assumed to take on the smallest length/time scales, whose details are of no relevance to the overall behavior. Here

r = \frac{v_{s}Δt}{Δx} (v_{s} = 1)  (11)

should be less than 1, by the Courant-Friedrichs-Levy (CFL) convergence condition^{4}. Fig. 2 is a simulation of the evolution of an initially straight triangular cusp between x = [0, 2π] (under PBC) from t = 0 to t = 1, with B = σ = 0, in comparison with the exact solution. One can see that (9) is stable and works quite well.

4A. Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996.

Why do I insist on using a regular x-mesh? Because it simplifies the computation of obstacledislocation interactions. Say there are M mesh-points and W obstacles, how can I determine which nodes are crossing which obstacle at a given time? If my nodes are specified by (x_{i}, y_{i}), i.e., they have two degrees of freedom, the calculation may take O(MW) efforts^{5} per step. But, if I use a regular x-mesh, x_{i} = ih, and make a not-too-restricting assumption that all obstacles also lay on integer multiples of h, x_{k}^{o} = n_{k}h, then instead of letting the dislocation seek out obstacles, I can let the obstacles seek out the dislocation, i.e., check if the kth obstacle is crossed by node i = x_{k}^{o}/h or not, which takes only O(W) effort to do per step.

Before moving on, let me add that there are two other ways a dislocation may overcome an obstacle in addition to “simple breaking through” indicated by (7). The first one is “Orowan loop formation”, which happens when the obstacles is impenetrable (dispersalparticle strengthening) by the dislocation, K/2ϵ > 1, but the external stress σ is large enough to make θ_{1} and θ_{2} approach π/2; then the two flanks merge ahead of the obstacle – a new dislocation is formed and marches on, leaving a little dislocation loop behind surrounding the obstacle. The actual process is quite complicated because of dislocation self-interaction, which is ignored in our treatment; so I won’t lose much by simply stipulating that whenever θ_{1} + θ_{2} > π/3 during the simulation, the dislocation is freed by Orowan loop shedding. The other mechanism is thermal activation, which means that even if K − F > 0, there is still a finite probability that an arrested segment can spontaneously overcome the energy barrier

E ≈ a · (K − F)  (12)

of the obstacle with the help of thermal fluctuation, much like tunneling in quantum mechanics. The rate of transition ν is on the order of

ν ≈ \frac{ν_{D}}{n} exp (−\frac{E}{k_{B}T} ) ,  (13)

where ν_{D} ≈ v_{s}/b is the Debye (characteristic atomic vibration) frequency of the solid and n is the number of atoms involved in this collective transition, which is ∼ a/b in our case. Thus, even if E > 0 and the segment can not be freed by force, there is still probability 1 − e^{−νt} that the trapped segment could be released by thermal activation. These two features can be easily incorporated into the simulation.

5Unless something clever, like dividing the cell into small bins and maintaining a link-list, is done. But I cannot afford that complexity in my pilot project.

The problem now is to generalize (1), which is only valid for “nearly straight” dislocations, to situations where |∂y/∂x| is O(1). But I am still going to adopt a regular x-mesh with spacing h due to above-mentioned efficiency considerations, thus comes by the name “abacus model”. So long as |∂y/∂x| does not approach ∞, the abacus discretization should work fine as (9). But as I find out, it is not as easy to come up with a stable scheme. Some seemingly plausible schemes works fine for smooth y(x, t), but is unable to cope with the cusp situation and the quickly becomes unstable. Some schemes may work fine but are implicit, requiring at least a tridiagonal linear system solver^{6}. I am looking for a robust scheme that is explicit, has good physics in its form and does not have to provide too high accuracy because neglecting dislocation self-interaction is already a sin.

To incorporate the correct dynamics, we must know the Lagrangian

L = K − V,  (14)

and minimize the action

A = ∫ L(t)dt.  (15)

L is itself a line integral over ds, and it is more convenient to discretize it that way first, as a function of {y_{i}(t), i = 1..M} where I stipulate that (ih, y_{i}(t)) is right on the dislocation at time t. Remember that (ih, y_{i}(t)) is not a real object but collectively they represent a curve, so it makes no sense to identify (ih, y_{i}(t)) as a physical object. Also, it makes no sense if the mesh point moves along the dislocation line because the curve they represent does not change.

In the constant line energy model, the elastic energy is simply

V_{e} = ∫ ϵds = ∫ ϵ\sqrt{dx^{2} + dy^{2}},  (16)

6It is worse than that: because the meshes wrap around (PBC), the matrix is actually tridiagonal plus two elements on the corners, which could be solved “quickly” by moving the two elements to RHS and do a self-consistent iteration. But I am too lazy to code it up.

which I discretize as

V_{e} ≈ ∑_{i} ϵ\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}} + O(h^{2}).  (17)

I tried to be clever once and discretized it as

V_{e} ≈ ∑_{i} ϵ\sqrt{h^{2} + (y_{i+1} − y_{i−1})^{2}/4},  (18)

thinking it may provide better accuracy. But the resultant numerical scheme was found to have a “crumpling” instability for short-wavelength fluctuations, because (18) does not provide enough tension to straighten up the curve. The reason is because (17) is a physicallybased discretization where we “sew” straight elastic strings between nodes, which though crude, provides real cohesion. But what physical model (18) represent, with its central difference look? Strings with detached ends? Clearly there exists modes which deform the nodes, such as y_{i} = (−1)^{i}, but which does not receives energy penalty from (18) at all, that is the cause of the “crumpling” instability.

The kinetic energy K is more complicated than V_{e}, which formally is

K = ∫ ds · mv_{n}^{2}/2,  (19)

where v_{n} is the normal, or “apparent” velocity of the dislocation. At the ith node, I approximate the normal vector n as

n_{i} ≈ \frac{[(y_{i−1} − y_{i+1})/2, h]}{\sqrt{h^{2} + (y_{i+1} − y_{i−1})^{2}/4}} = [−\frac{d_{i}}{s_{i}} , \frac{h}{s_{i}}],  (20)

where I define

d_{i} ≡ \frac{y_{i+1} − y_{i−1}}{2} , s_{i} ≡ \sqrt{h^{2} + d_{i}^{2}},  (21)

thus

v_{n}(ih) = v_{i} · n_{i} = (ẏ_{i})(n_{i}^{y}) ≈ (ẏ_{i})(\frac{h}{s_{i}}).  (22)

Also, let me approximate ds centered around node i to be

ds ≈ s_{i},  (23)

so,

K ≈ ∑_{i} ms_{i}(\frac{h}{s_{i}})^{2}(ẏ_{i})^{2}/2.  (24)

How dare I use central difference here? I can fend off by mumbling: “energetics is zerothorder, dynamics is first-order”, but it does not convince people. The actual reason is because it is quite difficult to make the effective mass go to zero in (24), but we can invent a mode in (18) to make things immediately go wrong (dV_{e} ≡ 0) without any initial resistance. With the correct (17) straightening out the curve, (24) is in fact a very good approximation. However,

(24) is not consistent with the physical picture depicted by (17) of nodes linked by straight segments, whose kinetic energy should be

K = ∑_{i} \frac{m}{2} \frac{h^{2}(ẏ_{i}/2 + ẏ_{i−1}/2)^{2}}{\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}}} .  (25)

I am convinced that (25) would lead to a stable, and better scheme; but computationally it is a nightmare, which as I mentioned, is implicit and involves tridiagonal linear system and iterative relaxations. So I have to give it up and use (24) instead – and let me note here that

\frac{∂s_{i}}{∂d_{i}} = \frac{d_{i}}{s_{i}} , \frac{∂s_{i−1}}{∂y_{i}} = \frac{d_{i−1}}{2s_{i−1}} , \frac{∂s_{i+1}}{∂y_{i}} = −\frac{d_{i+1}}{2s_{i+1}} ,  (26)

from (21).

So for an elastic curve with inertia, without dissipation, the equations of motion of “abacus” mesh points are

\frac{∂}{∂t} (\frac{∂L}{∂ẏ_{i}} ) = \frac{∂L}{∂y_{i}} .  (27)

Here

\frac{∂L}{∂ẏ_{i}} = \frac{∂K}{∂ẏ_{i}} = \frac{mh^{2}}{s_{i}} ẏ_{i},  (28)

and so

\frac{∂}{∂t} (\frac{∂L}{∂ẏ_{i}} ) = \frac{mh^{2}}{s_{i}} ÿ_{i} − \frac{mh^{2}}{s_{i}^{3}} d_{i}ḋ_{i}ẏ_{i}.  (29)

On the other hand,

\frac{∂L}{∂y_{i}} = \frac{∂K}{∂y_{i}} − \frac{∂V_{e}}{∂y_{i}} ,  (30)

with

\frac{∂K}{∂y_{i}} = \frac{mh^{2}}{2} (\frac{∂s_{i−1}^{−1}}{∂y_{i}} (ẏ_{i−1})^{2} + \frac{∂s_{i+1}^{−1}}{∂y_{i}} (ẏ_{i+1})^{2}) = \frac{mh^{2}}{2} (\frac{d_{i+1}(ẏ_{i+1})^{2}}{2s_{i+1}^{3}} − \frac{d_{i−1}(ẏ_{i−1})^{2}}{2s_{i−1}^{3}} ) , (31)

and

\frac{∂V_{e}}{∂y_{i}} = ϵ (\frac{y_{i} − y_{i−1}}{\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}}} + \frac{y_{i} − y_{i+1}}{\sqrt{h^{2} + (y_{i} − y_{i+1})^{2}}} ) .  (32)

What is the effect of external stress? Say y_{i} is increased by amount δy_{i}, then the area swept by the dislocation increases by hδy_{i}, so the external stress would have done work σb · hδy_{i} – thus the RHS of (27) should include an external stress term σbh.

What is the effect of phonon drag? Say that the line segment ds feels dragging force

dF = −ds · Bv_{n}n,  (33)

decomposed into the y-axis, that is

dF_{y} = −ds · Bv_{n}n_{y} = −Bv_{n}h,  (34)

thus, following (22), we should add a force term −B(ẏ_{i})(h/s_{i})h to RHS of (27), which is very reasonable since, if the dislocation slope is large, ∂y/∂x ≫ 1, even if ẏ_{i} is quite big, there should not be a significant phonon drag.

So let me put all the pieces together and the equation now reads

\frac{mh^{2}}{s_{i}} ÿ_{i} − \frac{mh^{2}}{s_{i}^{3}} d_{i}ḋ_{i}ẏ_{i} = \frac{mh^{2}}{2} (\frac{d_{i+1}(ẏ_{i+1})^{2}}{2s_{i+1}^{3}} − \frac{d_{i−1}(ẏ_{i−1})^{2}}{2s_{i−1}^{3}} ) +  (35)

ϵ (\frac{y_{i−1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}}} + \frac{y_{i+1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i+1})^{2}}} ) + σbh − B(ẏ_{i})(\frac{h}{s_{i}})h.

Dividing both sides by h^{2}/s_{i}, we have

mÿ_{i} − \frac{m}{s_{i}^{2}}d_{i}ḋ_{i}ẏ_{i} = \frac{ms_{i}}{2} (\frac{d_{i+1}(ẏ_{i+1})^{2}}{2s_{i+1}^{3}} − \frac{d_{i−1}(ẏ_{i−1})^{2}}{2s_{i−1}^{3}} ) +  (36)

ϵ\frac{s_{i}}{h^{2}} (\frac{y_{i−1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}}} + \frac{y_{i+1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i+1})^{2}}} ) + σb · \frac{s_{i}}{h} − Bẏ_{i}.

Let me define

σ′b ≡ σb · \frac{s_{i}}{h} + \frac{ms_{i}}{2} (\frac{d_{i+1}(ẏ_{i+1})^{2}}{2s_{i+1}^{3}} − \frac{d_{i−1}(ẏ_{i−1})^{2}}{2s_{i−1}^{3}} ) ,  (37)

and

B′ ≡ B − \frac{m}{s_{i}^{2}}d_{i}ḋ_{i} = B − \frac{md_{i}(ẏ_{i+1} − ẏ_{i−1})}{2s_{i}^{2}},  (38)

so now we have

mÿ_{i} = ϵ\frac{s_{i}}{h^{2}} (\frac{y_{i−1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i−1})^{2}}} + \frac{y_{i+1} − y_{i}}{\sqrt{h^{2} + (y_{i} − y_{i+1})^{2}}} ) + σ′b − B′ẏ_{i},  (39)

which, compared with (1), suggests a finite-difference scheme by just taking σ′ as σ and B′ as B in (9). Of course, σ′ and B′ are not constants here but depend on i and t_{n}, and in evaluating (37), (38) we need to use ẏ_{i±1}(t_{n}) which, unless one enjoys doing full-blown Newton’s iteration with lots of indices, can only be explicitly computed by forward difference approximation,

ẏ_{i±1}(t_{n}) ≈ \frac{y_{i±1}^{n} − y_{i±1}^{n−1}}{Δt} + O(Δt),  (40)

which ruin the second-order accuracy of (9) but match the first-order accuracy^{7} of (17). With cusps everywhere, self-interaction ignored and constant line energy model, this is probably as good as we deserve. (37)-(40), in combination with (9), and “arrest/free” rules (Orowan loop formation, simple breakaway, thermal activation), define the numerical scheme that I will use below to simulate precipitate-particle strengthening.

Before looking into the results, let me define a prototype system, which is Cu (A=63.546) strengthened by Co (A=58.933) or Fe (A=55.845) precipitate-particles of volume concentration 0.0001 − 0.05 and K/2ϵ = 0.05 − 1. The particles have size 5 − 50b and are coherent Figure 3: Overall view of a simulation: dash rectangle shows the PBC cell; small green dots are the obstacles, which turns red/big if in contact with dislocation.

7One may improve upon (17) though, by inferring curvature through interpolation, with special consideration to cusps.

with the matrix, with size misfit on the order of 1%. At T = 300K, Cu has shear modulus μ = 48 GPa and density ρ = 8920 kg/m^{3}. ρ is taken to be independent of T, from which we infer that the Buger’s verctor b = 2.56 A; but μ is in fact quite sensitive to the temperature. Using a rule of thumb that μ is reduced to half its T = 0K value at the melting point (for Cu, T_{m} = 1357.77K), we assume that

μ(T) = μ(0) (1 − \frac{T}{2T_{m}} ) ,  (41)

with μ(0) = 54 GPa. The point is, our fundamental units ϵ and v_{s} depend on μ, and so they will depend on T. For numbers, at T = 300K, ϵ ≈ 0.872 eV/A and v_{s} = 2187.8 m/s; and one reduced stress unit = μ(T)/2 ≈ 24 GPa, etc.

The code is written in C. In order to visualize the simulation in real time I also programmed a threaded X Window 2D animator. Fig. 3 shows an overall snapshot of the simulation with cell size 40960b × 10240b (10μm × 2.4μm). The precipitate-particle size a = 10b, which is also the x-mesh size Δx, c = 0.001 and K/2ϵ = 0.2. Fig. 4 blows up a tiny region, where one can see the individual dislocation nodes.

To determine the critical flow stress σ_{c}^{sim} in my simulation, I calculate the dislocation average velocity ⟨v⟩ under different loading, and plot its ratio to the “full speed” σb/B – its speed if there were no obstacles. If the dislocation is hopelessly stuck during the simulation^{8}, Figure 4: Blowup of a tiny region: the white points are dislocation mesh points linked by blue lines to guide the eye. The green and red dots are obstacles.

8As there is thermal activation, the dislocation will eventually be freed; but since E in (13) depends on

then ⟨v⟩ is taken to be 0. Fig. 5 shows ⟨v⟩/(σb/B) at T = 300K, a = 10b, c = 0.001, K/2ϵ = 0.1, versus external stress in [ϵ/b^{2}] = [μ/2]. We can clearly see a threshold behavior, with threshold stress σ_{c}^{sim} ≈ 0.8×10^{−4}μ, which compares well with the theoretical threshold of σ_{c}^{sim} = 1 × 10^{−4}μ of (3), predicted by Friedel. Fig. 6 shows the behavior of λ, the average separation between nearest obstacles impeding the same dislocation, in comparison with the theoretical result (4). We see that again (4) is a highly accurate prediction; when λ/λ_{c} drops below 1, the dislocation loses its mobility.

Fig. 7 verifies the σ_{c} ∝ (K/2ϵ)^{3/2} scaling law for Cu at T = 300K and c = 0.001. The agreement between theory and simulation is seen to be excellent. Fig. 8 verifies the σ_{c} ∝ c^{1/2} scaling law for Cu at T = 300K and K/2ϵ = 0.1. The agreement is quite good for c > 0.001, but there seems to be a systematic deviation when c is very small. I attribute it to the inertia effect, since if the average separation between the obstacles is large, the dislocation can build up some momentum during its free flight to make the next jump easier. The Friedel theory is a static one, thus that effect is neglected. This deviation should be more pronounced at low T since B ∝ T.

loading, it takes very long; and soon after it is freed it will be trapped again. So the “stuck” behavior is quite easy to identify.

Figure 5: Average dislocation velocity at T = 300K, a = 10b, c = 0.001, K/2ϵ = 0.1. Figure 6: Average separation of nearest obstacles impeding the dislocation λ, divided by theoretical threshold (4), λ_{c}.

Figure 7: Verification of the σ_{c} ∝ (K/2ϵ)^{3/2} scaling law. Figure 8: Verification of the σ_{c} ∝ c^{1/2} scaling law.