Structure and mobility of the 12<111>{112} edge dislocation in BCC iron studied by molecular dynamics

Available online at www.sciencedirect.com

Acta Materialia 57 (2009) 1416–1426 www.elsevier.com/locate/actamat

Structure and mobility of the 12 < 1 1 1 > f1 1 2g edge dislocation in BCC iron studied by molecular dynamics G. Monnet a,*, D. Terentyev b b

a EDF-R&D, MMC, Aveneue des Renardie`res, 77818 Moret-sur Loing, France SCK-CEN, Nuclear Material Science Institute, Boeretang 200, B-2400 Mol, Belgium

Received 13 August 2008; received in revised form 12 November 2008; accepted 13 November 2008 Available online 10 January 2009

Abstract In this paper, we carried out atomistic calculations to investigate in detail the core structure and motion mechanism of the < 11 1 > f1 1 2g edge dislocation in a-iron. First, molecular statics simulations are used to characterise the dislocation-core structure in the framework of the Peierls–Nabarro model. It is shown that the accommodation of the distortion due to the insertion of the extra half-planes is not equivalent in the planes above and below the dislocation slip plane and that the relative atomic-displacement profile in the dislocation-core region is asymmetrical. Then, molecular dynamics simulations are used to study the mechanism of the dislocation motion at different temperatures. At low temperature, the dislocation is found to move by nucleation and propagation of kink-pairs along its line. Independently of temperature, when loading is performed in the twinning direction, the critical stress is found to be lower than the one corresponding to the antitwinning loading direction. Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

1 2

Keywords: Iron; Edge dislocation; Kink pair; Molecular dynamics

1. Introduction The activity of slip systems, and specially the < 1 1 1 > f1 1 2g type, is a thorny question on plastic deformation of Body-Centred Cubic (BCC) materials [1,2]. In BCC materials, the {1 1 2} planes are often associated with the {1 1 0} planes to accommodate plastic deformation [3,4], in contrast to Faced-Centred Cubic (FCC) alloys that exhibit massive slip on only one crystallographic slip system. This concurrent activation of slip systems of different crystallographic nature leads to a number of consequences. (i) The respective Peierls stresses are a priori different. (ii) The thermal activation does not proceed with the same intensity in both systems. The spread of activation domains therefore depends on temperature. (iii) Impurities, especially interstitial atoms, may generate different pinning forces on dislocations belonging to different slip systems.

1 2

*

Corresponding author. Tel.: +33 1 60736473; fax: +33 1 60736889. E-mail address: [email protected] (G. Monnet).

(iv) Beyond the question of activation, the determination of the cross slip parameters between the two slip systems is a significant challenge in the understanding of deformation mechanisms. (v) The possible activation of the family of the 12 < 1 1 1 > f1 1 2g slip systems increases the number of slip systems to 24, which is twice the number of slip systems in the FCC lattice. This rises the question of the selection of active slip systems at a given loading state. All these features make the determination of the orientation domains for the activation of each slip system complex and difficult to address in experiment, especially when the critical stress depends on the orientation of load. In addition, the unavoidable presence of impurities is expected to modify or even to prevent completely the activation of some slip systems. Furthermore, the lattice friction in the 1 < 1 1 1 > f1 1 2g slip system is known to be dependent 2 on the loading direction. When the stress provokes shear of the dislocation core in the Twining Direction (TD), the critical stress is always found to be less than the critical stress in the Antitwining Direction (AD) [5,6]. The origin

1359-6454/$34.00 Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2008.11.030

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

of this asymmetry has been implicitly ascribed to the asymmetry of the staking fault energy along the ½1 1  1 direction in the (1 1 2) plane, see for example [7]. Taking into account difficulties in experimental testing of BCC single crystals, one may see the reason behind the numerous uncertainties concerning the mechanical behaviour of BCC materials. Nowadays, a variety of widely accepted computer simulation methods exist to tackle physical problems related to phenomena occurring at a space- and time-scales that are not available for experimental assessment. The so-called multi-scale modelling approach begins with the description of fine atomic structures of defects and their properties at the nano scale and continues with plasticity and fracture mechanics treatment at the mesoscopic scale. The latter methods, called dislocation dynamics (DD), deal with the dislocation arrays and require a precise knowledge about mobility of dislocations and mechanisms involved in their interactions. In particular, it is necessary to have information about the mobility of all dislocation characters in dif1 > f0 1 1g ferent slip systems. Investigations of the 12 < 1 1  slip system have already been reported [12,8]. Effort should 1 > f1 1 2g slip sysnow be put to characterize the 12 < 1 1  tem. So, we dedicate this work to the study of the 1 < 11 1 > f1 1 2g edge dislocation, whereas a forthcoming 2 paper will be devoted to the study of the screw dislocation mobility in the {1 1 2} plane. Here, we apply atomistic simulations, in particular molecular dynamics, to characterize the dislocation motion mechanisms and evaluate the corresponding dislocation mobility. From the point of view of atomic-scale modelling, only a few investigations were devoted to study the 1 > f1 1 2g slip system [7,9,10], properties of the 12 < 1 1  where the asymmetry of the edge dislocation-core structure and the asymmetric effect of the applied stress on this structure were revealed. However, results reported in [9,7] were obtained using molecular static (MS) technique employing pair potentials for BCC materials, not specifically designed for iron. The results were found to depend strongly on the type (or version) of the potential used in those simulations, thereby questioning the reliability of the values for the Peierls stress reported in that pioneering work. Recently, a new empirical many body interatomic potential for bcc Fe, fitted to a number of properties obtained from ab initio calculations, has been produced [11]. This potential was also found to reproduce the compact core structure of the ½<1 1 1> screw dislocation in BCC iron as observed in first principle simulations [12]. Therefore we rely on the accuracy and credibility of this potential, which we apply here to describe the properties of the ½<1 1 1 > {1 1 2} edge dislocation in a-iron. First, we describe the core structure of the edge dislocations in the (1 1 2) and (1 1 0) slip planes and compare them in the framework of the Peierls–Nabarro model. In the subsequent section, the stress–strain curves for the loading in the TD and AD at 0 K and at finite temperatures are presented. Then, the kink-pair mechanism controlling the motion of the dislocation at low temperature and the

1417

related thermal activation process are analysed and discussed. Finally, we draw some conclusions about the activation of the 12 < 1 1 1 > f1 1 2g slip system in a-iron. 2. Atomic model of the edge dislocation The simulation model developed by Osetsky and Bacon [8] to study the edge dislocation in the (0 1 1) slip plane was modified to create the edge dislocation in the (1 1 2) plane. In what follows, we shall call the former (1 1 0) dislocation and the latter (1 1 2) dislocation, skipping the ½<1 1 1> specification. Both dislocations are created with a Burgers vector b ¼ 12 ½1 1 1. In the case of the (1 1 2) dislocation, the principal axes x, y and z of the simulated BCC crystal were oriented along the ½1 1 1; ½1 1 0 and ½1 1 2 directions, respectively, as shown in Fig. 1. Periodic boundary conditions were applied along x and y directions [8]. The box was divided into three parts along the z-axis. The upper and lower parts consisted of several atomic planes in which atoms were rigidly fixed in their original positions (block F) or displaced according to the applied loading conditions (block D), see Fig. 1. Atoms in the inner region were free to move (block M). A glide force on the dislocation was generated by the displacement of the block D in the x direction (AD) or in the x direction (TD), which corresponds to simple shear strain c. The corresponding resolved shear stress induced by the applied deformation was calculated as s = Fx/Axy, where Fx is the total force projected in the x direction and computed from all atoms in block D and Axy is the xy cross-section area of the box. As explained in the original work where the currently applied model was suggested [8], the dislocation interacts with its image via imposed periodic boundaries along the direction of the Burgers vector. Hence a periodic array of dislocations, whose density is determined by the size of the crystal, is treated in the present simulations. In this array, each dislocation interacts with its two images (located in the two opposite sides of the crystal) and the net force due to this interaction is zero. It is, however, to be noted that the core structure of the dislocation and/or Peierls stress can be affected by the self-interaction of the

a

b

z = [11 2]

[001]

D

z = [11 2]

x = [111]

x = [111] M F Twining

Antitwining

Fig. 1. (a) MD simulation box consisting of M: freely mobile atoms, F: fixed atoms and D: atoms subject to imposed displacement; (b) a sketch of atomic positions viewed in the y direction. Squares and circles denote atoms belonging to the two neighbouring ð110Þ planes. Two different loading conditions are shown: twinning and antitwinning directions.

1418

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

dislocation with its image if the dislocation density is too high. Thus, we have performed a set of test calculations to study the sensitivity of the Peierls stress and equilibrium (i.e. unloaded at 0 K) core structure to the size of the crystal. It has been found that they both become essentially size-independent if the crystal size along x and z direction exceeds 12 nm. Static simulations were performed in order to compute the structure, elastic energy and Peierls barrier of the (1 1 2) edge dislocation. The relaxation of the atomic motion was achieved using a combination of the conjugate gradient algorithm followed by quasi-dynamic relaxation (quenching). The strain increment was 106, which provides sufficient accuracy for the estimation of the stress [8]. The size p of the crystallite p pused for static simulations was 140 3/2  60 2  24 6 (lattice units), i.e. 34.7  24.3  16.8 nm3. This corresponds to the dislocation density of 1.7  1015 m2. The total number of freely mobile atoms was 1.2 million. MD simulations were performed to obtain stress–strain curves corresponding to the shear deformation under constant strain rate applied in the two opposite ½1 1  1 directions (see Fig. 1) at a finite temperature varied from 10 up to 200 K. The size of the crystallite used in MD simulations was the same as in the above-described MS simulations along y and z axes and 12.4 nm along x direction, which results in the density of mobile dislocations equal to 4.8  1015 m2. The amount of freely mobile atoms was 450,000. The process of dislocation motion was then studied in detail using visualisation tools. The dislocation line was identified using atomic disregistry analysis, i.e. by analysing the interatomic distance in the ½1 1  1 direction of atoms laying in the (1 1 0) plane to define the core atoms. The maximum relative displacement along the ½1 1  1 direction corresponds to the position of the dislocation line. In addition, high potential energy atoms and atoms with a low coordination number were recorded as well. The interatomic potential for bcc Fe developed by Ackland et al. [11] was used. The integration of Newton’s equations was performed using a constant time step equal to 5 and 1 fs for simulations at low (1–50 K) and high (above 50 K) temperatures, respectively. All calculations were done in the microcanonical NVE ensemble, where particle number, system volume and total energy are conserved if the work of external forces is taken into account. No additional temperature control was applied; the maximum increase of the temperature over the whole period of the simulation was approximately 1.5 K, observed at 10 K for an accumulated strain of 2%. Temperature, T, and applied strain rate, e, were varied from 10 to 200 K and 107 to 108 s1, respectively.

to the spreading plane, i.e. the slip plane of the edge dislocation. In order to do so, dynamic simulations have been used to apply < 1 1 1 > f1 1 2g deformation to a perfect crystal at 1 and 400 K temperatures and the corresponding stress–strain response was measured, as shown in Fig. 2. The shear modulus is estimated to be 71 and 63 GPa for the two temperatures, respectively. According to the anisotropic-elasticity theory (see for example [13]) the lattice resistance to shear in the ½1 1 1 direction in the ð1 1 2Þ and ð1 1 0Þ plane is identical and corresponds to the shear modulus determined as 1/3 (C11  C12 + C44). At 0 K, the values of C11, C12, C44 are equal to 243, 145, and 116 GPa, respectively [11]. The value of the shear modulus expected from the theory is therefore 71.3 GPa, which agrees well with the results of our and other MD simulations [8]. The model crystal thus reproduces correctly the anisotropic properties of the system. However, the decrease of the shear modulus at 400 K amounts to 11% in MD simulations, while in polycrystal iron this decrease was experimentally estimated to be 7.4% [14]. Other experimental investigations [15] provided values of elastic constants at 400 K that allows us to calculate the shear modulus at this temperature to be 68.3 GPa. Both direct estimations and experimental measurements of elastic constants in polycrystals suggest that the value of the shear modulus is underestimated in MD simulations. This discrepancy most likely should be attributed to the inaccurate reproduction of the temperature dependence of the lattice parameter by the employed here potential [11]. 3.2. Dislocation-core structure and the Peierls model Inspired from the Peierls–Nabarro (PN) model [16,17], modified in [18] and later on generalized in [19], a convenient way to represent the core structure of the edge dislocation is to build the distribution of displacements in the atomic planes on both sides of the slip plane. However, in the case of the (1 1 2) dislocation, the BCC lattice implies the addition of three (1 1 1) half-planes on the one side of the slip plane in order to create the edge dislocation, see Fig. 3.

150 Stress (MPa)

100 1K 400 K 50

3. Static calculations 3.1. Elastic response of the model crystal Before investigating the dislocation-core structure, it is of interest to evaluate the shear modulus corresponding

shear(%) 0 0.0

0.1

0.2

Fig. 2. Elastic response on <1 1 1 > {1 1 2} shear applied on the perfect crystallite (at e = 106 s1), estimated using MD simulations at 1 and 400 K.

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

The ð1 1  1Þ plane is not a mirror plane in the BCC lattice, which suggests that the (1 1 2) dislocation has an asymmetrical core structure. This aspect has already been addressed in literature [10] and will be recalled in the following. In the original treatment of the PN model, one assumes that the absolute value of the displacement u(x) in the upper (1 1 2) plane (uu(x)) is equal to that in the lower (1 1 2) plane (ul(x)) adjacent to the dislocation slip plane. This is why the two displacement components are not distinguished in the literature, for a review see [20]. In general, only the differential displacement across the slip plane, i.e. the value of the function uu(x)  ul(x) is reported [8,10]. This function increases inevitably from 0 to b. Following the construction of the Peierls dislocation, the disregistry U(x), i.e. the shift of atomic ð1 1  1Þ planes across the slip plane, is given by the function U(x) = ½ b  2|u(x)|. However, as will be seen in the following, our results show that the absolute values of uu(x) and ul(x) are not equal. In order to characterise completely the core structure of this dislocation, we report here the method to be used to estimate the two displacement functions from MS simulations. The method used to create the edge dislocation [8] leads to the atomic configuration depicted in Fig. 3b. The periodic boundary condition in the x direction implies expansion of the upper halfcrystal and compression of the lower-half-crystal by b/2. This condition imposes the values for displacements in the upper and lower planes: ul(0) = uu(0) = 1=4 b, and uu(Lb) = ul(Lb) = 1=4 b, where Lb denotes the simulation box dimension in the x direction. In this way, every atomic layer accommodates one half of the total displacement magnitude b. This condition is necessary to deduce the reference atomic positions, needed to apply the Peierls model, see Fig. 3a. The difference between these positions and

a [112]

b/2

[111]

b/3

b

Lb

Periodic boundary

Fig. 3. (a) The original Peierls unrelaxed edge dislocation constituted of two half-crystals with different number of atomic planes and (b) the unrelaxed configuration of the model crystal used to create the edge dislocation [8]. The dashed lines indicate the disregistry of atomic layers across the slip plane in the x direction.

1419

those in the relaxed dislocation structure obtained in MS simulations provides the two displacement functions. Fig. 4a shows the evolution of these displacement functions. Contrary to the one-atom lattice, considered in the PN model, the absolute values of the displacements at a given x are not equal. In the upper layer the variation of displacement is abrupt while it is smooth in the lower layer. This leads to an ambiguity in the determination of the dislocation-core width. Furthermore, the relative displacements, obtained by subtracting the displacement of atom i from that of atom i + 1, are asymmetrical, as shown in Fig. 4b. In addition, in the upper layer, the differential displacement profile exhibits a local minimum at the centre of the dislocation. This aspect cannot be explained in the concept of the PN model, and it indicates strong anisotropy in atomic forces in the twining and antitwinning directions. The difference in the absolute value of the displacements is in fact in contradiction with the PN model. The distributions of infinitesimal dislocations in both layers are different, meaning that the local equilibrium equation of the PN model cannot be applied for the (1 1 2) dislocation. On the other hand, the atomic displacements, which are responsible for the stress sxy in the upper and lower planes, are expected to be balanced by forces generated by the atomic misfit across the slip plane, resulting from the local disregistry or misfit U(x). The latter should then be evaluated by the equation U(x) = ½ b + uu(x)  ul(x). In the relaxed atomic structure, U(x) varies from U(0) = 0 to U(Lb) = b and the periodicity of the crystal along x direction is recovered. One of the major modifications to the PN model [21–23] consists in the use of the Generalised Staking Fault Energy (GSFE), c, measured in the direction of the Burgers vector in the slip plane to replace linear elasticity. According to this approach, the local misfit stress smis in the x–y plane is given by: @cðtÞ ; ð1Þ smis ðxÞ ¼  @t t¼UðxÞ where t is the magnitude of the staking fault in the ½1 1  1 direction of the (1 1 2) plane. The GSFE is responsible for the decrease of the dislocation-core size [23] and can be accurately computed using ab initio (first-principles) calculations. The comparison between the GSFE computed on the ð1 1 0Þ and the (1 1 2) planes is presented in Fig. 5a. The comparison shows that, for the asymmetrical GSFE, the dislocation-core structure is asymmetrical too. In fact, the ‘‘atomic arrangement” or the lattice discreteness is accounted for through the GSFE. This is why the compression in the upper layer or tension in the lower layer cannot balance the asymmetrical resistance to shear represented by the derivative of the GSFE, see Eq. (1). The last feature of the dislocation-core structure investigated in this work is the evolution of the differential misfit DU(x) when an external stress is applied. According to the initial PN model and other subsequent developments [16– 18,21], the application of an external stress provokes a dis-

1420

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

b 0.20

a 0.30 0.20

0.15

uu (x)

ul (x)

0.10

Δuu (x)

0.10

0.00 70

75

85 x(b) 90

80

-0.10

− Δul (x) 0.05

x(b)

-0.20 0.00

-0.30

70

75

80

85

90

Fig. 4. Displacement functions (in units of b) of atomic layers immediately above (upper layer) and below (lower layer) the slip plane. (a) absolute value and (b) relative displacements.

a 0.08

b 0.3

γ (eV/A2)

ΔΦ(b)

(112)

(110)

0.06 0.2

(112)

0.04 (110)

0.1

0.02

x(b) 0.00 0.0

x(b) 0.0

0.2

0.4

0.6

0.8

1.0

-6

-3

0

3

6

Fig. 5. Effect of the asymmetry of GSFE on the dislocation-core structure. (a) Energy of the staking fault in the ð110Þ and ð112Þ planes and (b) the differential misfit in the ð 110Þ and ð112Þ dislocation-core structures.

placement of all atoms constituting the dislocation core. The magnitude of displacement is independent of the atomic position and represents a rigid translation of the whole core. In other terms, the misfit profile is assumed to be independent of the applied stress. Only its centreof-mass is translated. To the best knowledge of the authors, no experimental or theoretical evidence so far have been provided to confirm or disprove this assumption. By using atomistic simulations, Bulatov et al. [20,23] have shown that the latter assumption leads to an overestimation of the Peierls stress. However, they did not provide any evidence on the modification of the disregistry profile with the applied stress. In addition, in our case, it is of interest to follow the evolution of the asymmetry of the core struc-

a 0.20

ture of the (1 1 2) edge dislocation when the loading proceeds in the Twinning (TD) and Antitwinning (AD) Directions. As will be shown later in this paper, the Peierls stress measured using MS technique amounts to 640 and 700 MPa in the TD and AD, respectively. In order to reveal the effect of the stress on the differential misfit shown in Fig. 5b, two static simulations corresponding to 600 MPa loading in the TD and AD were carried out. The results depicted in Fig. 6a show clearly that, not only does the centre-of-mass of the (1 1 2) dislocation move with the stress, but also its displacement profile and distribution. In Fig. 6b, we plot the difference between the atomic positions before and after loading, di(x), in the two directions. The atomic translation is not uniform in

b 0.10

ΔΦ (b)

δi(b) Twinning load

0.15

0.05 0.10 unloaded TD loaded AD-loaded

0.05 0.00 -5.0

0.00

x(b)

x(b) -2.5

0.0

Antitwinning load

2.5

5.0

-0.05 -6.0

-4.0

-2.0

0.0

2.0

4.0

Fig. 6. (a) Evolution of the differential misfit under 600 MPa shear stress applied in the twinning and antitwinning directions; (b) modification of atomic positions due to load in the two directions in comparison with their initial positions in the unloaded relaxed state.

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

either case. A strong asymmetry can be clearly seen in the figure. Consequently, the stress effect on the core structure of the dislocation is not limited to a simple translation of its centre-of-mass. On the other hand, the magnitudes of the displacement of the centre-of-mass of the (1 1 2) dislocation in the TD and AD are not equal, even though the stress applied is the same. Thus, the absolute value of the displacement of the centre of the dislocation and modification of the core structure depend not only on the absolute value of the stress, but also on the loading direction. For the applied stress of 600 MPa, the centre-of-mass of the dislocation moved by 0.42b in the TD, while it was displaced only by 0.17b in the AD. As a conclusion of this section, we can confirm that the different physical features of the dislocation-core structure, i.e. the displacements in the upper and lower layers, the local mechanical equilibrium condition and the GSFE cannot be connected through the PN model in the case of the (1 1 2) edge dislocation. A forthcoming paper is proposing a set of equations enabling the prediction of the two unknown displacement functions uu(x) and ul(x).

tion. The pre-logarithmic factor is a complex function of elastic moduli. In the isotropic-elasticity picture, the edge dislocation energy is given by: E¼

lb2 R ln ; 4pð1  mÞ ro

The energy stored in the cylinder with the dislocation along its axis has been calculated as a difference between the potential energy of atoms in the relaxed crystal with and without a dislocation. We thus avoid the delicate problem of determining the atomic-level stress and strain tensors. The additional energy corresponds therefore only to the elastic strain energy, E, generated by the dislocation. This energy was estimated for the ð 1 1 0Þ and (1 1 2) edge dislocations as well as for the ½<1 1 1> screw dislocation. Fig. 7a displays the distortion energy as a function of the cylinder radius R in the logarithmic scale. It can be seen that above a certain radius, Rc, and below a certain radius Rl, the curve exhibits a linear tendency with ln(R), as expected from the elastic theory. By definition, the dislocation-core radius is Rc and can be evaluated from the figure. Its value is slightly above 2b. Below this radius, the atomiclevel deformation is so high that the elastic theory cannot be applied anymore, nor can the line tension approxima-

a 30

b 15

(a)

E (eV/nm)

20

10

10

5

{110} edge {112} edge screw

Ln ( R /b )

Ln (R /b )

0 -1

1

3

ð2Þ

where l is the shear modulus, m is Poisson’s ratio and ro is the inner cut-off radius. From Fig. 7a, we can estimate ro by extrapolating the logarithmic part of the curve down to zero energy. For the (1 1 2) dislocation its value is 0.4b while it is close to 0.5b for the (1 1 0) dislocation. Both edge dislocations exhibit the same slope of the logarithmic part, see Fig. 7b. The energy of the (1 1 2) edge dislocation is then higher than that of the (1 1 0) dislocation. The value of ro for the ð1 1 0Þ dislocation is the same as that obtained in [24] by analysing data from [8], where a different interatomic potential (from [25]) was used. Here, we note that this value of ro for the (1 1 2) and ð1 1 0Þ edge dislocations holds also for the screw dislocation. This is unexpected, since the core structure of the ½<1 1 1> screw dislocation is known to be compact in iron [12] unlike that of the edge dislocation shown in Fig. 5b. It is of interest to evaluate the elastic moduli entering Eq. (2), because this would allow a realistic elastically-isotropic picture of dislocation to be used and thus to be benefited from advances of dislocation theory developed for isotropic materials [26,27]. Since two moduli (the shear modulus and the Poisson’s ratio) intervene in Eq. (2), we need some additional information. This is why we have computed the elastic strain energy stored around the screw dislocation, shown in Fig. 7b, using the same method as for the edge dislocation. The linear part of the strain energy enables us to deduce directly the shear modulus, which is found to be close to 62 GPa. As expected, the latter is not equal to the value of l obtained by the ½1 1  1 ð1 1 2Þ shear shown in Fig. 2. This difference is due to the fact that the estimation of the dislocation strain energy requires the integration of the energy density around the dislocation line, involving several elastic moduli, while a [1 1 1] (1 1 2) shearing involves only the C44 modulus. Considering that the value of the shear modulus for the edge and screw dislocations are the same, the ratio between the elastic energy of the screw and edge dislocations is equal to (1  m). The

3.3. Dislocation energy

E (eV/nm)

1421

5

0 -1

0

1

2

3

Fig. 7. Elastic energy stored in a cylinder around the dislocation line. (a) the case of (1 1 2) dislocation in the periodical crystal and (b) comparison between the strain energies of the two edge dislocations and the screw dislocation.

1422

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

corresponding value of Poisson’s ratio was hereby found to be 0.49. The latter is quite large and close to the theoretical limit of stability of the edge dislocation (1  2m = 0) [28], established in isotropic media. However, the values of l and m are in very good agreement with those obtained using the approach proposed by Bacon et al. [27] giving an effective isotropic picture of dislocations in anisotropic medium. The effective isotropic picture is thus consistent with the atomic-level predictions for the dislocation strain energies. The distortion energy stored in a cylinder of radius larger than 30b, is larger than the prediction of Eq. (2) as can be seen from Fig. 7a. This feature was first underlined by Osetsky et al. [8] and it is the result of the boundary conditions applied to create the edge dislocation. Eq. (2) represents the energy of the dislocation within an infinite elastic medium. The dislocation model applied here, employing periodic boundary conditions in the x direction, is expected to reproduce an edge dislocation repeated infinitely in the image boxes (periodic array of dislocations). The interaction energy (of the dislocation with its image) is not zero especially at large distances from the dislocation core. In particular, in the periodic region far from the dislocation core, the atoms experience the strain field of the dislocation in the adjacent image box. On the other hand, the fixed atomic layers in the y direction do not accommodate the strain field of the edge dislocation. This generates repulsive forces between the fixed atomic layers and the dislocation, leading to an additional increase of the crystal energy. These two features are therefore responsible for the deviation of the strain energy of the dislocation shown in Fig. 7a. In the constructed crystal the total dislocation energy is 30 eV per nanometer, which is equivalent to the energy of a single dislocation introduced in a crystal that is 10 times larger. However, a recent investigation [24] has shown that during the bowing of the pinned dislocation, the associated additional energy does not depend on the simulation box size but simply equals to the additional core energy 7–8 eV/nm. We therefore believe that the size of the MD crystal used here is large enough to simulate the motion of the dislocation in both static and dynamics conditions. 4. Motion of the edge (1 1 2) dislocation 4.1. Dislocation motion at 0 K The MS technique was used to simulate the motion of the (1 1 2) edge dislocation at 0 K, from which the corresponding stress–strain curves were obtained (see Fig. 8). After a large applied strain, the dislocation experienced a rigid motion in the TD over the Peierls valley when the resolved stress reached 640 MPa. When load was applied in the AD, the dislocation started to move only at the resolved stress of 700 MPa. At 0 K, the lattice resistance is therefore larger in the AD than in the TD. This is the direct consequence of the asymmetry in the core structure

800

Antitwinning

Stress (MPa) 600

Twinning

400 200 shear(%) 0 0.0

0.4

0.8

1.2

1.6

Fig. 8. Stress–strain loading curves at T = 0 K. The loading is performed in the twining and antitwinning directions.

of the (1 1 2) dislocation. In the absence of thermal activation, the rigid motion in both loading directions requires a large quantity of mechanical work. The two estimated Peierls stresses are much lower than that for the ½ <1 1 1> screw dislocation [29] and, at the same time, much higher than that of the (1 1 0) edge dislocation [8]. Although the Peierls stress in the AD is about 700 MPa, the displacement of the centre-of-mass of the dislocation calculated at 600 MPa, see Fig. 6, amounts only to 0.17b. Thus, even at level of stress representing 87% of the Peierls stress (in the AD), the centre-of-mass of the dislocation experiences only a small shift in its potential valley. This is direct evidence from MS simulations of the asymmetry of the motion of the edge dislocation in the {1 1 2} slip system, in line with experimental observations [5,6]. The difference in the Peierls stress only amounts to 9% and decreases further at finite temperature, as will be shown in the following. 4.2. Effect of temperature on the stress–strain curves MD simulations of plastic shear in the TD and AD were performed at different temperatures at a fixed strain rate, s1, usingp the same 107p p crystal as in MS simulations (50 3/2  60 6  24 2  a03) and adjusting the equilibrium lattice unit and integration time step depending on the considered temperature. The applied simulation conditions (i.e. strain rate and dislocation density) are equivalent to impose a fixed average velocity of the dislocation close to 8.4 m/s. Although the fixed value of the strain rate is fairly high, which is unavoidable feature of MD simulations, the resulting average velocity is reasonable due to relatively large dislocation density (4.8  1015 m2). The latter represents the important physical parameter of the simulations. Stress–strain curves extracted from the simulations are presented in Fig. 9 for the two loading directions. The critical stress needed for the dislocation motion decreases drastically with temperature. In fact, it reaches one third of the Peierls stress already at 20 K. At 200 K, the maximum stress measured over the whole simulation period does not exceed 15 MPa, independently of the direction of load, as shown in Fig. 9. At low T, e.g. 20 K, there is a remarkable difference in the stress–strain curves obtained for TD and AD loads. In particular, applying

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426 400

400 stress (MPa)

stress (MPa)

(a) twinning

1423

(b) Antitwinning

300

300

20 K 20 K

200

200

50 K 100 0 0.0

77 K 100 K 200 K 0.5

1.0

1.5

shear (%) 2.0 2.5

50 K 77 K 100 K 200 K

100 0 0.0

0.5

1.0

1.5

2.0 (%) 2.5 shear

Fig. 9. Loading curves of the (1 1 2) dislocation in (a) the twinning and (b) antitwinning directions as a function of temperature.

AD deformation, the critical stress is clearly higher and the frequency of drops on the stress–strain curve is lower, than in the case of TD load. The larger drops, seen in the stress– strain curve for AD deformation, suggest that the number of successive jumps of the dislocation between the neighbouring Peierls valleys is larger than that of the dislocation moving in the TD. The reason for this will be revealed in the following section. 4.3. Mechanism of motion Using visualisation tools we have followed the motion of the dislocation core in detail. The procedure used to identify the dislocation-core atoms was described in Section 2. Following the motion of the dislocation core, it has been revealed that the dislocation may move via different mechanisms, depending on temperature. At 0 K, the (1 1 2) edge dislocation experiences a rigid motion by displacement of the entire dislocation line over the Peierls valley. At finite temperature, the dislocation is found to move via nucleation of kink-pairs and expansion of kinks along the dislocation line, in exactly the same manner as the ½<1 1 1> screw dislocation moves at low temperature in BCC Fe [12]. In Fig. 10, we show the loading curve for the AD deformation at 20 K on the left-hand side and the dislocation line shape during the jumps, starting at almost 0.5% of elastic shear strain, on the right-hand side. Every stress-drop corresponds to several kink-pair nucleation events. Although the stress decreases during the drop, successive events of kink-pair nucleation, following one another, occur after the recombination of the first kinkpair. The origin of this chain-like motion will be discussed and clarified later on. The configuration of atoms associated with the dislocation core at equally spaced time increments (2.5 ps) and starting from the shape of the dislocation at c = 0.5% are depicted in Fig. 10b. The stress-drop marked by a dashed rectangular in Fig. 10a corresponds to the displacement of the dislocation for a distance equal to four Peierls valleys. The passage of the dislocation from one to another valley proceeds via nucleation of kink-pairs, within a relatively small spacing, followed by the expansion of the kinks in opposite directions along the dislocation line.

The visualisation of the positions of the kinks at different moments (see example in Fig. 10b) makes it possible to evaluate the velocity of the kinks as a function of the applied stress. The latter is assumed to decrease linearly with time (during the drop), so that regarding Fig. 10a, the stress decreases from 260 down to 160 MPa at a constant rate of 5.5 MPa per 2.5 ps. The computed velocities are shown in Fig. 11 and are found to increase with stress, as expected. The stress dependence of the velocity of the left and right kinks seems to be the same within the calculation error. Assuming a linear dependence of the velocity on the stress, i.e. v = b(s  so)/B, where so is a threshold stress, one may deduce the viscous coefficient B. The best linear fit applied to the data presented in Fig. 11 provides the following values: 94 MPa for so and 3  105 Pa s for B. The latter is the same as was deduced from MD simulations used to model the (1 1 0) edge dislocation [8]. For more details see Fig. 15 in [8] where the viscous coefficient can be easily estimated using the linear part of the curve. The obtained coincidence can be explained taking into account that kinks are of mixed character [30,31] rather than pure screw type. On the other hand, at high enough temperature, the nucleation of kink-pairs is no longer the controlling mechanism [32]. Then, the dislocation mobility is controlled by the stress necessary to induce a collective motion of kinks. Therefore, at high temperature, the velocity of the whole dislocation, as well as the velocity of separated kinks, should be of the same order. Since the motion of the (1 1 2) edge dislocation at low temperature is controlled by the nucleation of the kinkpairs, one expects the dislocation jump rate to follow the same law as that for the screw dislocation, so that the dislocation velocity is proportional to the dislocation length and to the Boltzmann’s factor. One can thus apply the same treatment as proposed by Domain et al. [12], in which the following equation was derived to deduce the value of the critical stress from the stress–strain curve obtained from MD simulations:   kT Vs ln exp ; ð3Þ sc ¼ V kT where V is the activation volume and, k Boltzmann’s constant. Here, we note that the value of sc does not depend

1424

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

b

a 400

50 stress (MPa) 40

300

30 200 20 5

100 4

3

10 Position (ao)

6

0 0.0

0.5

1.0

0 1.5

Fig. 10. (a) Stress strain curve and dislocation position at T = 20 K and (b) configuration of atoms ascribes to the dislocation line at 2.5 ps time intervals.

2000 Kink velocity (m/s)

1500 1000 500 τ (MPa)

0 150

200

250

300

Fig. 11. Velocity of kinks as a function of the applied shear stress at T = 20 K. Circles and squares refer to the left and right kinks.

significantly on the activation volume. An approximate value of V is therefore enough to obtain a good estimate of the critical stress. In order to do so, we performed MD simulations at an intermediate temperature, T = 77 K, using different strain rates. Fig. 12 shows the evolution of the maximum stress as a function of the strain rate. Using the results of these simulations, the activation volume can be calculated using the following equation: V ¼

kT lnðe2 =e1 Þ s 2  s1

ð4Þ

From Fig. 12, we can see that the evolution of the stress is a logarithmic-like function, which suggests that the activation volume is not strongly altered in the range of stres-

120 Maximum stress (MPa)

100

80 6 -1

strain rate(10 s ) 60 0

20

40

60

80

100

Fig. 12. Evolution of the maximum stress as a function of the strain rate at T = 77 K.

ses corresponding to the motion of the (1 1 2) edge dislocation in the chosen simulation conditions. Eq. (4) provides a value of the activation volume close to 4b3 when evaluated using e1 = 107 and e2 = 108 s1. This value is small compared to the one measured experimentally at low temperature in iron [6,33]], where the flow stress is known to be controlled by the nucleation of the kink-pairs on screw dislocations. The difference between the two activation volumes arises from the difference in lattice friction for the screw and edge dislocation. At the same temperature and stress, the area swept by the dislocation with the higher Peierls stress is correspondingly higher. Assuming that the value V = 4b3 does not vary strongly with temperature, Eq. (3) can be used to evaluate the critical stress as a function of temperature. The results are shown in Fig. 13. The value of the critical stress at 0 K is equal to the Peierls stress computed from MS simulations presented above. From Fig. 13, one sees that the critical stress decreases strongly with temperature. Lattice friction is higher in the AD at low temperature. Above 50 K, the lattice friction is identical in the TD and AD. At 200 K, the edge dislocation can achieve measurable motion at very low stresses, independently of the direction of load. 4.4. Discussion The investigation of atomic positions and displacements in a crystal containing the (1 1 2) edge dislocation shows that the core structure is fundamentally different from that of other edge dislocations, for instance the (1 1 0) edge dislocation and the 12 < 1 1 0 > f1 1 1g edge dislocation of the FCC lattice. The core of the (1 1 2) dislocation was found to have asymmetrical structure, see Fig. 4, in the upper layer as well as in the lower layer. At this, the amplitude of atomic displacements in the upper layer is not equal to that in the lower layer. In addition to the translation of the centre-of-mass of the dislocation, the application of an external stress (below the Peierls stress) leads to a strong modification of the core structure. The Peierls–Nabarro model cannot provide an accurate description of the ½<1 1 1> edge dislocation in the (1 1 2) slip plane. The

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426 800

τ c (MPa) 600 400

Antitwinning direction Twinning direction

200 Τ (Κ)

0 0

50

100

150

200

Fig. 13. Evolution of the critical stress as a function of the temperature for both loading directions: twinning and antitwinning. The critical stress is obtained at a fixed velocity.

asymmetry is the natural consequence of the asymmetry of the generalised staking fault energy in this plane along the <1 1 1> direction. It is also expected to be the origin of the direction-dependent friction stress (see Fig. 8), widely accepted in the literature. All features shown in this work indicate that the glide of the (1 1 0) edge dislocation is easier than that of the (1 1 2) edge dislocations. On the one hand, the Peierls stress in the (1 1 2) plane is higher in the AD than in the TD. At finite temperature, the motion of the (1 1 2) edge dislocation in both directions occurs through the nucleation and propagation of kink-pairs, as in the case of screw dislocations [12]. However, the velocity of kinks on the screw dislocation was found to be close to the sound velocity, whereas it was lower in the case of the (1 1 2) edge dislocation (see Fig. 11), which is an important observation. The reason for this could be the level of the stress applied, which is much higher in the case of the screw dislocation. Unlike the case of the screw dislocation [12], the sequential nucleation of kink-pairs found in the present study is correlated to the annihilation of the two kinks, generated by the previous kink-pair nucleation event. In other words, the first kink-pair requires a large level of stress to nucleate (see Fig. 10). When one of the propagating kinks crosses the periodic boundary, imposed in MD simulations, it moves towards the other kink. The subsequent annihilation results in the energy release (heat), close to the double formation energy of a separated kink. If one neglects the interaction energy and formation entropy, this energy corresponds to the energy of a dislocation segment of a length b, which can be easily deduced from Fig. 7 by taking an average of the energy of the screw and edge dislocation. It was calculated to be about 0.5 eV per kink, meaning that the energy released due to the annihilation is close to 1 eV. One half of this energy dissipates into heat resulting in a local increase of the atomic kinetic energy, and one half goes to the potential energy of the crystal. If the corresponding kinetic energy is distributed over a small volume of about 100 atoms, DT is close to 50 K. The influence of this local heating on the events of kink-pair nucleation depends on two factors: the ratio of DT/T; and the activation energy of the process, entering Boltzmann’s factor. At

1425

high temperature, the ratio DT/T is low so that the influence of the annihilation event of the two kinks on the nucleation probability is small and the nucleation of the two successive kink-pairs can be considered as uncorrelated processes. At low temperature, say 20 K, the nucleation probability is tremendously increased. Hence, the probability for the nucleation of the subsequent kink-pair increases in the region where the two pre-existing kinks annihilated. This explains why at low temperature a sequence of kinkpair nucleation-correlated events was observed, even though the stress continued to decrease during the drop (see Fig. 10). This chain process, however, stops when the local increase in temperature no longer compensates the rise of the activation energy provoked by the decrease of the applied stress. When the kink-pair nucleation stops, the continuous shearing of the model crystal makes the stress increase again, until the next kink-pair (this time not correlated to the previous kinks) will be nucleated. According to the above discussion, only the loading part coming before the nucleation of the first kink-pair (uncorrelated) is representative of the lattice friction of the material. As a simple consequence of this finding, in experimental conditions, where the stress is a smooth function, the motion of edge dislocations at low temperature is expected to occur by succession of avalanches. When the dislocation is pinned at its extremities, the edge dislocation finishes its avalanches by acquiring a significant curvature, resulting from accumulation of signed kinks on each extremity. From Fig. 7, we can see that the total energy of the (1 1 2) dislocation (i.e. the sum of the core and elastic energies) is larger than the one of the (1 1 0) dislocation. Although the difference amounts only to 5%, associated with the core contribution, the bowing of the edge dislocation in the (1 1 2) plane requires larger stress than in the (1 1 0) plane. Since the dislocation bowing process was found to be controlled by the core energy [24], as a first approximation, the stress necessary to bow out the (1 1 2) edge dislocation should be larger by 5% than that necessary to bow the (1 1 0) dislocation. Thus, one may expect that the critical stress to overcome obstacles pinning the dislocations moving in the (1 1 2) plane would be higher than that for the (1 1 0) plane, since the unpinning from obstacles always requires some curvature of dislocations. The temperature dependence of the friction stress, presented in Figs. 9 and 13, in the TD and AD, indicates that the activation of the 12 f1 1 2g < 1 1 1 > slip system is unlikely to occur at very low temperature. The lattice resistance was, however, found to decrease rapidly with temperature. The decrease should be even larger in experimental conditions, where the strain rates are lower than those in MD simulations by a few orders of magnitude [12]. Therefore, at high enough temperature, say larger than 100 K, the activity of edge dislocations in the {1 1 2} and the {1 1 0} planes is expected to be comparable. This conclusion holds for both loading conditions, i.e. in the twinning and antitwinning directions. The characterisation of the slip activity

1426

G. Monnet, D. Terentyev / Acta Materialia 57 (2009) 1416–1426

in the 12 < 1 1  1 > f1 1 2g systems is thus suspended until the complete description of the mobility of the screw dislocation in the {1 1 2} plane is assessed. Prior to the conclusions, we would like to give a comment about the sensitivity of the obtained results to the interatomic potential. Over the last ten years the most common cohesive model used to simulate bcc Fe has been the potential developed by Ackland et al. [26]. The newer version employed here predicts principally different core structure of the ½<1 1 1> screw dislocation, while giving a better description of the properties of self interstitials in comparison with the potential from Ref. [26]. In addition, it has been recently shown that the Peierls stress of the ½<1 1 1 > {1 1 0} edge dislocation is about 80 MPa with the potential employed here [34], that is about three times higher than for the model from Ref. [26]. No essential difference in the core structure of the ½<1 1 1 > {1 1 0} edge dislocation was revealed when comparing these two potentials [34] and the difference in the Peierls stress was unexplained. We have carried out a set of static calculations using the older interatomic potential to find out a possible influence on the already presented results. It turned out that the main statements, i.e. the asymmetry of the core structure and non-equivalence of the twinning and antitwinning loads, remain unchanged. The obtained Peierls stress, however, was found to be systematically lower (by about 100 MPa) with the potential from Ref. [26]. The relative difference of about 60 MPa between the Peierls stresses for the two opposite loading directions also remained unchanged. We therefore see that the two potentials give slightly different estimations of the Peierls stress for the (1 1 2) edge dislocation. One can, however, ensure that the difference becomes insignificant in simulations performed at a finite temperature, as has been found in the case of the (1 1 0) edge dislocation [34]. 5. Conclusions The following conclusions are drawn based on the results obtained in this work. 1. The ½<1 1 1 > (1 1 2) edge dislocation exhibits asymmetrical core structure. External applied stress is found to modify the displacement profile within the dislocation core, therefore, the standard Peierls–Nabarro model cannot be applied directly. 2. The origin of the asymmetry in the core structure is directly related to the asymmetry of the general stacking fault energy of the {1 1 2} plane in the <1 1 1> direction. As a consequence, the Peierls stress for the (1 1 2) edge dislocation depends on the particular direction of load. The Peierls stress in the twining direction is lower than in antitwinning direction.

3. The motion of the (1 1 2) edge dislocation was observed to occur via nucleation and propagation of kink-pairs, similar to the mechanism of motion of the screw dislocation, if the temperature does not exceed 100 K. The velocity of kinks on the line of the (1 1 2) edge dislocation was calculated and found roughly proportional to the applied stress. At higher temperature, the motion is controlled by the drift of thermal kinks, found in large concentration along the dislocation line. 4. The Peierls stress for the (1 1 2) edge dislocation exceeds significantly that of the (1 1 0) edge dislocation, whereas it is about twice as low as the Peierls stress for the screw dislocation. By increasing temperature, the critical stress needed to move the (1 1 2) edge dislocation decreases drastically and above 100 K it is of the same order as that of the (1 1 0) edge dislocation. The difference in the critical stresses in the AD and TD also becomes negligible.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

Christian JW. Metall Trans A 1983;14:1237. Taylor G. Prog Mater Sci 1992;36:29. Jaoul B, Gonzalez DJ. Mech Phys Solids 1961;9:16. Vogel FL, Brick RM. Trans AIME 1953;197:700. Spitzig WA, Keh AS. Acta Metall 1970;18:611. Kuramoto E, Aono A, Kitajima K. Scripta Metall 1979;13:1039. Yamaguchi M, Vitek V. J Phys F Metall Phys 1975;5:11. Osetsky Yu N, Bacon DJ. Model Simul Mater Sci Eng 2003;11:427. Yamaguchi M, Vitek V. J Phys F Metall Phys 1975;5:1. Masuda K, Kobayashi K, Sato A, Mori T. Phil Mag 1981;43:19. Ackland GJ, Mendelev MI, Srolovitz DJ, Han S, Barashev AV. J Phys Condens Matter 2004;16:1. Domain C, Monnet G. Phys Rev Lett 2005;95:215506. Martin HS. Elasticity. Burlington (USA): Elsevier Academic Press; 2005. Ghosh G, Olsson GB. Acta Mater 2002;50:2655. Dever DJ. J Appl Phys 1972;43:3293. Peierls RE. Proc Phys Soc 1940;52:23. Nabarro FRN. Proc Phys Soc 1947;59:256. Wang JN. Acta Mater 1996;44:1541. Ngan AHW. J Mech Phys Solids 1997;45:903. Bulatov VV, Cai W. Computer simulations of dislocations. Oxford: University press; 2006. Christian JW, Vitek V. Rep Prog Phys 1970;33:307. Joos B, Duesbery MS. Phys Rev Lett 1997;78:566. Bulatov VV, Kaxiras E. Phys Rev Lett 1997;78:4221. Monnet G. Acta Mater 2007;55:5081. Ackland GJ, Bacon DJ, Calder AF, Harry T. Phil Mag A 1997;75:713. Bacon DJ, Barnett DM, Scattergood RO. Prog Mater Sci 1978;23:51. Bacon DJ. In: Bilby BA, Miller KJ, Willis JR, editors. Fundamentals of deformation and fracture. C.U.P; 1985. p. 401. De Wit G, Koehler JS. Phys Rev 1959;116:1113. Chaussidon J, Fivel M, Rodney D. Acta Mater 2006;54:3407. Duesbery MS, Xu W. Scripta Mater 1998;39:283. Dorn JE, Rajnak S. Trans Metall Soc AIME 1964;230:1052. Hirth JP, Lothe J. Theory of dislocations. 2nd ed. New York: Wiley; 1982. p. 376. Spitzig WA. Mater Sci Eng 1973;12:191. Terentyev D, Bacon DJ, Osetsky Yu N. J Phys Condens Matter 2008;20:445007.