s e b a s t i á n e c h e v e r r i r e s t r e p o
A T O M I S T I C A N D A R T I F I C I A L I N T E L L I G E N C E S I M U L A T I O N S O F G R A I N B O U N D A R I E S A N D D I S L O C A T I O N S
A T O M I S T I C A N D A R T I F I C I A L I N T E L L I G E N C E S I M U L A T I O N S O F G R A I N B O U N D A R I E S A N D D I S L O C A T I O N S
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universtiteit Delft,
op gezag van de Rector Magnificus Prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,
in het openbaar te verdedigen op maandag 02 april 2015 om 12:30 uur
door
s e b a s t i á n e c h e v e r r i r e s t r e p o
Dit proefschrift is goedgekeurd door de promotor: Prof. dr. B.J. Thijsse Samenstelling promotiecommissie: Rector Magnificus Prof.dr. K. Albe Prof.dr. S. Cottenier Prof.dr.ir M.G.D. Geers Prof.dr.ir E. van der Giessen Dr.ir. M.H.F. Sluiter
Prof.dr.ir. T.J.H. Vlugt
voorzitter
Technische Universität Darmstadt Universiteit Gent
Technische Universiteit Eindhoven Rijksuniversiteit Groningen Technische Universiteit Delft Technische Universiteit Delft
ISBN 978-90-9028921-2
© 2015, Sebastián Echeverri Restrepo
All rights reserved. No part of this book may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, without prior permission from the copyright owner.
C O N T E N T S
i n t r o d u c t i o n 1 i m i n i m i z a t i o n m e t h o d s 5 1 a t o m i s t i c r e l a x a t i o n o f s y s t e m s c o n t a i n i n g p l a s t i c i t y e l e m e n t s 7 1.1 Computational approach . . . 8 1.2 Algorithm benchmarking . . . 111.3 Average and local energy convergence . . . 13
1.4 Local stress convergence . . . 16
1.5 Average stress convergence . . . 18
1.6 Conclusions . . . 19
ii dislocations 23 2 g e n e r a t i o n o f d i s l o c a t i o n s 25 2.1 Removal of planes . . . 26
2.2 Application of a displacement field . . . 26
2.3 Atomistic relaxation . . . 28
2.4 Conclusions . . . 28
3 p e i e r l s s t r e s s 29 3.1 Methods for the calculation of the Peierls stress . . . 30
3.2 Simulation set-up . . . 32
3.3 Measurement of the Peierls stress . . . 33
3.4 Analysis . . . 34
3.5 Conclusions . . . 42
iii grain boundaries 49 4 a t o m i s t i c s i m u l a t i o n s o f g r a i n b o u n d a r i e s 51 5 g r a i n b o u n d a r y g e n e r a t i o n : t h e “ b r u t e f o r c e ” m e t h o d 55 5.1 Crystal orientations . . . 55
5.2 Microtranslations . . . 56
5.3 Atom removal . . . 59
5.4 Grain boundary selection . . . 59
5.5 Application of the algorithm . . . 61
5.6 Conclusions . . . 64
6 g r a i n b o u n d a r y g e n e r a t i o n : g e n e t i c a l g o r i t h m s 67 6.1 Energy function . . . 67
6.2 General description of the algorithm . . . 69
6.4 Application of the algorithm to an unknown grain boundary . . . 75
6.5 Conclusions . . . 78
7 p r e d i c t i o n o f g r a i n b o u n d a r y e n e r g i e s 79 7.1 Artificial neural networks . . . 80
7.2 Input data . . . 82
7.3 Training . . . 84
7.4 Results . . . 85
7.5 Conclusions . . . 87
iv grain boundaries & dislocations 93 8 g r a i n b o u n d a r y - d i s l o c a t i o n i n t e r a c t i o n 95 8.1 Grain boundary selection . . . 96
8.2 Computational approach . . . 97
8.3 Grain boundary width . . . .100
8.4 Dislocation spread . . . .104
8.5 Conclusions . . . .108
v a p p e n d i x 115 a m i n i m i z a t i o n a l g o r i t h m s 117 a.1 Dynamic relaxation methods . . . .117
a.2 Static relaxation methods . . . .120 b c r y s t a l o r i e n t a t i o n s e l e c t i o n 121 s u m m a r y 127 s a m e n v a t t i n g 131 a c k n o w l e d g e m e n t s 135 p u b l i c a t i o n s 137 c u r r i c u l u m v i t a e 139
I N T R O D U C T I O N
1Metallic interfaces are a key ingredient in controlling the strength, ductility,
reliability and lifetime properties of metal-based structural and functional ma-terials and devices. This holds for bulk mama-terials, where many properties are controlled by the behaviour at and through physical boundaries (precipitate hardening, Hall-Petch effect, etc.), and is even more so for small metallic com-ponents, like those used in Microelectromechanical systems (MEMS) where the mechanical properties are largely dominated by the metallic interfaces present (e.g. Grain Boundaries (GBs) and oxide layers). Likewise, solder interconnects are prone to failure along interfaces.
Studies of GBs, dislocations and their interaction, inevitably, have to be done at the atomic level, where the dynamics of individual atoms are expected to provide the answers that will help resolve the unknowns at larger scales.
Within the framework of the MuMIM (Multiscale Metallic Interface Mod-elling) project, which aims at constructing a generic multiscale technique for the quantitative description of metallic interfaces, this thesis aims at generating a framework that allows the simulation and understanding of GBs, dislocations and their interaction at the atomic scale.
A method that allows to deal with atomistic systems that are big enough to capture the behaviour of GBs and dislocations is Molecular Dynamics/Statics (MD/MS). Here, in order to maintain a reasonably good representation of the 1 Some extracts taken from the original description of the MuMiM (Multiscale Metallic Interface
properties of the material, the interaction of the atoms is defined via interatomic potentials based on the Embedded Atom Method (EAM).
MS allows to model the behaviour of atomistic systems in quasi-static condi-tions. but it has the mayor drawback that, since a full relaxation of the system has to be done after each deformation step, the calculations can become com-putationally expensive. Currently, there exist a variety of methods to perform this relaxation, but not much attention has been paid in the literature to the se-lection of the optimal choice for treating atomistic systems containing plasticity elements.
To tackle this issue, in Partian analysis and comparison of different minimiza-tion algorithms that can be used for molecular static simulaminimiza-tions is presented. The objective is to find an appropriate and efficient method for modelling GBs and dislocations at the atomic scale.
The study of single dislocations via MD/MS starts by the generation of an atomistic system containing such a defect. Once a dislocation is inserted, MD/MS can be used to calculate its static and dynamic properties, among which is per-haps the most fundamental characteristic of dislocations: the Peierls stress.
In the specific case of aluminium, available values of the Peierls stress in the literature show a spread of around two orders of magnitude. To provide some clarity in this topic, Partiideals with the generation and the properties of single edge dislocations. Here a detailed description of a simple method for inserting dislocations into face-centred cubic crystals is introduced and a study of the motion of partial dislocations and its relation with the Peierls stress is carried out.
Following the line of thought used in the case of dislocations, the simulation of GBs also requires as an initial step the generation of a system containing the desired GB. It turns out that, due to the complex order of the atoms in interfaces, the generation of realistic GBs is more delicate than just putting two single crystals together; GBs with equal orientations can have different energies.
To shed some light in this topic, Partiiitreats the subject of GBs at the atomic scale. First, two different methods for the generation of GBs in atomistic sim-ulations are shown, one of which is based on genetic algorithms. And second, artificial neural networks are used for the prediction of GB energies as a func-tion of their misorientafunc-tion.
A few attempts have been made to address the interaction between dislo-cations and a GB, but mainly on a phenomenological basis. This plasticity-interface interaction takes place at the level of individual dislocations, which can be blocked by, absorbed in, transmitted through or nucleated from GBs.
Partivgathers the results obtained and the methods developed in the previ-ous parts, and combines them to perform an analysis of the interaction between
a GB and a single dislocation. This analysis is focused on extracting quantitative information that is useful for larger scale simulation methods
Part I
1
A T O M I S T I C R E L A X A T I O N O F
S Y S T E M S C O N T A I N I N G
P L A S T I C I T Y E L E M E N T S
1Molecular Statics (MS) simulations are widely used for the study of materials.
For example, many properties of dislocations –which play a major role in the deformation behaviour of crystalline materials– were recently calculated using this method. Such properties include core structure [1–4], motion [5–7], Peierls stress [8–11], the interaction with voids [12] and interstitials [13], and more [14–17]. The MS technique is usually carried out by applying a series of discrete deformations to a sample, each being followed by a full relaxation of the system. This relaxation requires the use of a minimization routine to obtain equilibrium after each deformation step.
In most cases the relaxation routine is based on the Conjugate Gradient (CG) method [4,7,10,13,16], and the potential energy of the system is used as the diagnostic parameter to determine if convergence of the relaxation has been reached [18]. Surprisingly, there does not seem to be much discussion about these choices. Although there exist a few articles where minimization methods for Molecular Dynamics (MD) simulations are benchmarked [19, 20], hardly any concern is given to the question of which parameter is best to determine convergence in MS simulations. As the potential energy of a configuration is
1 Adapted from: Echeverri Restrepo S., Sluiter M. H. F. and Thijsse B. J. Atomistic relaxation of systems containing plasticity elements. Computational Materials Science, 73, 154-160 (2013).
relatively easy to compute and must be in a (local) minimum for equilibrium to be reached, it is natural to let the energy values drive the relaxation.
The current chapter will show that applying these common methods without much deliberation may waste considerable computing time, and, even worse, may lead to erroneous results and distorted physics. The efficiency of the CG is compared with that of six other relaxation algorithms and it is shown that CG, together with rapid thermostat-controlled quenching, is an order of magnitude less efficient than two modern methods (FIRE [19] and L-BFGS, to be discussed below).
Additionally, the influence of “plasticity elements” –structural features that play a role in plastic behaviour such as dislocations and grain boundaries– in the relaxation procedure is discussed in terms of the difference of the convergence behaviour of the energy and the average stress, and of the importance of local effects to assure convergence. It is also shown that the choice of the energy by itself as a convergence criterion for structural relaxation is incorrect if one wants to obtain equilibrium also in the atomistic stress field.
Part of this work can be seen as an extension of the results reported in the paper [19] introducing the FIRE method. The current comparison of energy relaxation and stress relaxation is –to the best knowledge of the authors– new to the materials modelling community.
This chapter is organized in the following way: In Section1.1a description of the computations that were used for the analysis is given. Section 1.2presents a benchmark of seven different relaxation routines. Sections 1.3 to1.6 present results on the convergence of the energy and the stress, both locally and glob-ally, leading to conclusions about the behaviour of plasticity elements and the role that they play during relaxation. The relaxation algorithms themselves are summarized in AppendixA.
1.1
computational approach
An MS simulation, in which the response of a system to an imposed deforma-tion is studied, usually consists of the following steps:
1: Generate the desired initial atomic configuration
2: Deform the system incrementally, and impose a constraint to prevent the system from returning to the earlier state
3: Fully relax the system, while maintaining the constraint (some atoms may be constrained in one particular direction only)
4: Return to 2 until the desired full deformation is obtained
Step 3 is the most time consuming task of the algorithm, therefore it is im-portant to choose an efficient method. In this step, the system is relaxed until
the nearest local energy minimum state is reached. If atomic velocities are con-cerned in the minimization algorithm, they will all be essentially zero at the end of the relaxation (0 K temperature). If only atomic positions are concerned in the algorithm, this issue is immaterial; at the end of the deformation it may yet be desirable to set the velocities to zero for subsequent operations.
In this work, the problem of structural relaxation –after the application of a deformation increment in an MS simulation– is analysed using seven iterative minimization routines, divided into two groups. The first group consists of five “dynamic” methods that need to be linked to a molecular dynamics integrator.
Along with solving the system’s equations of motion they minimize the energy by damping the movement of the atoms according to a particular strategy. These five methods are:
• Cooling to 0 K [21]
• Uphill Freezing (UF)
• Uphill Freezing xyz (UFxyz) • Quick-min [20]
• Fast Inertial Relaxation Engine (FIRE) [19]
The second group consists of two “static” methods that minimize the inter-atomic forces by iteratively repositioning the atoms according to an optimiza-tion scheme:
• Conjugate Gradient (CG) [22]
• Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS)
In these static methods the atomic velocities are not considered. The algo-rithms of the seven methods are given in AppendixA.
For all minimization routines applied in this work, the calculation of the forces between the atoms was performed using an interatomic potential for cop-per developed in reference [23] based on the embedded atom method (EAM) [24]. This potential has been shown to give correct values for the stable and unstable stacking fault energy [25], making it an adequate choice to simulate the nucle-ation properties of dislocnucle-ations [26].
For the comparison of the minimization methods, a study was performed using a Face-Centred Cubic (FCC) crystalline system containing an edge dislo-cation and a discretely applied shear stress (System I, Figure i.1). In order to obtain comparable results, the following parameters were chosen to be equal for all the relaxations:
• Initial positions of the atoms.
• Initial velocities set to zero (for the dynamic methods).
• Equations of motion integrated using the velocity-Verlet algorithm
imple-mented in the MD code CAMELION [27].
• Fixed time step of 30 fs for the dynamic methods (except FIRE which has
an automatic time step control limited, in this case, to a maximum of 21.27 fs).
• Periodic boundary conditions along ¯112 and [110], the directions of the
dislocation line and the Burgers vector, respectively.
• Free surfaces in the direction perpendicular to the slip plane 1¯11. • No pressure control imposed during the minimization process (box size
fixed).
After the benchmarking of the different methods was done, two other FCC systems were built to further study the behaviour of the relaxation: one contain-ing a perfect crystal, System II, and one containcontain-ing a grain boundary, System III (Figure i.1), both subject to the same boundary conditions as the system containing the dislocation.
1.1.1 System set-up
A simulation box with dimensions 1.77 × 26.65 × 25.04 nm3containing a single copper FCC crystal was generated with the directions ¯112, 1¯11 and [110] pointing along the global x, y and z axes, respectively. A single edge dislocation was inserted by removing two adjacent (110) half planes [28]. Periodic boundary conditions along the x and z directions were prescribed, corresponding to the direction of the dislocation line and its Burgers vector, respectively. Free surfaces were used along the y direction. The system was composed of 93600 atoms after the insertion of the dislocation.
An initial MD simulation was run at a fixed temperature of 10 K while main-taining the pressure of the system at 0 Pa along the two periodic directions using a Berendsen thermostat and barostat [21]. As expected, the inserted dislocation splits into two partial dislocations separated by a stacking fault ribbon. The sim-ulation was run until equilibrium in the total energy (kinetic + potential) of the system was reached. This equilibrated system at 10 K was subsequently relaxed using each of the seven techniques described above. The same final configura-tion was obtained for all the methods.
For the generation of System I this relaxed configuration was taken as starting point. In order to prepare for applying shear, the first five 1¯11 atomic planes on the top (top atoms) and the bottom (bottom atoms) (Figurei.1) were harmon-ically attached to points in space that will be called “anchor points”. Initially, the anchor points coincide with the bottom and top atoms, so that the harmonic forces are null. The application of a shear strain increment to the system was done by rigidly displacing the top anchor points and atoms by 0.013 42 Å along the [110] direction, which can be considered as the first step of a full shear stress application. After the displacement of the top atoms, the system was relaxed using the seven techniques while keeping the positions of all the anchor points fixed.
A similar approach was followed for the generation and relaxation of Sys-tems II and III, containing a single crystal and a bicrystal with a grain bound-ary, respectively. The single crystal was oriented in the same way as System I, is composed of 94080 atoms, and the dimensions of the simulations box that con-tains them are 1.77 × 26.65 × 25.04 nm3. The grain boundary of System III, is described using the interface-plane scheme [29] as (174), (714), 0 and was gen-erated using the method explained in reference [28]; it is composed of 110376 atoms and the size of the simulation box is 1.77 × 31.13 × 25.15 nm3 . The ap-plied initial shear and the boundary conditions for the relaxation of these two systems were the same as for System I.
1.2
algorithm benchmarking
In order to compare the seven selected algorithms, they were applied to System I, and the convergence of the potential energy of the system per atom, E was monitored. All methods converged to the same final value of energy, but the convergence rate varied widely among them. Here we measure the convergence rate in number of iteration steps. This is somewhat flattering for CG, which –as only method– requires the forces on the atoms to be calculated at least twice per step, because a line search is needed. Note that monitoring the energy of the system and the total magnitude of the forces is almost equivalent. In fact, if the total sum of the forces is equal to zero, classical mechanics ensures that the system potential energy is at its minimum.
The convergence of the seven methods is shown in Figure i.2, and summa-rized in Table i.1, where Nit
¯E is defined as the number of iterations needed to reach the final value within 8.7 × 10−10 eV/atom. This value is a rather arbi-trary choice, based on another system of units, which was found to work well in practice. Note that for systems that are larger than the ones considered here, the atoms far away from the plasticity element will be negligibly influenced by
I II III
x y
z
Figure i.1:Initial configurations for Systems I, II and III to compare the relaxation
meth-ods. Shear strain is applied by displacing the top atoms and top anchor points along the z direction while keeping the bottom atoms and bottom anchor points fixed. System I contains a single edge dislocation, System II is a perfect crystal and System III contains a grain boundary. Atoms are coloured ac-cording to their local symmetry type [30]: light blue represents Face-Centred
Cubic (FCC), yellow Hexagonal Close-Packed (HCP) and dark blue atoms that are neither FCC nor HCP. Bottom and top atoms are highlighted with a darker colour [31].
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 10 100 1000 ∆ E − [10 −7 eV/atom] Steps Cooling to 0 K UF UFxyz Quick−min FIRE CG L−BFGS
Figure i.2:Energy relaxation of System I, containing a dislocation (Figure i.1), after
rigidly displacing the top atoms and top anchor points along the z direction by 0.013 42 Å to generate a shear stress. Steps 10 to 2000. Energy deviation from converged value per atom vs step.
it. It would make sense to apply the 8.7 × 10−10 eV/atom convergence criterion only to the atoms in the region having a similar proximity to the plasticity ele-ment as the systems used here. The fastest two methods are L-BFGS and FIRE, followed by UFxyz, UF, Quick-min, CG and Cooling to 0 K.
1.3
average and local energy convergence
To perform a more detailed analysis of the relaxation behaviour the method FIRE is chosen because of its fast convergence rate and simplicity of implemen-tation. Additionally, to increase the number of examples and to be able to un-derstand the influence of plasticity elements Systems II and III are also subject to relaxation.
Figures i.3a andi.3b show the evolution of the average energy and the stan-dard deviation of the atomic energies, respectively. The first two columns of Tablei.2show the number of steps needed to reach convergence of these two quantities, which are defined to have converged when they reach their final value within 8.7 × 10−10 eV/atom. Figuresi.3a andi.3b show, first, that despite being structurally different, the average energy in the three systems converges
Table i.1:Number of iterations needed for convergence of the average energy (energy
per atom) Nit E
of System I for the seven different minimization methods.
Method Nit¯E Cooling to 0 K 53604 UF 10350 UFxyz 9273 Quick-min 21334 FIRE 460 CG 31914 L-BFGS 458
Table i.2:Number of iterations needed for convergence of the average energy (energy
per atom)Nit¯E, the standard deviation of the atomic energies NitStd.Dev.E and the average shear stressNit
¯τyz of Systems I, II and III. The minimization
method used is FIRE.
System Nit
¯E NitStd.Dev.E Nit¯τyz
I Dislocation 460 400 2560
II Single Crystal 490 530 1180 III Grain Boundary 700 1680 3800
in a similar way and, second, that there does not seem to be a one-to-one rela-tion between the evolurela-tion of the average energy and the standard deviarela-tion of the atomic energies, especially for System III.
The important conclusion can be drawn that reaching convergence in terms of average energy does not necessarily mean that the system is completely re-laxed; Figurei.3b illustrates this: the standard deviation of the atomic energies converges with a rather unstable behaviour. The reason for this is that the aver-age energy does not automatically take into account the full influence of local variations and “outliers” that have not reached equilibrium even after there is no noticeable change in the average energy. Therefore, using the standard devia-tion of the atomic energies as an addidevia-tional convergence indicator to the average may very well be a reliable method to safeguard against premature stopping of convergence algorithms.
−1 0 1 2 3 4 5 6 7 8 10 100 1000 ∆ E − [10 −7 eV/atom] Steps I II III a) −6 −5 −4 −3 −2 −1 0 1 2 3 10 100 1000 ∆ Std. Dev. E [10 −7eV/atom] Steps I II III b) −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 10 100 1000 ∆τ − yz [MPa] Steps I II III c)
Figure i.3:Relaxation of Systems I, II and III. Steps 10 to 2000. Deviation from converged
values. a) Energy per atom (average energy) ¯E vs step. b) Standard deviation of the of the atomic energies (Std.Dev.E) vs step. c) Average shear stress ( ¯τyz)
1.4
local stress convergence
In many types of computations the energy of the system is not the key parame-ter. For example, if one wants to determine the Peierls stress, the main interest lies in the local shear stress on the attached top and bottom atoms (see Figurei.1) after a full relaxation of the system (τtop, τbottom). For the atoms in the top layers of the current configurations, τtop is calculated by summing the [110] components of the forces acting on them –generated by the harmonic forces from being attached to the anchor points–, P
i
Ftopiz , and dividing the result by the area of the 1¯11 plane A(1 ¯11):
τtop = P i Ftopiz A(1 ¯11) (i.1)
Similarly, τbottom is calculated from the forces on the bottom atoms from their anchor points. As said earlier, it is usually assumed in the literature, with-out much discussion, that as a by-product of the minimization of the energy of a system, the stresses will also reach their equilibrium values. However, we will now show that this is not automatically the case.
Figurei.4illustrates the behaviour of τtopand τbottomduring the relaxation of Systems I, II and III using FIRE. As in the case of the average energy, in spite of the presence of a dislocation and a grain boundary in Systems I and III, respectively, the behaviour is very similar in the three systems.
Initially the stress on the top atoms is high, and it gradually decays to a con-stant value, while the bottom atoms have the opposite behaviour. The reason is that, at the beginning of the relaxation, only the top atoms are displaced, gen-erating a region of stress concentration on the plane where they interact with the free atoms. This generates a wave that propagates downward reaching the bottom atoms after about 500 iteration steps, which can be deduced from the increase of −τbottom in Figurei.4. The FIRE algorithm takes care of the nec-essary damping. Furthermore, if compared with Figure i.2, one also sees that τtopand τbottomdo not behave in the same way as the energy throughout the relaxation. This is explained by the fact that here local quantities are measured, namely the stress values at the bottom and top extremes of the system. Addi-tionally, there is an influence from the way in which the strain is applied. This is the issue of the difference in energy and stress convergence just hinted at. It will be addressed in more detail in the following section.
−5 0 5 10 15 20 25 30 100 1000 τyz (local) [MPa] Steps τtop −τbottom I −5 0 5 10 15 20 25 30 100 1000 τyz (local) [MPa] Steps τtop −τbottom II −5 0 5 10 15 20 25 30 100 1000 τyz (local) [MPa] Steps τtop −τbottom III
Figure i.4:Evolution of the stress during relaxation after rigidly displacing the top atoms
and top anchor points along the z direction by 1342 Å to generate a shear stress. Systems I, II and III, contain a dislocation, a single crystal, and a grain boundary, respectively (Figurei.1). The minimization method used is FIRE. Shear stress on the bottom (τbottom) and top atoms (τtop) vs step. Steps
1.5
average stress convergence
Unlike the results obtained in the previous section –where it was shown that the evolution of the stress during a relaxation varies depending on the place where it is measured, and that no apparent influence of the presence of a grain boundary or a dislocation is detected– measurements of the average shear stress during relaxation show that the presence of plasticity elements has an important influence on the behaviour of the stress. Here, the average stress is calculated as the average of the atomic stresses of all the atoms in the system; the atomic stress
σiαβ
is calculated using the virial definition without the kinetic portion [32]:
σiαβ= − 1 2Ωi X j Fijαrijβ (i.2)
In the previous equation Ωi is the atomic volume of atom i, Fijα is the α component of the force on atom i from atom j, rijβ is the β component of the vector from atom i to atom j and α, β ∈ x, y, z.
Figure i.3c shows the deviation to convergence of the average shear stress in the system ( ¯τyz= ¯σyz). Additionally, the third column of Table i.2 shows the number of steps needed to obtain convergence of ¯τyz within a range of 7.2 × 10−8MPa. This value is again a rather arbitrary choice, based on another system of units, which was found to work well in practice.
The convergence of Systems I and III follows a similar trend: first there is a plateau where ¯τyz remains constant and then there is a sudden drop that quickly leads the system towards equilibrium. Although this behaviour might seem puzzling, it can be explained by the discrete way in which the deformation is applied and the presence of plasticity elements: at the beginning of the relax-ation, the top atoms of the systems are rigidly displaced in a discrete motion generating a wave that propagates towards the bottom atoms. Initially, this wave travels through a perfect crystal generating only elastic deformation and, thus, no change in ¯τyz. After approximately 300 relaxation steps, the wave reaches the plasticity element –a dislocation in System I, a grain boundary in System III– generating plastic deformation and, thus, “dissipating” part of the stress and causing the rapid decrease of ¯τyz. It is also interesting to notice that the ori-entational difference between systems I and III does not affect the convergence as much as the difference between the defects that they contain. The behaviour of System II is simpler to explain: since it is a perfect FCC crystal, there is only elastic deformation and, therefore, ¯τyz remains constant (there are only small fluctuations around a constant value).
A comparison between ¯E (Figure i.3a) and ¯τyz (Figure i.3c) makes evident that, although they are both interdependent global quantities –they are both de-termined by the positions of the atoms–, their behaviour during relaxation can be very different. Again, the present results suggest that to control convergence in molecular static algorithms, it is better to use two (or even more) indicators rather than one.
1.6
conclusions
An analysis of the relaxation behaviour of atomistic systems containing plas-ticity elements –grain boundaries and dislocations– after the application of a deformation increment in an MS simulation was performed.
For this, first, seven relaxation methods were compared in terms of energy convergence in a system containing a dislocation and subjected to an externally applied displacement. All the methods eventually converged to the same energy value; the fastest algorithms are FIRE and L-BFGS, being up to approximately 60times faster than the widely used conjugate gradient method.
After selecting FIRE as preferred method due to its fast convergence and sim-plicity of implementation, the relaxation of two other systems –one containing a single crystal and one a grain boundary– was performed. Although no strong influence of the plasticity elements was detected in the convergence of the aver-age energy, the averaver-age stress and the standard deviation of the atomic energies evolved at much slower rates than the average energy. It was thus seen that the energy convergence can be reached long before (local) stresses and (local) structural relaxations reach convergence. This is a potential danger in relaxation practice and should be kept in mind carefully.
It was also noticed that there is not a one-to-one relationship between the convergence of the average energy and that of the standard deviation of the atomic energies. In some cases the standard deviation takes many more itera-tions to converge. This effect was especially strong for the most complex system studied, where a grain boundary is present.
Finally, an analysis of the evolution of the average stress showed that the pres-ence of plasticity elements, combined with the boundary conditions used, can explain its relaxation behaviour. Here no direct relation with the evolution of the energy was seen, and it was also found that many more iterations are needed to obtain stress equilibrium (up to a factor of approximately five compared with the energy). It is thus recommended not to use only the average energy, but to also consider the stress and the standard deviation of the atomic energies as more accurate parameters to determine the convergence of the relaxation of such systems.
B I B L I O G R A P H Y
[1] Duesbery M S, Vitek V and Bowen D K 1973 Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 332 85–111
[2] Fang Q and Wang R 2000 Physical Review B 62 9317–9324
[3] Gröger R, Bailey A and Vitek V 2008 Acta Materialia 56 5401–5411 [4] Monnet G and Terentyev D 2009 Acta Materialia 57 1416–1426
[5] Bulatov V V, Yip S and Argon A S 1995 Philosophical Magazine A 72 453–496 [6] Chaussidon J, Fivel M and Rodney D 2006 Acta Materialia 54 3407–3416 [7] Groh S, Marin E B, Horstemeyer M F and Bammann D J 2009 Modelling and
Simulation in Materials Science and Engineering 17 075009 [8] Sydow B 1999 Computational Materials Science 15 367–379
[9] Olmsted D L, Hardikar K Y and Phillips R 2001 Modelling and Simulation in Materials Science and Engineering 9 215–247
[10] Osetsky Y N and Bacon D J 2003 Modelling and Simulation in Materials Sci-ence and Engineering 11 427–446
[11] Pasianot R C and Moreno-Gobbi A 2004 Physica Status Solidi (B) 241 1261– 1268
[12] Monnet G 2007 Acta Materialia 55 5081–5088
[13] Tapasa K, Osetsky Y and Bacon D 2007 Acta Materialia 55 93–104
[14] Pestman B J, De Hosson J T M, Vitek V and Schapink F W 1991 Philosophical Magazine A 64 951–969
[15] Rao S 1998 Philosophical Magazine A 77 231–256
[16] Bulatov V, Richmond O and Glazov M 1999 Acta Materialia 47 3507–3514 [17] Osetsky Y 2000 Journal of Nuclear Materials 283-287 784–788
[19] Bitzek E, Koskinen P, Gähler F, Moseler M and Gumbsch P 2006 Physical Review Letters 97 170201
[20] Sheppard D, Terrell R and Henkelman G 2008 The Journal of chemical physics 128 134106
[21] Berendsen H J C, Postma J P M, van Gunsteren W F, DiNola A and Haak J R 1984 The Journal of Chemical Physics 81 3684
[22] Press W H, Flannery B P, Teukolsky S A and Vetterling W T 1992 Numerical Recipes in Fortran 77: The Art of Scientific Computing 2nd ed (Cambridge University Press)
[23] Mishin Y, Mehl M, Papaconstantopoulos D, Voter A and Kress J 2001 Phys-ical Review B 63 224106
[24] Daw M and Baskes M 1984 Physical Review B 29 6443–6453
[25] Spearot D, Jacob K and McDowell D 2007 International Journal of Plasticity 23 143–160
[26] van Swygenhoven H, Derlet P M and Frø seth A G 2004 Nature materials 3 399–403
[27] Thijsse B J 2011 Camelion
[28] Echeverri Restrepo S and Thijsse B J 2011 MRS Proceedings 1224 1224–GG05– 03
[29] Randle V 1993 The measurement of grain boundary geometry Electron mi-croscopy in materials science series (Institute of Physics)
[30] Lazi´c I, Klaver P and Thijsse B 2010 Physical Review B 81 1–13
[31] Li J 2003 Modelling and Simulation in Materials Science and Engineering 11 173–177
[32] Allen M P and Tildesley D J 1989 Computer Simulation of Liquids Oxford science publications (Oxford University Press, USA)
Part II
2
G E N E R A T I O N O F
D I S L O C A T I O N S
Most atomistic studies of dislocations require the artificial insertion of disloca-tion(s) into the initial system. To this end several approaches have been pro-posed [1–7].
Due to the fact that the objectives of this work include being able to insert dislocations into bicrystals that can be used to make a connection with a compu-tational method that works at a higher scale –two Dimensional Discrete Disloca-tion plasticity (2D-DD) [8]– the following two requirements need to be fulfilled:
• Only edge dislocations are considered
• The crystal that will contain the dislocations needs to be periodic along the dislocation line−→ξ
Because of its simplicity, a method similar to the one used in references [1,6] is adopted, but modified in such a way that, instead of a pair, only one edge dislocation is generated. The procedure of dislocation insertion can be divided into three stages: removal of planes, application of a displacement field and atomistic relaxation.
2.1
removal of planes
Edge dislocations are line defects characterized by the lack of atomic half-planes in an otherwise perfect crystal. These planes are perpendicular to the Burgers vector (−→b) [9] and contain the dislocation line −→ξ. Depending on the crys-talline structure of a material, a different number of planes have to be removed (or inserted) in order to generate an edge dislocation.
In face centred cubic (FCC) materials the edge dislocation with the lowest energy and, thus, the most common type is−→ξ =h211iand −→b = 12h110i. The reason for this is that the shortest lattice vector is 1
2h110i. The atomic planes perpendicular to−→b have a stacking sequence · · · ABABAB · · · where A and B are the two possible layer positions in this direction [10]. In order to generate a perfect dislocation two half planes (one A and one B) have to be removed (see Figureii.1a-b).
Similarly, in body centred cubic crystals (BCC) materials, where the most common dislocation is defined by −→ξ = h211i and −→b = 12h111i [2, 11], and the stacking sequence in the direction of −→b is · · · ABCABCABC · · · , three half planes need to be removed.
2.2
application of a displacement field
Starting with a perfect crystal, the necessary half-planes are removed (see ii.1 a-b) and a displacement field predicted by the elastic isotropic theory of dislo-cations, centred in the core of the dislocation, is applied to all the atoms in the system (see Figure ii.1b-c) in order to get a first approximation to realistic atomic positions around the dislocation. This displacement field is given by the following equations [11]: uy = − b 2π " 1 − 2ν 4 (1 − ν)ln z2+ y2+ z 2− y2 4 (1 − ν) z2+ y2 # (ii.1) uz= b 2π " tan−1y z + zy 2 (1 − ν) z2+ y2 # (ii.2)
where uα is the displacement to be applied to each atom in the α direction; bis the magnitude of the Burgers vector (in this case−→b = 12h110i); y and z are the original coordinates of the atoms, and ν is the Poisson’s ratio of the material.
Removal of planes Application of a Displacement field Atomistic relaxation a) b) d) c) x y z
Figure ii.1:Generation of a single dislocation for MD simulations. a) Perfect FCC crystal
with h211i pointing towards −→x, h111i towards −→y and h110i towards −→z. a-b) Removal of two h110i half-planes. b-c) Application of the displacement field predicted by isotropic elastic theory of dislocations [11]. c-d) Relaxation of
It is important to note here that the crystal has to be oriented in such a way that −
→x is parallel to−→ξ, −→y is perpendicular to the slip plane and −→z is parallel to−→b.
2.3
atomistic relaxation
The final step needed for the generation of a dislocation for atomistic simula-tions, consists on the relaxation of the system (Figure ii.1d). For this, several methods can be used (see Appendix A) and the idea is to allow the atoms to re-accommodate to find the positions that minimize the energy of the system.
In the case of aluminium, copper and FCC materials in general, after relax-ation, the inserted dislocation splits into two partial dislocations separated by a stacking fault ribbon. This is consistent with experimental results and the reac-tion can be written in terms of the burgers vector of the initial dislocareac-tion (−b→1) and that of the two generated partials (−b→2,−b→3):
−→ b1→ −→ b2+ −→ b3 1 2h110i → 1 6h211i + 1 6121 (ii.3) Figure ii.1d shows the detail of the separation of a perfect dislocation into two partials in the case of aluminium modelled with the interatomic potential developed by Mishin et al. [12].
2.4
conclusions
In the present chapter, a simple and versatile method to insert edge disloca-tions into FCC crystals was presented. This method will be utilized throughout the remainder of this work as a point of departure for the analysis of single dislocations and of their interaction with grain boundaries.
3
P E I E R L S S T R E S S
Apart from a few special cases, the plastic behaviour of crystalline materials is determined by the presence and movement of dislocations. Atomistic sim-ulations have been widely used to study their behaviour at time and length scales that are not (yet) achievable using experimental methods. Thanks to these simulations, new insights have been obtained about their structure [13–19], mo-tion [1,2,4,6,14,15,18,20–24], interaction with point [2,13,25–30] and planar defects [5,31–49], and more.
In previous decades, several studies of the motion of dislocations across Peierls barriers have appeared in the literature. Two representative examples are the work of Benoit et al. [50], in which experimental internal friction data were explained in terms of kink pair formation and on the basis of an assumed pres-ence of interstitials in the slip plane between the two partials; and the paper by Schoeck and Püschl [51], which explains high and low values of τpby coupled partials having an idealized separation distance of an integer or half-integer multiple, respectively, of the period of the Peierls potential. The idealized sepa-ration distance is the equilibrium sepasepa-ration distance in the hypothetical case of a vanishing Peierls potential. The analysis in [51] is executed for a screw disloca-tion in a simple cosine potential and assuming linear elastic interacdisloca-tion between the partials.
With the aim to enhance the analytical approach of Schoeck and Püschl to more practical conditions, the work presented here is done in full atomistic
de-tail. For this, Molecular Statics (MS) simulations are performed to calculate τp and analyse the behaviour of a single edge dislocation in aluminium modelled with the Embedded Atom Method (EAM) [52] potential developed by Mishin et al. (MFMP) [12], which has been widely used for the simulation of disloca-tions [40,53–55]. Sample conditions in the present work are such that no kink formation was observed, nor were interstitials or other point defects detected in the stacking fault strip, dealing thus with dislocations in their purest form.
The results are compared with another EAM potential for aluminium devel-oped by Ercolessi and Adams (EA) [56], for which more information about the values of τp that it produces is available in the literature. The comparison be-tween the values of τp for the two potentials yields a difference of over two orders of magnitude, which might seem surprising in view of the fact that both potentials were designed to model aluminium in a realistic way. This difference will be explained based on a different movement of the partial dislocations.
The present chapter is organized in the following way: Section 3.1 summa-rizes the existing methods for the calculation of τp existing in the literature; Section 3.2presents the process used in this work for the insertion of an edge dislocation in an aluminium crystal and the procedure used for the application of a shearing strain to force the dislocation to move. In Section3.3the method-ology used for the selection of an appropriate system size and for the measure-ment of τp is explained. Finally, in Section3.4, the obtained values for τp and their mismatch with each other and with values available in the literature is accounted for via a direct analysis of the movement of the partial dislocations
3.1
methods for the calculation of the peierls
s t r e s s
Various approaches have been proposed to measure τp using atomistic simu-lations, each one of them having its own advantages and disadvantages. Here, four different types of methods are mentioned and briefly explained with the aim of justifying the choice made. More detailed reviews can be found in refer-ences [2] and [57].
The first and simplest type (Group I) [24,58,59] consists of inserting a single dislocation in an otherwise perfect crystal and selecting all the atoms outside a given radius from the dislocation line. All the selected atoms are then displaced using an elastic field that mimics the long range strains generated by the dislo-cation and are then fixed to their new positions. Finally the movement of the dislocation is enforced by the application of a uniform strain field to the system. In this approach, periodicity is only maintained along the dislocation line.
The main disadvantage of this group of methods is that when the dislocation starts to move, the applied elastic field on the atoms that are away from the dislocation is no longer correct, which can lead to an overestimation of the value of τp [57].
The second type (Group II) [60, 61] is similar to the previous one but with more sophisticated ways of dealing with the boundary conditions. Here, a dis-location is inserted into a perfect crystal while enforcing periodicity along the dislocation line. Then, the influence of the boundary forces is dynamically elimi-nated by by determining their contribution and correcting for their effect [57,62] or by using flexible boundary conditions [63].
The main point in favour of these methods is that the size of the system can be considerably reduced while still obtaining a correct core structure and value of τp, but the choice of a small system hinders the simulation of long range interactions [2]. Another method that can also be included in this group is the quasi-continuum method, which takes into account the boundary forces by linking the system to a continuum model [64].
The third group (Group III) contains the methods that use periodic simulation boxes along the direction of the Burgers vector and of dislocation line [2,24,29, 30,65]. For these methods, a dislocation is inserted and the first layers of atoms at the free boundaries (perpendicular to the glide plane) are kept fixed at their positions. It is important to note that the forces generated by these immobilized atoms do not have a component in the glide direction [57].
A disadvantage of this group of methods is that the dislocation “feels” the presence of its periodic images. Nevertheless, if the system is big enough in the glide direction, the core size can be correctly represented [63] and, since there is symmetry in the periodic images, the net force acting on the dislocation caused by the periodicity cancels out.
The last group of methods (Group IV) uses a fully periodic system [6,66]. For this group, the selection of the boundary conditions obliges the insertion of dislocations to be done in pairs, this being the cause of one of its major flaws: all the quantities that can be measured for one of the dislocations are corrupted by the stress (strain) field generated by the presence of the other one.
Due to the simplicity of the methods and the fact that it has been successfully used to determine quantitative properties of dislocations [2], a methodology that can be classified as belonging to Group III is adopted here, taking spe-cial care of the influence of the boundary conditions and making sure that the relaxation after every deformation step is fully achieved [67].
a)
b)
c)
d)
z
y
x
Figure ii.2:Construction of an edge dislocation in a Face-Centred Cubic (FCC) crystal. a)
Removal of two h110i half planes. b) Relaxation of the system using Molec-ular Dynamics (MD). c) Insertion of two semi-rigid slabs of atoms. d) MD relaxation of the system keeping the atoms in the inserted slabs attached to harmonic springs to straighten the system.
3.2
simulation set-up
For the set-up of the simulation, an FCC crystal is generated in such a way that the directions ¯112, 1¯11 and [110] point towards the global −→x, −→y and −→z axis, respectively. Periodic boundary conditions are imposed along −→x and −→z, while the boundaries along −→y are left as free surfaces. In order to insert a single edge dislocation, two (110) half planes are removed (Figureii.2a) and a displacement field –corresponding to an edge dislocation– calculated using isotropic elastic-ity theory [10] is applied to all the atoms. The system is subsequently relaxed using the method FIRE [68] while maintaining the pressure set to zero along the directions of the dislocation line (−→x) and of the Burgers vector (−→z).
As a result of the relaxation the inserted dislocation splits into two Shock-ley partials separated by a stacking fault ribbon (Figureii.3). Additionally, the whole system becomes bent due to the difference in the number of planes in the upper and lower part (Figureii.2b). This bending causes the sliding plane of the dislocation to be non-straight and, therefore, has an influence on the movement of the dislocation and on its Peierls stress. In order to straighten the dislocated crystal, two slabs of atoms are added (Figure ii.2c): one to the top and one to the bottom along the non periodic direction (−→y). These slabs contain atoms at-tached by harmonic springs to “anchor points" which are located at the same initial positions of the newly inserted atoms in the slabs. All the system is then further relaxed, becoming almost perfectly straight (Figureii.2d). The y position of the inserted slabs is set in such a way that the average force measured by the springs in the −→y direction is zero.
x y z
Figure ii.3:Relaxed 12h110i {111} edge dislocation –using the MFMP potential– in an
otherwise perfect FCC crystal. After relaxation the dislocation splits in two partials separated by a stacking fault. Grey atoms are FCC, lighter atoms Hexagonal Close-Packed (HCP) and darker atoms are neither FCC nor HCP. Only part of the configuration is shown.
3.3
measurement of the peierls stress
The shear stress is measured in the bottom and top atoms –τbottom[110] and τtop[110]respectively–, using the following relation:
τtop= P Ftop[110] A(1 ¯11) , τbottom= P Fbottom[110] A(1 ¯11) (ii.4)
where P Ftop[110] and P Fbottom[110] are the sums of the components of the forces exerted by the harmonic springs in the direction [110] on the top and bottom atoms, respectively, and A(1 ¯11) is the area of the 1¯11 plane.
After checking that the initial configuration is fully relaxed, the anchor points attached to the top atoms are incrementally displaced in the direction of the Burgers vector [110] –generating a shearing strain (γxz)– in order to force the dislocation to move. This displacements are done by increments of 2.14 × 10−4Å, each of them being followed by a full energy minimization of the system using FIRE [68] and the MFMP potential. The relaxation is only stopped after the values of τbottom[110] and τtop[110] converge to the same magnitude. This deformation-relaxation process is repeated several times in order to force the dislocation to overcome the Peierls barrier.
Table ii.1:Values of the Peierls stress (τp)determined for different system sizes. MFMP potential. System Size Å3 Atoms τp[MPa] 1 7.39 × 142.97 × 257.52 318 784 0.0130 2 7.39 × 234.50 × 347.19 703 088 0.0151 3 7.39 × 294.58 × 409.08 1 040 056 0.0148
To verify that there is no influence of the boundary conditions the calculation of τp is done for three different system sizes, as shown in Table ii.1. Note that since the length of the system along the dislocation line has no influence on the value of τp, it is kept constant.
In Figures ii.4a, ii.4b, and ii.4c the results of the simulations for the three systems are shown. In this work, τp is calculated as the maximum shear stress (maximum absolute value of τbottom[110]and τtop[110]). The values are shown in Table ii.1. Based on these values, System 2 –for which τp = 0.0151 MPa– is chosen as a sufficiently large system to study in more detail the evolution of the dislocation as it overcomes the Peierls barrier. Note that although the displacement is applied through the anchor points, the strain is measured using the positions of the attached atoms.
Using the same method and system size, τp is also determined for the alu-minium EA potential. A plot of the shear stress vs the applied shear strain is shown in the lower part of Figureii.6b. Note that, despite the fact that the only difference between this simulation and the one used to generate Figure ii.4b (repeated in extended form in the lower part of Figureii.6a) is the interatomic potential, The shape of the stress curves is not the same. Additionally, for the EA potential the calculated value of τpis found to be 2.24 MPa, more than two orders of magnitude larger than for MFMP.
The results obtained with the EA potential are of the same order of magni-tude as other MD calculations in the literature using this same potential; they also agree with recent calculations done with orbital-free density functional the-ory [24], see Tableii.2.
The reason for the large difference between the values of τp for the two po-tentials will be discussed in the following section.
3.4
analysis
Tableii.2presents a a compilation of experimental and simulated values of τp available in the literature, where the wide spread in the data can be seen. For
−0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0 2 4 6 8 10 12 14 16
Shear Stress [MPa]
Shear Strain (×10−6) System 1 a) −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0 1 2 3 4 5 6 7
Shear Stress [MPa]
Shear Strain (×10−6) System 2 b) −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0 1 2 3 4 5
Shear Stress [MPa]
Shear Strain (×10−6) System 3
c)
Figure ii.4:Shear stress vs. shear strain strain curves –obtained using the MFMP
potential– for systems 1, 2 and 3 containing a dislocation. The only differ-ence between the systems is their size (see Tableii.1).
Table ii.2:Overview of values of the Peierls stress (τp)reported in the literature.
Author Method τp[MPa] Type
Current Work EAMEAM (EAa(MFMPc) b) 0.01512.24 EdgeEdge Fang and Wang [16] EAM (EA) ≈ 2.4 Edge Fantozzi et al. [69] Internal Friction ≈ 3.2 Screw Kosugi and Kino [70] Internal Friction ≈ 0.96 Screw Olmsted et al. [57] EAM (EA) ≈ 1.914-18 EdgeScrew
Shin and Carter [24] OFDFTd
6.2 Screw 9.9 Screw 355 Screw 1.6 Edge Lu et al. [71] PN-based 3.2, 24 Edge 256, 88 Screw
EAM (EA) 82 Screw
Srinivasan et al. [23] EAM (EA)EAM (BAMe) 113, 225 EdgeEdge Takeuchi [72] Internal FrictionMechanical test ≈ 32≈ 0.3 ScrewScrew Wang [73] Theoretical model ≈ 0.81 Edge Wang and Fang [20] EAM (EA) ≈ 3.2 Edge aEAM, embedded atom method [52]
bMFMP, interatomic potential developed by Mishin et al. [74] cEA, interatomic potential developed by Ercolessi and Adams [56] dOFDFT, orbital-free density functional theory [24]
eBAM, interatomic potential developed by Baskes et al. [75]
simulations, this spread appears to be especially big when different interatomic potentials are used. The results obtained in the current work and those pre-sented in Tableii.2show that the variations are of approximately two orders of magnitude.
In order to understand this spread the most natural thing to begin with is to look at the unrelaxed generalized stacking fault curves [76] of the potentials used, since their maximum slope is a measure of the theoretical shear strength of the lattice [77] and, indirectly, of the Peierls barrier [23,78]. These curves are presented in Figureii.5for the MFMP and the EA potentials. The ratio of their maximum slopes is 1.39, not at all large enough to explain such a big ratio in τp. In addition, the slope of the curve is higher for MFMP, which does not match with the fact that in the current work τpis higher for EA. A similar conclusion
0 50 100 150 200 250 0 0.5 1 1.5 2 Energy/Area [mJ/m 2 ] Shear Displacement [Aº ] MFMP EA
Figure ii.5:Unrelaxed generalized stacking fault curves for the MFMP and EA potentials.
Note that the maximum slopes of these curves, which are measures of the Peierls stress, do not differ enough to account for the difference of two orders of magnitude obtained for τp.
can be drawn from the comparison between EA and the potential developed by Baskes et al. (BAM) [23,75].
A second possible explanation is presented in reference [23], where τp is calculated with the BAM and the EA potentials. In their work the difference in τp is explained by the fact that the ground state of the dislocations is different: for EA the dislocation spontaneously splits into two Shockley partials, while in the case of BAM the dislocation core remains compact. The authors argue that this compactness generates an increase of τp. This hypothesis of the increase of τp with decreasing dislocation core width is supported by the classical Peierls-Nabarro model, τp= µexp −4πζ b , (ii.5)
where µ is the shear modulus (≈ 32 GPa) [12], ζ the core radius and b the magnitude of the Burgers vector. Nevertheless, in the current case, neither of the potentials produce a compact dislocation core; the initial separations be-tween the partial dislocations, i.e. in the presence of the Peierls potentials, are 9.1 Å and 16.4 Å for MFMP and EA, respectively (coincident with the data from
reference [23]). Here, the position of the cores of the partial dislocations are ap-proximated by the coordinates of the atoms that have the highest tensile stresses along −→z. More importantly, contrary to the predictions of the Peierls-Nabarro model, the dislocation with a wider separation of the partial dislocations (EA) produces a higher value for τp.
Since the two previous arguments fail to explain the difference in the relative order and of the two orders of magnitude between the values of τp for the MFMP and EA potentials, a detailed analysis of the initial equilibrium positions and the motion of each of the partial dislocations was carried out.
As mentioned earlier, the initial equilibrium separation of the partial dislo-cations for the EA potential is 16.4 Å and there are ten planes between the ref-erence planes that are used to approximate the position of the partials. This gives an average distance of approximately 1.49 Å along the ¯112 direction for the planes in the tensile region generated by the dislocation within the stack-ing fault, which is bigger than the equilibrium atomic bulk separation of 1.43 Å along the same direction. This means that the deformation is not only concen-trated in the dislocation cores of the partials, but it spreads throughout the width of the stacking fault. A similar analysis for the MFMP potential yields an initial equilibrium separation of 9.1 Å, five planes between the reference planes and an average atomic separation of approximately 1.52 Å along ¯112 within the stacking fault.
A plot of the displacements (u) of the leading (ul)and trailing (ut)partial dislocations, and their separation (dl−t)as a function of the applied shear strain is presented in the top parts of Figures ii.6a and ii.6b for the MFMP and EA potentials, respectively.
An schematic representation of the movement of the two partial dislocations is shown in Figure ii.7b, where the partial dislocations are drawn as point masses, the staking fault between them as a spring and the Peierls potential as a sinusoidal path. For the EA potential, the initial separation of the partial dislocations remains constant as the dislocation travels through the Peierls land-scape in a discontinuous stick/slip mode, making “jumps” of approximately 1.43 Å, following the initial build-up of the necessary stress bias (2.17 MPa), seen in Figureii.6b. This is the situation that in [51] is called “strong coupling” with an idealized separation distance of the partials being an integer multiple of the Peierls period. The hypothetical possibility that the two partials are fully uncoupled and move simultaneously simply because they experience the same shear stress level at the same moment can be discarded. The presence of the stacking fault strip between the partials, with its associated energy, is already proof by itself for the existence of coupling.
The situation with the MFMP potential is more complex. As shown in Fig-ure ii.6a the separation between the partials does not remain constant, in
con-−0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0 2 4 6 8 10 12 14 16 18
Shear Stress [MPa]
Shear Strain (×10−6) MFMP τxz 0 2 4 6 8 10 Displacement [A º] dl−t ul ut a) 0 0.5 1 1.5 2 0 20 40 60 80 100
Shear Stress [MPa]
Shear Strain (×10−6) EA τxz 0 2 4 6 8 10 12 14 16 18 Displacement [A º] dl−t ul ut b)
Figure ii.6:Upper part of the figures: displacement of the leading (ul)and trailing (ut)
partial dislocations and the distance between them (dl−t)as a function of the
applied shear strain. Lower part of the figures: shear stress vs. shear strain as the dislocation moves through the Peierls landscape. Results presented for the MFMP (a) and EA (b) potentials.
i) a) γxz ii) iii) iv) v) b) γxz
Figure ii.7:Schematic model of the movement of two partial dislocations linked by a
stacking fault represented by two point masses and a spring, respectively, in an FCC material. The applied shear strain (γxz)points to the right. In the
case of the MFMP potential (a) the two partial dislocations do not move at the same time, and their motion is assisted by the compression/extension of the stacking fault. For the EA potential (b), the two partials move in phase and the distance between them remains constant. Stages (i) to (v) are explained in the main text. Note that, although it is not shown in this schematic view, the two potentials produce a different separation of the partial dislocations.
trast with the EA potential: the leading and trailing partials do not move in phase. The trailing partial moves first and it is then followed by the leading par-tial, making “jumps” of approximately 1.52 Å. Figure ii.7a shows a schematic representation of their movement as they travel to the next stable position in the Peierls potential. Again one should note the difference between the continuous motion displayed in Figure ii.7a and the stepwise “frozen” motion displayed in Figureii.6a. In fact, the behaviour in Figureii.6a corresponds to the time se-quence (i), (iii), (v) of Figureii.7a, repeated indefinitely. It is again the “strong coupling” case [51], but this time with the idealized separation distance of the partials being a half integer multiple of the Peierls period.
For the MFMP case, the behaviour of the partial dislocations and the stacking fault between them can be explained in the following manner (see Figureii.7a):
• Initially the leading and trailing partials are in their equilibrium positions
(i). The stacking fault can be thought of as a pre-tensed spring.
• As the magnitude of the applied shear displacement increases, both
par-tials start feeling its effect. Since the trailing partial is at the same time being “pulled” by the stacking fault, it climbs the Peierls potential hill (ii) and reaches the next valley, while the leading partial stays trapped (iii). At this point, the partial dislocations find a new equilibrium, one that is different from (i). Here, the stacking fault acts like a compressed spring. In fact, this compressed spring exerts so much force on the leading par-tial that the applied stress needed to move the leading parpar-tial forward, becomes negative (see Figure ii.6a).
• As the magnitude of the shear displacement continues to increase, the
leading partial climbs over the next Peierls hill (iv) and subsequently goes down the valley (v). This happens because it is being “pushed” simultane-ously by the stacking fault. A new equilibrium position equivalent to the initial configuration (i) is reached. As said, for both potentials, the only equilibrium stages are those where both partial dislocations sit in a valley; the other stages are transient and are only shown for the sake of clarity. Note that in this model the energy needed to perform the caterpillar-like mo-tion (Figure ii.7a) is lower than that needed to move the two point masses in phase (Figure ii.7b). This phenomenon where the stacking fault acts like a tensed-compressed spring changing the dynamics of the motion of the disloca-tion explains to a large extent the difference in the calculated values of τpusing different potentials.
This analysis is different from that of reference [50], where it is proposed that the size of the stacking fault can be modified by the presence of point defects along the path of motion of a dislocation, and that the stress needed to move
a dislocation segment may in this way be reduced almost to zero [79]. It is also different to the analysis given in references [72], [23] and [24], where the compactness of the dislocation core takes all the credit for the so called Peierls stress controversy[23] generated by the wide spread in the values of τp for FCC metals reported in the literature.
Finally, the current analysis is set-up as the realistic atomistic version of the analytic analysis pioneered in reference [51]. It is observed here, for the two stud-ied potentials, two widely different cases of partial dislocation movement under applied stress, which agree with the integer and half-integer period separation cases of the analytic approach. The values of the Peierls stresses themselves vary widely: no less than two orders of magnitude.
3.5
conclusions
Simulations of dissociated single edge dislocations in Face-Centred Cubic (FCC) aluminium modelled with two different interatomic potentials using molecular statics were performed and analysed in terms of their Peierls stress τp and of the individual behaviour of their composing partial dislocations.
For the calculations, a single edge dislocation was inserted in an otherwise perfect FCC crystal by removing two half-planes. After letting the inserted dis-location relax and ensuring that the glide plane of the disdis-location was straight, incremental displacements were applied to subsets of atoms to impose a shear strain, and force the dislocation to move and overcome the Peierls barrier.
Values of 0.0151 MPa and 2.24 MPa were obtained for τp using the inter-atomic potentials developed by developed by Mishin et al. (MFMP) [12] and Ercolessi and Adams (EA) [56], respectively. This difference of two orders of magnitude –and the difference with other reported values in the literature– can-not be explained by the classical Peierls-Nabarro model, nor by the differences in the stacking fault energy curves or by the differences in the compactness of the core of the dislocations generated with the used potentials.
The explanation of such differences in the values of τpwas pioneered in the analytical work of Schoeck and Püschl [51]. Here, in full atomistic resolution, the explanation is expanded via a detailed analysis of the motion of the partial dislocations. In the two cases analysed, τpturns out to be very small –when the stacking fault compression and tension almost perfectly match the energy paths of the two partials– or in the 2 MPa range –when the stacking fault operates as a rigid bond between the two partials–.
It is concluded that depending on the idealized width of the stacking fault that separates the partials, the full edge dislocation may behave very differently, and result in a wide spread in the values of τp.
B I B L I O G R A P H Y
[1] Chang J, Cai W, Bulatov V V and Yip S 2002 Computational Materials Science 23 111–115
[2] Osetsky Y N and Bacon D J 2003 Modelling and Simulation in Materials Sci-ence and Engineering 11 427–446
[3] Bulatov V V and Wei C 2006 Computer simulations of dislocations Oxford Series on Materials Modelling (Oxford: Oxford Univ. Press)
[4] Liu X, Golubov S, Woo C and Huang H 2004 Materials Science and Engineer-ing A 365 96–100
[5] Jin Z H, Gumbsch P, Albe K, Ma E, Lu K, Gleiter H and Hahn H 2008 Acta Materialia 56 1126–1135
[6] Chang J, Cai W, Bulatov V and Yip S 2001 Materials Science and Engineering A 309-310 160–163
[7] Bulatov V V 2008 Crystal Plasticity from Dislocation Dynamics Materials Issues for Generation IV SystemsNATO Science for Peace and Security Series B: Physics and Biophysics ed Ghetta V, Gorse D, Mazière D and Pontikis V (Springer Netherlands) pp 275–284
[8] van der Giessen E and Needleman A 1995 Modelling and Simulation in Ma-terials Science and Engineering 3 689
[9] Frank F C 1951 Philosophical Magazine Series 7 42 809–819
[10] Hull D and Bacon D J 2011 Introduction to Dislocations 5th ed Materials science and technology (Oxford: Elsevier Science)
[11] Hirth J P and Lothe J 1982 Theory of dislocations 2nd ed (New York: Wiley) [12] Mishin Y, Farkas D, Mehl M J and Papaconstantopoulos D A 2011 MRS
Proceedings 538 535
[13] Foiles S and Medlin D 2001 Materials Science and Engineering: A 319-321 102–106
[14] Monnet G and Terentyev D 2009 Acta Materialia 57 1416–1426 [15] Gröger R, Bailey A and Vitek V 2008 Acta Materialia 56 5401–5411
