Computer simulation of carbon diffusion and vacancy–carbon interaction in α-iron

Acta Materialia 55 (2007) 1–11 www.actamat-journals.com

Computer simulation of carbon diffusion and vacancy–carbon interaction in a-iron K. Tapasa a, A.V. Barashev a

a,*

, D.J. Bacon a, Yu.N. Osetsky

b

Materials Science and Engineering, Department of Engineering, University of Liverpool, Brownlow Hill, Liverpool L69 3GH, UK b Computer Sciences and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6138, USA Received 21 March 2006; accepted 14 May 2006 Available online 13 November 2006

Abstract Static and dynamic properties of the interstitial carbon atom and vacancy–carbon atom complexes in a-iron are modelled by a molecular dynamics (MD) method using a pair interatomic potential for the iron–carbon interaction and two different many-body potentials for the iron matrix. The diffusion parameters of an interstitial solute in iron are obtained by MD simulation for the first time. The binding energy and migration energy of a vacancy–carbon complex are also obtained: the complex is immobile and has higher energy for dissociation than the carbon atom migration energy. The results are compared with recent ab initio calculations and experimental data from the literature. Experimental data on the recovery stages of electron-irradiated Fe–C are analysed using rate theory equations and found to be consistent with the ab initio calculations for diffusion of a vacancy and its interaction with carbon atoms.  2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Iron; Carbon interstitial; Carbon diffusion; Vacancy–carbon complex; Molecular dynamics

1. Introduction The properties of iron–carbon (Fe–C) solid solutions have been investigated extensively due to their particular relevance in steel technology. It is known from experiment that the interaction of carbon (C) interstitial atoms with lattice defects significantly affects the mechanical properties of ferritic steels, e.g. their strength [1], and the properties of point defects. Although it is widely accepted that vacancy (V) diffusivity is significantly decreased due to trapping by C atoms, there has been considerable uncertainty concerning the magnitude and form of the interaction. Vehanen et al. [2] deduced a value of 0.85 eV for the binding energy of a V–C complex, which, when added to a vacancy migration energy value of 0.55 eV in pure Fe, is consistent with an effective vacancy migration energy of 1.35 eV for steels. However, V–C binding energy values of 0.41 and 1.1 eV *

Corresponding author. Tel.: +44 151 794 5384; fax: +44 151 794 4675. E-mail address: [email protected] (A.V. Barashev).

have also been reported [3,4], the former estimate being close to the result of 0.47 eV obtained in recent ab initio calculations [5]. Vehanen et al. [2] argued that the value of 0.41 eV found by Arndt and Damask [3] must be associated with the binding of C atoms to V–C complexes and vacancy agglomerates. In addition, the experiments on V–C interaction are performed using electron irradiation to build up a measurable concentration of non-equilibrium vacancies, when self-interstitial atoms are also produced. The role of the self-interstitial atoms (I) in trapping C atoms is unclear. Little and Harries [6] considered the binding energy of 0.41 eV found in Ref. [3] to be more probably associated with a C atom–interstitial cluster complex, although the later experiments [2,4] indicated the C–I binding energy to be much smaller, i.e. 0.1 eV, and ab initio calculations actually give a negative energy in the range 0.2 to 0.4 eV for this pair [5]. Hence, it is desirable to use computer simulation methods to investigate the migration properties and defect interactions of C solute, but an Fe–C interatomic potential set is required and the development of such potentials has lagged behind those for pure

1359-6454/$30.00  2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2006.05.029

2

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

a-Fe and Fe alloys with substitutional solute, and only a few Fe–C potentials have been proposed to date. Johnson et al. [7] developed a pairwise potential set for the a-Fe–C system by fitting to the experimental value of 0.86 eV for the migration energy of C atoms in a-Fe [8], zero activation volume of migration [9] and a V–C binding energy of 0.41 eV. No C–C interaction was assumed in this model. This potential reproduces satisfactorily the experimental value of the energy of solution of a C atom in a-Fe relative to that in Fe3C [10] and gives the octahedral site as the most stable one for a C interstitial, consistent with experiment [11] and ab initio calculation [5]. There have been several attempts to improve the description of C atoms in Fe using empirical potentials. Rosato [12] rescaled the Fe–C potential of Johnson et al. [7] to use with a Finnis–Sinclair (F–S)-type many-body potential for Fe [13]. Again, the C–C interaction was ignored and, although this model gives the octahedral site as the most stable for C and provides a better description of the elastic properties of a-Fe, it yields 1.14 eV for the C migration energy, which is much too high. Ruda et al. [14] used embedded atom method (EAM)-type potentials fitted to ab initio data for metastable carbide FeC with B1 structure. The C–C interaction was described using the Tersoff [15] potential with no angular dependency. The equilibrium lattice constant, bulk modulus and cohesive energy of stable and metastable carbides were reproduced, but, in contradiction with experiment, the tetrahedral site is the most stable position for a C atom in a-Fe with this potential. Lee [16] has developed very recently a modified EAMtype potential set for Fe–C by fitting to experimental information and ab initio data for the dilute heat of solution of C, the location of a C interstitial, the migration energy of C in body-centred cubic (bcc) Fe, and to a first-principles calculation result for the cohesive energy of a hypothetical NaCl-type Fe–C. It produces a C–V binding energy of 0.9 eV, which is close to that predicted by Vehanen et al. [2] (i.e. 0.85 eV), but much higher than the value of 0.41 eV deduced by Arndt and Damask [3] and the ab initio value of 0.47 eV. Domain et al. [5] investigated properties of C atoms in a-Fe and their interaction with intrinsic point defects and small clusters by ab initio calculations based on the density-functional theory. It was found, in particular, that a C atom occupies an octahedral site and the activation energy of C atom jumps, i.e. the energy difference between octahedral and tetrahedral configurations, is ECm  0:90 eV (for a supercell containing 128 bcc lattice sites), which is 0.05–0.1 eV higher than the experimental values. Using a similar approach, but higher cut-off energy and density of k-points, Jiang and Carter [17] calculated this energy to be 0.86 eV. In the present work we investigate thermally activated diffusion of C atoms in bcc Fe and V–C interaction by the molecular dynamics (MD) method. Only molecular statics (MS) modelling (i.e. representing 0 K) of the interaction

between a C atom and some point defects has been performed with the potentials mentioned above. Since we are not aware of any potential set that accurately represents the Fe–C system, the approach adopted here is to study a ‘model’ system for an octahedral interstitial solute in a bcc metal matrix that has many similarities to the real alloy. We have used for this purpose the pair potential of Johnson et al. [7] for the Fe–C interaction together with either the many-body F–S-type potential for the Fe–Fe interaction of Ackland et al. [18] or the more recent EAM potential of Ackland et al. [19]. We also use the data obtained, together with results from the ab initio calculations referred to above, to review and interpret the earlier experimental information on C-point defect interactions. The paper is organised as follows. In Section 2, the calculation scheme and the potential sets are described. In Section 3, the MD results on C diffusion, comparison with ab initio and experimental data and a statistical analysis of the C atom jumps are presented. In Section 4, the results of MS and MD calculations of the interaction properties of a C atom with a vacancy are presented and compared with the ab initio data [5]. In Section 5, the results of annealing experiments are analysed on the basis of rate theory equations and the ab initio data. Conclusions are drawn in Section 6. 2. Calculation model Two atomic-scale techniques are used here to study static and dynamic properties of a C atom and V–C complexes in an Fe matrix. First, MD is used to model thermally activated motion of an interstitial C atom, and, second, MS relaxation using a combination of conjugate gradients potential energy minimisation and quasi-dynamic quenching is used to study different configurations at zero temperature and to calculate formation, binding and migration energies. In both cases the calculations were performed for a simulation box containing 2000 bcc lattice sites with constant volume and periodic boundary conditions. The F–S potential of Ackland et al. [18] and recently proposed EAM potential of Ackland et al. [19] were used for the Fe–Fe interaction, while the Fe–C interaction was described by the pair potential of Johnson et al. [7]. Both potential sets allow a correct description of the elastic properties of a-Fe, while retaining almost the same properties of a C atom as the original model of Johnson et al. [7]. This was verified by static simulation, as follows. The octahedral site was found to be the most stable position for a C atom. Both Fe potentials give the insertion energy of a C atom into a-Fe (difference in the crystal energy with and without C atom) of 1.23 eV, which is close to the corresponding energy of 1.34 eV calculated using the pair potential of Johnson et al. [7] for Fe–Fe interactions. The tetrahedral site is the saddle point for a C jump from one octahedral site to another, as with the potential set of Johnson et al. [7] and in the ab initio calculations by Domain et al. [5].

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

3

Table 1 Relaxation of Fe atoms around a C atom in a-Fe at 0 K Data from

Most stable configuration

Dd 1 =d 01 (%)

Dd 2 =d 02 (%)

This work, F–S This work, EAM Johnson et al. [7] Lee [16] Domain et al. [5]

Octahedral Octahedral Octahedral Octahedral Octahedral

+20.92 +23.00 +22.56 +28.70 +24.30

4.34 4.18 4.36 3.80 1.80

d 0i and Ddi are the distance between the C and an Fe atom at the ith nearest-neighbour distance before relaxation and the change of this distance after relaxation, respectively.

The energy barrier for C atom migration, ECm , was estimated by displacing the C atom and allowing it to relax in the plane perpendicular to the jump direction, whilst the Fe atoms relaxed fully. ECm was found as the difference in energy between octahedral and tetrahedral positions of a C atom at 0 K and is equal to 0.83 and 0.82 eV for the F–S [18] and EAM [19] Fe–Fe potentials, respectively, and thus close to 0.86 eV fitted by Johnson et al. [7] and the experimental values [8,20,21] that fall in the range 0.81–0.86 eV. The volume expansion due to a C atom is 0.3X, where X is the Fe atomic volume, and is smaller that the estimated 0.8X deduced from the experimental value of lattice parameter change with C composition [22]. However, the relaxation of Fe atoms neighbouring a C atom obtained by the combination of the Fe–C potential of Johnson et al. with either of the many-body Fe–Fe potentials is in good agreement with that found in the recent ab initio calculation by Domain et al. [5]. Fe atoms in the first nearest neighbour sites are repelled in Æ1 0 0æ directions, while those in the second nearest neighbour sites are attracted inwards along Æ0 1 1æ directions. Relative displacements of Fe atoms from lattice sites for several models are shown in Table 1. This analysis shows that both sets of the mixed pair and many-body potentials reproduce some basic properties of C atoms in Fe. 3. Diffusional properties of carbon atoms In this section we present results of MD modelling of C atom diffusion in Fe with the F–S-type potential of Ackland et al. [18]. The simulations were carried out for the temperature range from 1000 to 1800 K for up to 50 ns. (High temperatures are required to generate high statistics of C atom jumps in MD simulation, and we note that the potential gives the bcc crystal structure at all temperatures.) The thermal expansion of the lattice was taken into account by a corresponding increase of the lattice parameter, a0, to maintain zero pressure. Time integration was performed using the velocity Verlet algorithm with a variable time step, ts, which was controlled by fixing the maximum displacement of the fastest atom in the system at each step to be 0.005a0 [23]. The mean value of ts was from 0.6 to 0.8 fs depending on the temperature. During simulations, the jump direction, displacement and coordinates of

Fig. 1. MSD of a C atom at various temperatures.

the solute atom were recorded every time a jump occurred from one octahedral position to another. 3.1. Diffusion coefficient of a carbon atom Fig. 1 presents the results of MD simulations of the time dependence of the displacement squared of a C atom at different temperatures. This data allow the diffusion coefficient of a C atom to be calculated via the Einstein equation: DC ¼

hR2 i ; 6t

ð1Þ

where ÆR2æ is the mean square displacement (MSD) of a C atom for time t. One trajectory of the solute atom was simulated for each temperature and, following Osetsky [23], the MSD was estimated by decomposing it into segments of equal time, t, and taking the average over N = tsim/t segments, where tsim is the total simulation time: hR2 i ¼

N 1 X R2 : N i¼1 i

ð2Þ

Fig. 2 shows DC as a function of t. It can be seen that DC decreases with increasing t at small t, indicating a tendency of a C atom to return back to its original position. With increasing t this tendency disappears and the trajectory segments become uncorrelated [23]. The best estimate of DC was chosen to correspond to the value of t when DC saturates and its standard error is minimal. For example, at 1800 K the best estimate is DC = 1.9 · 107 m2/s (see dashed line in Fig. 2). This corresponds to t  100 ps, which is much greater than the mean time delay between sequential jumps, s, which is 2 ps for this temperature (see next section). Fig. 3 shows the temperature dependence of DC. An Arrhenius dependence is observed: DC ¼ DC0 expðbECm Þ; 1

ð3Þ

where b = (kBT) , kB being the Boltzmann constant and T the absolute temperature. The pre-exponential factor DC0

4

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11 Table 2 Comparison of the activation energy for C atom migration, ECm , and the pre-exponential factor in the diffusion coefficient, DC0, obtained by different methods Method

DC0 (107 m2/s)

EC m (eV)

Ref.

MD, F–S Static, F–S Static, EAM Static, pair potential

1.89 ± 0.01 – – –

0.71 ± 0.02 0.83 0.82 0.86

This work This work This work [7]

Ab initio, static

2.13 1.44

0.92 0.86

[5] [17]

20 6.20 1.67

0.86 0.83 0.81

[8] [20] [21]

Experiment: 243 K < T < 1043 K T > 623 K 243 K < T < 347 K

Fig. 2. Dependence of the estimate of diffusion coefficient on the duration of trajectory segments in the decomposition of a C atom trajectory at various temperatures. Error bars show the standard errors.

Arrhenius plots of the C diffusivity in a-Fe are known to exhibit positive deviations from linearity at high temperature [21]. This is probably the main reason for the scatter of the experimental data for DC0 presented in Table 2, since the treatments were performed for different temperature ranges, i.e. below 347 K in Ref. [21], above 623 K in Ref. [20] and from 243 to 1043 K in Ref. [8]. The activation energy for migration is similar, however. It is interesting to note that the ab initio [5,17] values of ECm are 0.05–0.1 eV higher than the experimental value [8,20,21]. Since these values also correspond to 0 K, it is tempting to speculate that the discrepancy is partly due to dynamic effects. Assuming these effects would lower the ab initio value by 0.1 eV, as they did in our simulations, it is reasonable to suppose that an improved potential should be fitted to the static value for C atom migration energy found in ab initio calculations, rather than the experimental value. 3.2. Analysis of carbon atom jumps

Fig. 3. Temperature dependence of the C atom diffusion coefficient.

is estimated to be 1.89 · 107 m2/s and the activation energy ECm is 0.71 eV. Table 2 presents comparison of the MD results with the ab initio calculations [5,17] and experimental data [8,20,21]. As can be seen, our estimate of DC0 is close to the experimental value derived from the low-temperature measurements (below 74 C) [21] and the ab initio [5,17] value estimated using harmonic transition-state theory [24] and the random walk formulation of interstitial diffusion in a bcc lattice [25]. Also, the dynamic value of ECm is lower by 0.1 eV than the experimental and static values, which were fitted in the potential of Johnson et al. [7]. There are some uncertainties in such a comparison with experimental data, however. In addition to experimental uncertainties, perhaps the most important one is that the

A C atom has four possible jump directions from one octahedral site to another, which, with respect to the direction of the previous jump, are forward, backward and aside (to the left and to the right). Information on this was recorded during simulations and special care was taken in identification of a solute atom jump. Application of the usual Wigner–Seitz cell method [23], in which a jump is counted when the migrating atom crosses the boundary between two octahedral sites, resulted in a very high jump frequency and a high proportion of backward jumps. This occurred due to the significant difference in properties (mass, atomic size and interaction) of C and Fe atoms. A C atom moves faster on average than Fe atoms, while Fe atoms need some time to relax around this solute every time it jumps to a new octahedral position. This time is longer than the reciprocal of the C atom attempt frequency and, as a consequence, in the majority of cases the C atom jumps backward before sufficient relaxation of Fe atoms occurs. To reduce the number of unsuccessful jumps we

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

imposed conditions for the C atom to be counted as resident in a new octahedral position by defining a successful jump as that when the C atom moved to a cube of linear size 0.1a0 centred on the new octahedral position. The frequency, mC, of jumps in each direction thus defined was treated using the Arrhenius equation mC ¼ mC0 expðbECm Þ:

ð4Þ

The mean time delay between two sequential jumps, s, which is the reciprocal of the total jump frequency forward mtot þ mbackward þ 2maside , was found to be 50, 19.2, C ¼ mC C C 4.8 and 1.9 ps at 1000, 1200, 1500 and 1800 K, respectively. The results for each type of jump are presented in Fig. 4. They show that the jump frequency obeys the Arrhenius dependence at all temperatures except 1000 K, the lowest temperature studied, for which the number of backward jumps is much higher than those in other directions. As a consequence, there is a temperature dependence of the correlation factor for C atom diffusion, fc, which is defined as the ratio of the C diffusion coefficient obtained in MD, DC, to that corresponding to the random walk, DRW C : fc ¼

DC ; DRW C

ð5Þ

where, since the jump distance of a C atom from one octahedral site to another is equal to half of the lattice parameter, DRW ¼ mtot C C

a20 : 24

5

relax around a light interstitial C atom every time it jumps towards another octahedral position. Strictly speaking, the jump definition, and therefore the residence conditions, should be temperature-dependent. Since there are no peculiarities in the C diffusion coefficient obtained at 1000 K via treatment of MSD (see Fig. 3) as compared with the other temperatures, we assume that the diffusion mechanism does not change within the temperature range considered; i.e. it is a random walk on the lattice of octahedral sites, as follows from the analysis of C atom jumps at higher temperature. To estimate an effective jump frequency at 1000 K we applied the resident requirement time (RRT) criterion [23]. According to this, a C atom must spend a longer time than the RRT in the vicinity of a new site in order to be counted as resident in that position. We found that for an RRT of 2.5 ns the number of backward jumps becomes approximately equal to those in other directions, resulting in an increase of fc from 0.66 (without RRT criterion) to unity. The corresponding frequency of C atom jumps is now related to the random walk of a C atom. As shown in Fig. 4, after application of the RRT criterion, the jump frequency obeys the Arrhenius dependence over the entire temperature region studied. Note also that the activation energy of the jumps is found to be the same as that for the diffusion of a C atom, which is another confirmation that a C atom executes random walk. 4. Properties of a V–C complex

ð6Þ

The correlation factor calculated using Eq. (5) is close to unity at 1200, 1500 and 1800 K, but is equal to 0.66 at 1000 K. This indicates that our definition of jumps fails at this temperature. As already discussed above, this is probably because Fe atoms do not have sufficient time to

Fig. 4. Temperature dependence of jump frequency of C atom in directions relative to that of the previous jump. The data for 1000 K are shown for two definitions of jumps: with and without the RRT criterion.

4.1. Static simulations The binding energy of a V–C complex, ECV , was calcub lated using static relaxation for both potential sets described in Section 2 and configurations shown in Fig. 5. (Binding energy is defined as the difference in the potential energy between systems of non-interacting and interacting defects: a positive value implies attraction and vice versa.) The results are shown in Table 3, together with the ab initio data of Domain et al. [5] calculated for a box size of either 54 or 128 bcc lattice sites. As seen from the table, qualitative agreement is obtained with the ab initio calculations, for there is a strong binding (0.51 eV for both

Fig. 5. Simulated configurations of C-vacancy complex. Open circle, Fe atom; filled circle, C atom; square, vacancy.

6

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

Table 3 C atom–vacancy binding energy ECV (eV) for configurations in Fig. 5 b Configuration

MD (this work)

Ab initio (Ref. [5]) 54 atoms

ECV b 1 2 3 4 5

+0.51 0.48 0.51 0.25 +0.16

dC–V a

(+0.51) (0.45) (0.38) (0.23) (+0.17)

0.36 0.71 0.00 0.20 1.21

(0.35) (0.72) (0.00) (0.20) (1.20)

128 atoms

ECV b

dC–V

EbCV

dC–V

0.44 0.09 0.33 0.11 –

0.40 0.64 0 0.24 –

0.47 0.01 – – –

0.41 0.64 – – –

The distance between a vacancy and a C atom, dC–V, is in units of the lattice parameter. a The data obtained are for F–S-type Fe–Fe potential of Ackland et al. [18]; the data in parentheses are for EAM-type Fe–Fe potential of Ackland et al. [19].

potential sets) when the C atom is in the nearest-neighbour octahedral site to a vacancy (configuration 1 in Fig. 5), whereas other configurations show repulsive interaction (except for configuration 5, not studied ab initio). Our result for configuration 1 is close to 0.41 eV, the experimental value obtained by Arndt and Damask [3] which was

Fig. 6. Activation energy of C atom jumps near a vacant site. Solid arrow shows the direction in which the atom was moved, while broken arrow shows the actual motion direction. All activation energies are found to be equal to ECm ¼ 0:82 eV. The calculations were performed with the EAMtype Fe–Fe potential of Ackland et al. [19].

used in the work by Johnson et al. [7] (see, however, discussion below), while different from 0.85 and 1.1 eV deduced by Vehanen et al. [2] and Takaki et al. [4], respectively. Both F–S and EAM potential representations give similar results. We applied the EAM-Johnson potential set to model statically the migration of a V–C complex. This potential gives 0.63 eV for the vacancy migration energy in a-Fe, which is closer to the experimental estimate of 0.55 eV by Vehanen et al. [2] than the value 0.78 eV with the F–S-type potential by Ackland et al. [18]. The results are presented in Figs. 6–8. As seen from Fig. 6(a), an attempt to move the C atom from its octahedral site to the vacant site resulted in movement to another octahedral site. The same movement occurs when the C atom is displaced towards one of the four nearest octahedral sites, as shown in Fig. 6(b). A successful jump occurs when the C atom is moved directly to one of the equivalent octahedral sites (Fig. 6(c)). The energy of this jump was found to be the same as for an isolated C atom, i.e. 0.82 eV. This is in contradiction with the ab initio calculations [5], which give 0.4 eV for the energy barrier from one octahedral site nearest a vacancy to another.

Fig. 7. Vacancy (V) jumps near a C atom.

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

7

a result the C atom motion is confined to one cell centred on the vacancy position of the maximum binding with the C atom. 4.2. MD simulations

Fig. 8. Energy profiles for the vacancy jumps 1 ! 2, 2 ! 5 and 2 ! 4 (see Fig. 7). The reference state for each curve is the initial configuration before the jump.

Table 4 Activation energy for vacancy jump from site i to site j shown in Fig. 7 and V–C binding energy in different configurations calculated using EAM-type Fe–Fe potential of Ackland et al. [19] Jump type, i!j

1!2 2!3 2!4 2!5 2!6 2!7 2!8 1!9

EV m (eV)

ECV (eV) at site j b

Direct jump (i ! j)

Reverse jump (j ! i)

0.49 0.79 0.63 0.95 0.63 0.63 0.68 1.60

0.28 0.61 0.47 0.73 0.47 0.50 0.48 0.64

+0.17 0.01 +0.01 0.04 +0.01 +0.04 0.02 0.45

The energy barriers for vacancy jumps near a C atom and binding energies of the complex in different configuration are shown in Table 4. As can be seen, the V–C interaction is short-ranged, up to the sites labelled 3 to 8 in Fig. 7, while the jump 1 ! 9 is associated with a high (1.60 eV) energy barrier. The energy profiles during vacancy jumps 1 ! 2, 2 ! 5 and 2 ! 4 are shown in Fig. 8 as examples. The reference state for each curve is the corresponding initial configuration before the jump. Note that for the jump 1 ! 2, there is a local energy minimum at about 0.4a0 distance, which has lower energy than when the jump is complete. For all the sites from 3 to 8 except site 5, the jumps back towards the C atom involve an energy barrier lower than for a jump away (0.63 eV). Hence, we expect a vacancy to return back more frequently to the C atom from these sites than jump away. The ab initio calculations [5] give 0.5 and 0.3 eV for the 1 ! 2 and 2 ! 1 jumps, which agree well with our results. Finally, according to these calculations, a V–C complex cannot migrate without dissociation, since the vacancy cannot go around the C atom and still remain bonded to it. As

This scenario for migration of the V–C pair was verified by MD simulations at temperatures of 800, 900, 1000 and 1200 K. The calculations were performed for 25–40 ns depending on the temperature. The behaviour of the C atom and vacancy was observed to be in accordance with zero temperature relaxation calculations described above. Namely, jumps 1 ! 2 of the vacancy and their reverse 2 ! 1 were the most frequent events. Less frequent were vacancy jumps from site 2 to sites 3 to 8. After these jumps, reverse jumps to the site 2 were most frequently observed, while jumps 1 ! 9 were not observed at all. C atom jumps were observed at high temperature in two situations: when the vacancy was in site 1, the closest position to the C atom, and when the complex was dissociated. Four dissociation events, defined as those when the vacancy position is out of the range shown in Fig. 7, were observed for 43 ns simulated time at 800 K, 5 for 23 ns at 900 K, 18 for 36 ns at 1000 K and 37 for 32 ns at 1200 K. Finally, the calculations confirmed the immobility of the V–C complex unless dissociated. 5. Rate theory analysis of experimental results It follows from the calculations described above that, at temperatures when C atoms are immobile, V–C complexes are also immobile. In addition, the activation energy for dissociation of a V–C complex is equal to 0.63 + 0.51 = 1.14 eV, the sum of the vacancy migration energy in pure Fe matrix and binding energy with a C atom, and thus higher than the migration energy of a C atom. This is consistent with the interpretation of experimental data proposed by Vehanen et al. [2]. Namely, at temperatures corresponding to vacancy migration, vacancies should interact with C atoms and form immobile V–C complexes, while during the stage corresponding to C atom migration, higher-order immobile complexes, presumably V–2C, must form. At higher temperatures, the reverse reactions should take place. However, recent ab initio calculations [5] show that V–2C complexes have stronger binding than V–C pairs. The binding energy of a V–C pair is ECV ¼ 0:44 and b 0.47 eV for supercells of 54 and 128 lattice sites, respectively, while the corresponding values for the binding of a V–C pair with an additional C atom are EC–VC ¼ 0:80 b and 1.03 eV. This suggests that the vacancy recovery stage observed at 480 K should be attributed to dissociation of higher-order complexes, such as V–2C, rather than V–C pairs. We use these data below to verify whether they reproduce correctly the recovery stages observed in Fe–50 appm C alloys irradiated with electrons to a fluence of 1022 e/m2 and heated at a rate of 20 K per 30 min.

8

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

5.1. Equations We consider temperatures above 200 K and assume that our system contains mobile vacancies and C atoms, immobile V–C and V–2C complexes and dislocations. The formation of higher order complexes is neglected. The following set of equations for the concentrations of unpaired (free) vacancies, C 0V , free C atoms, C 0C , V–C complexes, CVC, and V–2C complexes, CV2C, describe the evolution of the system: oC 0V þ ¼ qD DV C 0V þ ðR V–C  RV–C Þ; ot oC VC þ  þ ¼ ðR V–C  RV–C Þ þ ðRC–VC  RC–VC Þ; ot oC V2C þ ¼ ðR C–VC  RC–VC Þ; ot oC 0C þ  þ ¼ ðR V–C  RV–C Þ þ ðRC–VC  RC–VC Þ: ot

ð7aÞ ð7bÞ ð7cÞ ð7dÞ

The total vacancy concentration is C V ¼ C 0V þ C VC þ C V2C :

ð7eÞ

In these equations, qD is the dislocation density, DV is the vacancy diffusion coefficient in pure Fe, and R+ and R are the rates of complex creation and dissociation reactions, respectively. The latter can be expressed through the diffusion properties of mobile species and the binding energy of complexes as

0 0 Rþ V–C ¼ zV–C mV C V C C ; þ tot RC–VC ¼ zC–VC mC C VC C 0C ;

ð8aÞ ð8bÞ

V–C R V–C ¼ mV C VC expðbEb Þ;

ð8cÞ

R C–VC

ð8dÞ

¼

mtot C C V2C

expðbEC–VC Þ; b

2 2 where mtot C ¼ 24DC =a0 and mV ¼ 8DV =a0 are the jump frequencies of C atoms and vacancies in the Fe matrix, respectively, while the z factors account for the details of complex creation/dissociation events. This equation set can be used to study the evolution of the vacancy population at temperatures when interstitials are already removed from the system. Below we justify the choice of parameters entering Eqs. (7) by fitting the numerical calculations to the experimental data. The results of the following analysis are summarised in Table 5, and are discussed below sequentially stage by stage. The calculations for the best fit parameters are also illustrated in Fig. 9 together with experimental measurements [2].

5.2. Stage III: vacancy migration and formation of V–C complexes Stage III is associated with the motion of vacancies. We assume that by this stage the microstructure consists of isolated vacancies, C atoms and dislocations. The sink strength of C atoms for vacancies, i.e. the coefficient before

Table 5 Calculation parameters and results for Fe–50 appm C alloys heated at a rate of 20 K per 30 min DV0 (106 m2/s)

EV m (eV)

zV–C

TIII (K) Integration, Eqs. (7)

Stage III: vacancy migration, formation of V–C complexes (DV0 are the best fit to self-diffusion coefficients [20] in a-Fe at 510 and 600 C for EV f ¼ 2:0 eV) 5.2 0.55 10.0 194 1.0 208 0.2 220 23.0 0.65 10.0 220 1.0 234 EbCV (eV)

qD (1014 m2)

Stage C: carbon migration, formation of V–2C complexes (DC0 = 1.44 · 107 m2/s, EC m ¼ 0:86 eV, zC–VC = 1.0) 0.44 0.17 1.00 0.47 0.45 1.00 0.50 1.00 1.20 EbC–VC (eV)

Analytical, Eq. (10)

Experiment [2]

189 203 214 214 229

220

TC (K)/CV (appm) Integration, Eqs. (7)

Analytical, Eq. (12)

Experiment [2]

350 345 350 346 350 350

345

350 / 3.3

/ / / / / /

3.3 2.0 3.3 2.8 3.4 3.3

TVC (K) Integration, Eqs. (7)

Stage VC: vacancy recovery, dissociation of V–2C complexes (All other parameters are in bold type above) 1.03 – 0.80 500 0.75 487

Analytical, Eq. (17)

Experiment [2]

545 487 474

480–490

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

Fig. 9. Calculation results on the dependence of concentrations of free vacancies and V–C and V–2C complexes during annealing experiments with a heating rate of 20 K per 30 min. The calculation parameters are those in bold type in Table 5.

DV C 0V in the reaction rate described by Eq. (8a), is equal to zV–C 8C C =a20 > 5  1015 m2 (for CC = 50 appm and zV–C P 1), and is thus higher than that of dislocations qD 6 1014 m2. Hence, in this approach, the removal of free vacancies from the system is dominated by the formation of V–C complexes. We note this does not explain the decrease of vacancy concentration observed during stage III, i.e. around 220 K. It might be possible that some immobile or trapped interstitial atoms, in single or clustered form, still remain in the system by this stage and provide necessary sinks for vacancies. We do not study this and, in order to check whether the vacancy diffusion stage is described correctly, start the calculations at a temperature below this stage, using 4 appm for the initial value of vacancy concentration, which was estimated by Vehanen et al. [2] just above 220 K. To make an analytical estimate of the temperature of the stage, we reduce Eqs. (7) to the following single equation: oC 0V ¼ zV–C mV C 0V C C ; ot

ð9Þ

where we take C 0C ¼ C C , since concentration of C atoms is an order of magnitude higher than the vacancy concentration. Eq. (9) has the solution C 0V ¼ C V0 expðzV–C mV C C tÞ, where CV0 is the vacancy concentration before the stage. We define the temperature corresponding to stage III, TIII, as that when the concentration of free vacancies is reduced by e times over time Dt, taken to be equal to 30 min, i.e. that of the heat treatment of a sample in a 20 K temperature interval in the experiments of Vehanen et al. [2]. For this definition, it can be readily obtained that k B T III ¼

EVm ; lnð8zV–C C C DV0 Dt=a20 Þ

ð10Þ

9

where DV0 and EVm are the pre-exponential factor in the diffusion coefficient and the vacancy migration energy, respectively. To perform analytical estimates using Eq. (10) and numerical integrations of Eqs. (7), we need DV0. This was obtained with values EVm ¼ 0:55 and 0.65 eV, as proposed by Vehanen et al. [2] and predicted by ab initio calculations [26], respectively, as the best fit for self-diffuV sion coefficients, DSD in a-Fe of V ¼ DV expðbE f Þ, 22 2 20 2 2 · 10 m /s at 510 C and 10 m /s at 600 C (lowest temperatures in Fig. 13.1 in Ref. [20]). We tried values for the vacancy formation energy, EVf , from 1.8 to 2.3 eV and found that EVf ¼ 2:0 eV, which coincides with the ab initio calculations [26] and experimental value [27], and DV0 shown in Table 5 gives excellent fit to the experimental data. Then, with this choice of DV0, we fitted zV–C to the experimental value for the temperature of 220 K for stage III. The results are shown in Table 5 for the analytical estimate using Eq. (10) and numerical integration of Eqs. (7). It is seen that to reproduce the experiment, the choice EVm ¼ 0:55 eV requires zV–C to be much lower than unity, while, according to the physical significance ofzV–C, it is approximately equal to the number of sites where a complex is bound and must be higher. The best fit was obtained with EVm ¼ 0:65 eV and zV–C = 10, which seems better that 0.55 eV. Thus, the ab initio value EVm ¼ 0:65 eV for the vacancy migration energy reproduces satisfactorily the experimental data for this stage with a reasonable choice of zV–C = 10, EVf ¼ 2:0 eV and DV0 = 2.3 · 105 m2/s, which are also supported by the ab initio calculations [26] and experiment [27]. The results of numerical integration of Eqs. (7) performed with these values are shown in Fig. 9. 5.3. Stage C: carbon migration and formation of V–2C complexes The stage associated with the mobility of C atoms is analysed in a similar way. In this case, the concentration of V–C complexes decreases due to the formation of V–2C complexes according the following equation: oC VC ¼ zC–VC mtot C C VC C C ot

ð11Þ

and the temperature TC of this stage can be estimated via the following equation: kBT C ¼

ECm : lnð24zC–VC C C DC0 Dt=a20 Þ

ð12Þ

To describe this stage, we decided to use DC0 = 1.44 · 107 m2/s and ECm ¼ 0:86 eV. These were obtained in ab initio calculations of Jiang and Carter [17] and are close to experimental values (see Table 2). Then, with zC–VC = 1.0, Eq. (12) predicts TC = 345 K for the temperature when C atoms become mobile, which is in excellent agreement with the observations, while numerical integration of Eqs. (7) gives 350 K for the temperature of this stage if the dislocation density is low enough. The dislocation density enters

10

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

the analysis of this stage because V–C complexes may dissociate, if the binding energy is low enough, and, as a consequence, some fraction of vacancies may annihilate at dislocations. The higher the dislocation density, the lower the vacancy fraction that would remain in the system after this stage. Thus, the fitting parameters for these calculations included the vacancy concentration observed immediately after this stage as reported by Vehanen et al. [2]. The calculation results for different EbCV and qD are shown in Table 5. Note that the ab initio calculations suggest that V–2C complexes are also seen by positrons, as mentioned as a possibility by Vehanen et al. [2]. We further note that, since samples were cold-rolled [2], the dislocation density is expected to be relatively high. This is why we decided that EbCV ¼ 0:47 eV obtained by ab initio calculations for a 128-site supercell is a better choice than 0.44 eV, since it requires higher dislocation density: qD = 4.5 · 1013 m2. The calculation results for this choice of parameters are shown in Fig. 9. 5.4. Stage VC: dissociation of V–2C complexes and complete removal of vacancies from the system Complete removal of vacancies from the system to dislocations observed at TVC  480–490 K requires dissociation of both V–C and V–2C complexes. In order to obtain an analytical expression for the temperature of this stage, we use the detailed balance approximation for the concentrations of the complexes, which requires the rates of complex creation and dissociation to be equal to each other, i.e. R+ = R and hence quasi-steady-state concentrations: oCVC/ot = oCV2C/ot = 0. With this assumption, it is readily obtained that C VC ¼ zV–C C 0V C 0C expðbEV–C b Þ; C V2C ¼

zC–VC C VC C 0C

expðbEC–VC Þ: b

ð13Þ ð14Þ

We note that most of the vacancies are in complexes with two C atoms due to the higher stability of V–2C complexes. Hence the vacancy conservation equation is simplified to C V ¼ C 0V þ C VC þ C V2C  C V2C . Then, with the aid of Eqs. (12) and (13), and taking C 0C  C C , the concentration of free vacancies can be obtained as C 0V 

zV–C zC–VC C 2C

CV : exp½bðEC–VC þ EC–V b b Þ

ð15Þ

The removal of vacancies to dislocations is described by the first term on the right-hand side of Eq. (7a): oC V ¼ qD DV C 0V : ot

ð16Þ

Eqs. (16) and (15) allow us to write an explicit expression for TVC in the following form: k B T VC ¼

EVm þ EbC–V þ EC–VC b : lnðqD DV0 Dt=zV–C zC–VC C 2C Þ

ð17Þ

The calculation results are presented in Table 5 and Fig. 9. ¼ 1:03 eV, the ab initio value for a We see that EC–VC b

128-site supercell [5], overestimates the temperature of this ¼ 0:80 eV, which is for a 54-site superstage, while EC–VC b cell, gives a reasonable value. An even lower value would seem better and EbC–V ¼ 0:75 eV reproduces the stage temperature perfectly. This value also seems reasonable, since typical methodological errors of ab initio calculations are around 0.1 eV. We note that higher binding energy requires lower concentration of C atoms in solution. It is known that C atoms interact strongly with lattice imperfections, such as dislocations and grain boundaries, and so their concentration may be reduced during the C atom migration stage. To fit the temperature of 480 K of the stage with EC–VC ¼ 1:03 eV, an about 20 times reduction of C concenb tration is needed. The dislocations of density qD  1014 m2 with line trap density a1 0 and grain boundaries of radius RGB = 30 lm with surface trap density a2 0 correspond to volume density of traps qD a20 =2  4  106 and 3a0/ 2RGB  105, respectively. These values seem to be too small to provide a sufficient reduction of C concentration. Finally, we conclude that the ab initio data reproduce satisfactorily the experimental data.

6. Conclusions This paper has reported the first simulation study of the dynamics of the behaviour of an interstitial solute and a vacancy–solute complex in Fe. The aim has been to model C in a-Fe, and the simulations have been performed using a combination of two different many-body potentials for the Fe matrix [18,19] and an empirical pair potential of Johnson et al. [7] for the Fe–C interaction. The following conclusions are drawn: (1) The pre-exponential factor in the diffusion coefficient of a C atom is close to the experimental value derived from the low-temperature measurements (below 347 K) [21] and ab initio [5,17] value. The dynamic value of the migration energy of a C atom is lower by 0.1 eV than the static and experimental values. (2) This suggests that the ab initio values [5,17] for the activation energy of C atom migration, which are 0.05–0.1 eV higher than the experimental value, may be consistent after all with experiment if dynamic effects are taken into account. (3) Analysis of C atom jumps shows that it is possible to describe C diffusion as a random walk on the lattice of octahedral sites if an RRT criterion is used for the definition of jumps. (4) MD simulations show that V–C complexes are immobile and have dissociation energy higher than the C atom migration energy. This confirms the interpretation of recovery stages proposed by Vehanen et al. [2] that, during the vacancy migration stage, immobile V–C complexes must form and during the C atom migration stage, higher order complexes, presumably V–2C, must form. However, recent ab initio calcula-

K. Tapasa et al. / Acta Materialia 55 (2007) 1–11

tions [5] indicate that the temperature stage corresponding to complete removal of vacancies from the system is associated with dissociation of higher order complexes rather than V–C pairs. A rate theory analysis of the annealing experiments [2] has been performed. This shows that the results of ab initio calculations [5] for the vacancy migration energy (0.65 eV), the binding energy for a V–C pair (0.47 eV) and the binding energy of an additional C atom with a V–C pair (0.80 eV) reproduce satisfactorily the recovery stages obtained experimentally. Acknowledgements K.T. would like to thank the Science Service Division of the Ministry of Science and Environment, Thailand, for providing a studentship grant. The research was supported by a research grant from the UK Engineering and Physical Sciences Research Council and grant PERFECT (F160CT-2003-508840) under programme EURATOM FP-6 of the European Commission. References [1] Cottrell AH, Bilby BA. Proc Phys Soc Lond A 1949;62:49. [2] Vehanen A, Hautojarvi P, Johansson J, Yli-Kauppila J, Moser P. Phys Rev B 1982;25:762. [3] Arndt RA, Damask AC. Acta Metall 1964;12:341. [4] Takaki S, Fuss J, Kugler H, Dedek U, Schults H. Radiat Eff 1983;79:87.

[5] [6] [7] [8]

[9]

[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

11

Domain C, Becquart CS, Foct J. Phys Rev B 2004;69:144112. Little EA, Harries DR. Met Sci J 1970;4:188. Johnson RA, Dienes GJ, Damask AC. Acta Metall 1964;12:1215. Wert CA. Phys Rev 1950;79:601; Doremus RH. Trans Am Inst Min (Metall) Eng 1960;218:596; Stanley JK. J Metals 1949;1:752; Wagenblast H, Damask AC. J Phys Chem Solids 1962;23:221; Fujita FE, Damask AC. Acta Metall 1964;12:341. Bosman AJ, Brommer PE, Rathenau GH. Physica 1957;23:1001; Bosman AJ, Brommer PE, Eijkelenboom LCH, Schinkel CJ, Rathenau GH. Physica 1960;26:533; Bass J, Lazarus DJ. Phys Chem Solids 1962;23:1820. Smith RP. Trans Metall Soc AIME 1962;224:105. Williamson GK, Smallmann RE. Acta Crystallogr 1953;6:361. Rosato V. Acta Metall 1989;37:2759. Marchese M, Jacucci G, Flynn CP. Philos Mag Lett 1988;57:25. Ruda M, Farkas D, Abriata J. Scripta Mater 2002;46:349. Tersoff J. Phys Rev Lett 1988;61:2879. Lee B-J. Acta Mater 2006;54:701. Jiang DE, Carter EA. Phys Rev B 2003;67:214103. Ackland GJ, Bacon DJ, Calder AF, Harry T. Philos Mag A 1997;75:713. Ackland GJ, Mendelev MI, Srolovitz DJ, Han S, Barashev AV. J Phys Condens Matter 2004;16:S2629. Brandes EA, editor. Smithells metals reference book. London: Butterworth; 1983. Da Silva JRG, McLellan RB. Mater Sci Eng 1976;26:83. Pearson WB. Lattice spacings and structures of metals and alloys. London: Pergamon Press; 1958. p. 919. Osetsky YN. Defects Diffus Forum 2001;188:71. Vineyard GH. J Phys Chem Solids B 1957;3:121. Wert C, Zener C. Phys Rev 1949;76:1169. Becquart CS, Domain C. Phys Rev B 2003;65:024103. DeSchepper L, Segers D, Dorikens-Vanpraet L, Dorikens M, Knuyt G, Stals LM, Moser P. Phys Rev B 1983;27:5257.