Simulation of Diffusion in FCC NiFe Binary Alloys Using Kinetic Monte

Article pubs.acs.org/JPCC

Simulation of Diffusion in FCC NiFe Binary Alloys Using Kinetic Monte Carlo Method Dominic R. Alfonso*,† and De Nyago Tafen†,‡ †

National Energy Technology Laboratory, U.S. Department of Energy, Pittsburgh, Pennsylvania 15236, United States URS Corporation, P.O. Box 1959, Albany, Oregon 97321, United States

‡

ABSTRACT: The use of the atomistic kinetic Monte Carlo method was explored to examine the vacancy-mediated diffusion in fcc NiFe binary alloy. Relevant energetic and kinetic parameters were calculated from density functional theory, and these were used to parametrize a pair interaction model for determination of environment-dependent diffusion barriers. Kinetic Monte Carlo simulations were performed to compute the tracer diffusivities as functions of composition and temperature. Calculations for pure Ni and for Ni with dilute amount of Fe yielded results that compare well with experiments. With increasing amount of Fe, the model predicts a slight reduction in the diffusivity of Fe. This trend was examined on the basis of the neighbor pair interactions between atoms.

I. INTRODUCTION Face-centered-cubic (fcc) nickel (Ni) alloys are highly important metallic materials encountered in a wide range of applications. These materials are prominently employed in the manufacture of engine for commercial and military aircraft and for other challenging environments including power plants.1−3 In view of their high heat resistance and toughness, they are extensively used in the turbine sections of the engine where elevated temperatures are maintained during operation. A specific generation of nickel alloys that includes inexpensive minor iron element was proposed to have suitable properties for many large-scale industrial applications.4 In particular, the combination of low cost, enhanced workability, and high creep strength makes them candidate structural materials for advanced ultrasupercritical power plants. In order to achieve enhanced efficiency in high temperature operations under oxidizing conditions, aluminum (Al) can be added. Its inclusion in the alloy is a prerequisite for the formation of a thermally grown external alumina (Al2O3) layer or scale. Al2O3 is an ideal protective layer since it has the lowest oxygen mobility among all the oxide ceramics. This property coupled with its slow growth rate makes alumina an effective material for protection of the underlying alloy from continued aggressive oxidation.5 The inhibition of further internal oxidation is crucial to the lifetime of the alloy from which the pertinent power plant components are going to be made. To ascertain whether a stable external Al2O3 layer can be achieved and at what rate this can occur, knowledge of the diffusion behavior of atoms is of key importance. The diffusion rates of Al and oxygen, both in the metallic host and within the oxide phase, are expected to strongly influence the growth of the alumina layer.6 Despite this importance, there is relatively little information available in view of experimental limitations © 2014 American Chemical Society

such as lack of suitable isotopes for tracer measurements. Predictive computer simulations offer an alternative route and are well suited, for example, for systematic understanding of fundamental trends related to the chemical composition of the alloy. Fundamental modeling of diffusion in alloys is challenging since this phenomenon largely consists of interplay of ensemble of thermally activated atomic hops occurring for more than microseconds and longer. Such time scales are several orders of magnitude higher than that of atom vibrating around its equilibrium positions. This time scale issue precludes the use of standard molecular dynamics, even in a brute force manner, to simulate diffusion events at realistic temperatures. The kinetic Monte Carlo (KMC) method is a relatively more practical technique for modeling processes that require demanding simulation times. By treating diffusion as stochastic transitions between locally stable states, simulation can be carried out with large system sizes and time scales.7−11 In the past several years, it has become possible to elucidate fundamental diffusion mechanisms from ab initio calculations, yielding important quantities such as migration activation energies and attempt frequencies. As ab initio methods become more widely accessible, there is an increasing interest in applying KMC framework as it can be used to access macroscopic time scale while using kinetic parameters obtained from ab initio calculations to enhance the predictive quality. Herein, we explore the use of such a hybrid approach in which we combine results from ab initio calculations with the KMC method. We implement this approach in the study of Ni Received: June 10, 2014 Revised: August 13, 2014 Published: September 4, 2014 22221

dx.doi.org/10.1021/jp5057607 | J. Phys. Chem. C 2014, 118, 22221−22228

The Journal of Physical Chemistry C

Article

system can be propagated step by step until the target time is reached. Diffusion in typical metallic alloys is mediated by dilute concentration of vacancies.22,23 In this work, metal atoms are assumed to move by exchanging their positions with a vacancy at the nearest-neighbor site in the lattice. The rate ki associated with hop i follows the Arrhenius rate expression24

and Fe diffusion in fcc NiFe alloy. While characterization of the diffusivity of Al and oxygen diffusion in this system is our ultimate goal, looking at Ni and Fe diffusion would be a sensible initial step as it allows us to make initial assessment of the viability of our approach. It should be noted that for a case where Fe is present as an impurity it is possible to calculate diffusion coefficients employing approximate analytical expression founded on dilute diffusion model.12,13 However, such an approach breaks down for a situation where Fe is present in concentrated amount. That is, for realistic alloys where the migration barriers depend on the local environment, the diffusion coefficients have to be evaluated numerically.14 It is also appropriate to mention here previous experimental investigations of diffusion in NiFe alloy. Early on, Shinyaev15,16 and Guiraldenq17 carried out diffusion studies of dilute Fe in Ni using polycrystalline samples. The investigations yielded scattered Fe tracer diffusion coefficient data within the temperature range considered. Subsequent studies based on single Ni crystal using radiotracer method in conjunction with serial sectioning technique were carried out by Bakker et al.18 The utility of these efforts could be more beneficial if they were not confined to limited temperature range. For example, the single crystal measurements of Bakker et al. were limited to temperature range 1473−1673 K.18 Looking at the process at lower temperature, however, is often difficult to achieve from an experimental standpoint. The computational approach adopted here would not only fill these gaps but also could provide atomic level insights into the key aspects of the diffusion process. In what follows, a detailed description of the technical details of the KMC model is given in section 2. We described the atomistic processes and procedure for the determination of the relevant kinetic parameters incorporated into the KMC model. In addition, we discussed the parametrization of the pair interaction model used to represent the environment dependence of the activation barrier. In section 3, simulation results for self-diffusion in pure Ni and dilute Fe diffusion in Ni are first presented for validations purposes. Our investigation of the diffusion in fcc NiFe alloy at various concentration is also included in this section. Finally in section 4, we provide a summary of key findings gleaned from our theoretical work.

ki =

j−1 i=1

+(0) qvib =

(

) 1 − exp( ) hν

exp −0.5 k Ti B

hν − k Ti B

(4)

with νi corresponding to vibrational mode. The term Ei in eq 3 is the activation energy of the hop, which depends on the local environment of the vacancy and diffusing atom, and is given by Bronsted−Evans−Polanyi-type expression25−27 Efinal − E init (5) 2 where E0i is the intrinsic activation energy. Here, it corresponds to the vacancy-mediated diffusion barrier of hopping atom in the pure Ni matrix. The last two terms Einit and Efinal are the system energies before and after the hop. They are determined using pair interactions according to the equation Ei = Ei0 +

E init(final) =

∑ ∑ εjk(i) i

(6)

j>k

ε(i) jk

where is the ith neighbor pair interaction energy for atoms j and k. In this work, pair interactions up to second neighbor were considered. We used the following relationship to determine the pair interaction energy between the same types of atoms ε(i) AA. 2

EAcohesive = −∑ i

ni (i) εAA 2

(7a)

n ⎞ ⎛n (1) εAA = −EAcohesive /⎜ 1 + 2 α⎟ ⎝2 2 ⎠

(7b)

(2) (1) εAA = αεAA

(7c)

where α = 1/2 is adopted. The terms n1 and n2 are the number of first- and second-neighbor which for fcc crystal are n1 = 12 and n2 = 6, while Ecohensive is the cohesive energy of A metal A (= Ni and Fe). The pair interaction energies between Ni and a vacancy, ε(i) NiVac, were determined from the Ni vacancy formation energy 28

(1)

where the overall rate K represents the sum of rates of all possible hopping processes. After the system configuration is adjusted, the time is advanced stochastically according to t → t + Δt where Δt is given by Δt = −ln(σ2)/K

∏ i

j i=1

(3)

where kB, h, and T are the Boltzmann constant, Planck constant, and temperature. The terms q+ and q0 represent the vibrational partition function of the transition and initial states given by

2. COMPUTATIONAL APPROACH 2.1. Kinetic Monte Carlo Simulation. The KMC method implemented here is based on the Bortz−Kalos−Lebowitz (BKL) algorithm.19−21 In this scheme, the simulation cell is scanned, and a list containing all possible hopping processes and their associated rates ki at the current system configuration is created. A random number σ1 ∈ [0,1] is generated, and the hop j that is selected satisfies the following condition

∑ ki ≤ σ1K < ∑ ki

+ ⎛ E ⎞ kBT qvib exp⎜ − i ⎟ 0 h qvib ⎝ kBT ⎠

2 Vac E Ni =

∑− i

(2)

The term σ2 is another random number between 0 and 1. From the updated configuration, a new process table is recreated and the whole procedure starts again. Following this scheme, the 22222

ni (i) (i) εNiNi + niεNiVac 2

(8a)

n (1) n (2) ⎞ ⎛ Vac (1) ⎟ /(n + βn ) εNiVac = ⎜E Ni + 1 εNiNi + 2 εNiNi 2 ⎝ ⎠ 1 2 2

(8b)

(2) (1) εNiVac = βεNiVac

(8c)



Article

information. In this work, five movable images were used, and an initial chain of images was constructed between the initial and final states using linear interpolation between the two end points. The transition state was approximated by the image of highest energy. In the calculations, the coordinates of the Ni matrix and diffusing atom were allowed to relax. The vibrational frequencies, νi, in eq 4 were calculated by fixing the host atoms in their relaxed positions and computing the normal-mode frequencies of the diffusing atom (all the 3N degrees of freedom were considered) using the harmonic oscillator approximation. Experimental cohesive energy for Ni was used to fit the values of ε(i) NiNi since the one calculated from DFT is less precise. In the case of ε(i) FeFe, the rescaled cohesive energy for the fcc structure was employed which was determined by taking the DFT predicted energy difference between the fcc and the most stable bcc phase and then adding it to the experimental bcc phase cohesive energy. The vacancy formation energy and the binding energy that were used to parametrize the pair interaction energies between a metal atom and a vacancy and between two types of metal atoms were obtained from DFT calculations. The vacancy formation energy EVac Ni in eq 8 is calculated by the following expression

where β = 1/228 and EVac Ni is the Ni vacancy formation energy. The pair interactions between solute Fe and a vacancy ε(i) FeVac and between Ni and solute Fe ε(i) NiFe were estimated using the relation28 b(i) (i) (i) (i) (i) EXY = εNi Y + εNiX − εXY − εNiNi ,

i = 1, 2

(9a)

Eb(i) XY

where is the binding energy between two entities X and Y in the Ni host. That is, the binding energy between Fe atoms in the Ni host together with eq 9a was used to determine ε(i) NiFe (i) b(i) (i) (i) εNiFe = (E FeFe + εFeFe + εNiNi )/2

(9b)

while eq 9a, the binding energy between Fe and vacancy, and (i) ε(i) NiFe from eq 9b were used to obtain εFeVac (i) b(i) (i) (i) (i) εFeVac = −E FeVac + εNiVac + εNiFe − εNiNi ,

i = 1, 2 (9c)

The evaluation of Ni vacancy formation energy and binding energies are discussed below. Using the theory of random walk, the tracer diffusion coefficient which is a useful quantity to characterized mobility can be extracted from KMC simulations through ∞

D* = lim

t →∞

1 ⟨∑ [ r (⃗ t + Δt ) − r (⃗ t )]2 ⟩ 6t t = 0

Vac E Ni = EA [(N − 1) + Vac] −

(10)

In eq 10, t is the simulation time and the vector r(⃗ t + Δt) − r(⃗ t) represents the change in position after a time interval Δt. The angular brackets indicate an expectation value which is obtained by averaging over long-time intervals. In this work, 5 × 108 KMC steps were performed using 5000 sampling passes to collect data on the trajectory of the atoms. 2.2. Ab Initio Computational Details. First-principles spin-polarized DFT calculations were carried out as implemented in the Vienna Ab Initio Simulation Package (VASP) code.29 This implementation includes total energy and atomic force calculations. The local density approximation (LDA) was used to calculate the exchange-correlation energy.30 The electron−ion interaction was described by the projectoraugmented wave (PAW) method.31 The Kohn−Sham one electron valence eigenstates were expanded in terms of planewave basis sets with a cutoff energy of 400 eV. The k-point sampling of the three-dimensional electronic Brillouin zone of the periodic supercells was performed using the Monkhorst− Pack scheme.32 In the determination of the intrinsic activation energy in eq 5, the host Ni lattice was represented by a 3 × 3 × 3 fcc supercell of 108 sites with a 2 × 2 × 2 k-point grid. For selfdiffusion in the Ni lattice, we consider the exchange of a Ni atom with first-neighbor vacancy along the ⟨110⟩ direction. For Fe diffusion, a Ni atom in the lattice was replaced by a Fe atom, and a similar hop to first-neighbor vacancy was examined. The strategy in locating the transition states and calculating the corresponding activation energy consisted of two steps. First, we identified the process local minima (i.e., initial and final states), and then we determined the transition state configuration along the path that the system is likely to follow to its transition from one of the identified local minimum to another. The climbing image nudged elastic band (CI-NEB) method for finding saddle points and minimum energy paths between known two states was employed.33 A discrete representation of the path was generated, with the points (movable images) along the path being relaxed using first derivative based

N−1 E Ni[N ] N

(11)

where ENi[(N − 1) + Vac]] is the total energy for supercell containing (N − 1) atoms and a vacancy while ENi[N] is the total energy for perfect bulk with N atoms. The binding energy between ith neighbor X and Y in the Ni host is determined using b(i) EXY = E Ni[(N − 2)XY ∞] − E Ni[(N − 2)XY (i)]

(12)

where ENi[(N − 2)XY∞] and ENi[(N − 2)XY(i)] are the energies of Ni host containing two infinitely separated X and Y and ith neighbor X and Y. We used a 3 × 3 × 4 fcc supercell with a 2 × 2 × 2 k-point grid to calculate these two terms with two infinitely separated X and Y approximated by eight neighbor XY. As indicated previously, the activation energies for hopping are approximated by the Bronsted−Evans−Polanyi-type expression given in eq 5 whose predictive quality needs to be ascertained. To roughly assess its validity, additional calculations were carried out. The barriers of 14 different vacancy mediated diffusion events consisting of atomic Fe or Ni hopping in the vicinity of 1−3 Fe atoms were calculated by DFT. The host Ni lattice was represented by a 3 × 3 × 3 fcc supercell of 108 lattice sites. The barriers were determined using the CI-NEB algorithm as described previously. A correlation plot between the DFT calculated barriers and those from eq 5 is given in Figure 1. The level of agreement was found to be very good. On the average, the estimated barriers from eq 5 are within ∼3% of the DFT calculated ones.

3. RESULTS AND DISCUSSION 3.1. Ab Initio Benchmark Calculations. Bulk Ni adopts an fcc structure, and its equilibrium lattice parameter was obtained by calculating the energy of unit cells with lattice constants of ∼±10% of the experimental value. They were then fit to the Murnaghan’s equation of state34 yielding a bulk lattice constant of a = 3.42 Å, which is in good agreement with the value aexp = 3.52 Å found in the literature.35 The theoretical 22223



Article

Similarly for Fe diffusion in the Ni fcc lattice, the hopping of the solute to an adjacent vacancy was considered. The movement of the Fe atom toward the first-neighbor vacant site is essentially similar to that in Ni except that the hopping barrier is found to be slightly lower by 0.09 eV (Figure 2b). Our calculations concur with previous DFT works12,36 which also predicted a small decrease in Fe hopping barrier (∼0.1 eV) in comparison to Ni. They argued that the compressibility of elemental Fe provides a physical basis for understanding this difference. That is, the slightly bigger Fe forms compressible Fe−Ni host directional bonds, stabilizing the transition state and resulting into the lowering of the barrier. 3.2. Self-Diffusion of Ni. In the remainder of the paper, the KMC approach described in section 2 was used for the study of Ni and Fe diffusion in fcc NiFe alloy. The computed kinetic parameters adopted in the simulations are tabulated in Table 1.

Figure 1. Correlation plot between DFT calculated barriers and those obtained from eq 5. Diffusion event represented by diamond, square, and triangle corresponds to diffusion of either Ni or Fe in the vicinity of one, two, and three Fe atoms, respectively.

Table 1. Calculated Kinetic Parameters and Pair Interaction Database for Diffusion Process in Ni Host

value was found to be within ∼3% of the measured one. For Ni in the ferromagnetic state, the calculated local magnetic moment is 0.58 μB, which is in very good agreement with the experimental value of 0.61 μB.35 In a 3 × 3 × 3 108 atom Ni supercell, a Fe atom introduced in a substitutional site distorted the first-neighbor Ni atoms by a very small amount (0.004 Å). The calculated local magnetic moment of Fe is 2.49 μB. To calculate the vacancy formation energy in pure Ni, a Ni atom was removed from a 3 × 3 × 3 supercell. Upon relaxation, the Ni−Ni bonds in the first-neighbor shell around the vacancy are slightly reduced by about 1.5%. The vacancy formation energy was predicted to be EVac Ni = 1.63 eV. DFT calculations by Krcmar et al.12 and Janotti et al.36 based on 64-atom supercell and using the same LDA functional yielded values of 1.71 and 1.70 eV, which compare well with our prediction. The minimum-energy pathway (MEP) of Ni hop to the firstneighbor vacant site in the fcc lattice was determined, and the results are summarized in Figure 2. A 3 × 3 × 3 supercell was also used to study this mechanism. The diffusing Ni atom moves roughly linearly from the initial state to the final state recreating a vacancy. The saddle point occurs when the atom is about halfway to the target vacant site. Examination of the MEP in Figure 2b indicates that the diffusion barrier is 1.22 eV. A similar DFT-LDA calculation by Krcmar et al.12 using a 64-atom supercell gave an estimated value of 1.15 eV, which compares well with the prediction here.

E0i (eV)

process Ni + Vac → Vac + Ni Fe + Vac → Vac + Fe pair interaction first neighbor (eV) ε(i) NiNi ε(i) FeFe ε(i) NiFe ε(i) NiVac ε(i) FeVac

−0.59 −0.55 −0.62 −0.18 −0.20

1.22 1.13 second neighbor (eV) −0.30 −0.28 −0.28 −0.09 −0.08

from eqs 7b, 7c eqs 7b, 7c eq 9b eqs 8b, 8c eq 9c

We first applied it to the case of self- diffusion in Ni. The selfdiffusion coefficients of Ni were evaluated for temperature range T = 700−1700 K. The Ni lattice was modeled using a cubic box with linear dimension L = Na, where a is the lattice parameter and N = 10. This corresponds to 4000 lattice sites in the simulation cell. Periodic boundary conditions were then applied to simulate the infinite crystal. A vacant site was introduced in the lattice to promote the hopping of the atoms. One vacancy in the 10 × 10 × 10 simulation cell corresponds to a higher vacancy concentration compared to the experimental one at temperature T. To make a direct comparison with experiments, the predicted diffusion coefficients were scaled according to the expression37,38 D=

C Vac,real D* * C Vac

(13)

Figure 2. (a) Hopping pathway and (b) predicted MEP of Ni and Fe diffusion to first-neighbor vacant site in host Ni. 22224



Article

where C*Vac is the vacancy concentration in the simulation box and CVac,real = exp(−ΔGVac,real/kBT) is the real equilibrium vacancy concentration of Ni.39 The term kB is the Boltzmann constant while ΔGVac,real corresponds to the Gibbs free energy of vacancy formation in Ni given by ΔGVac,real = ΔGwith vacancy − ((N − 1)/N)Gno vacancy. The terms Gwith vacancy, Gno vacancy, and N refer to the free energy of the Ni crystal with vacancy, free energy of the perfect Ni crystal and the number of atoms in the supercell. The Gibbs free energy is defined as G = Eel + EZPE + Uvib − TSvib with Eel being the electronic energy calculated from DFT. The vibrational contributions to the internal energy, EZPE and Uvib, are defined as EZPE = 0.5∑ihvi and Uvib = R∑i[(hvi/kB)/(ehvi/kBT − 1)] while the entropic term is given by Svib = R∑i[(hvi/kBT)/(ehvi/kBT −1) − ln(1 − e−hvi/kBT)]. Configurational entropy was not considered since the resulting vacancy concentration is very small (vacancy fraction is ∼10−11− 10−4 in the T = 700−1700 K range). These three terms are calculated using DFT-predicted vibrational frequencies vi. The predicted Ni self-diffusion coefficients at different temperatures are shown in Figure 3 (labeled as KMC). The

As seen in Figure 3, each set of experimental data lies on a straight line for the temperature range considered showing that Ni diffusivity is an Arrhenius function of temperature as well. The calculated diffusion coefficients are in reasonable agreement with experiments with difference of about an order of magnitude. The difference could be due to the calculated equilibrium vacancy concentration of Ni, CVac,real, that was used for scaling (see eq 13). The calculated CVac,real in logarithmic scale as a function of 1000/T was plotted in and shown in the inset of Figure 3. Included in this plot are the experimental values of Scholz previously determined using differential dilatometry measurements.47 It can be seen that the DFTLDA calculations noticeably underestimates the vacancy concentration in comparison with experiment. According to our calculations, the entropic term in ΔGVac,real contributes by only ∼3 × 10−4 eV/K in the 700−1700 K temperature range. As noted in previous studies, the discrepancy between theory and experiment is due to the inaccurate estimation of the vacancy formation energy (the dominant component of ΔGVac,real) attributed to the shortcoming of the uniform density based assumption of LDA to describe vacancy internal surface.48 That is, the uniform density based assumption of DFT breaks down in terms of describing a material with vacancy since it introduces a steep variation of electronic density near the vacant site. It has been noted that even the use of slowly varying density based generalized gradient approximation (GGA) does not address this issue as LDA tends to provide a better description of internal surface due to the wellknown cancellation effect.48−50 The KMC predicted self-diffusion coefficients based on the experimental vacancy concentration are plotted in Figure 3 (labeled KMC′). For reference, the self-diffusion coefficients from the Einstein equation were also calculated under such condition. Again, our KMC prediction is in excellent agreement with the values obtained from such approach. In this case, an improvement in the agreement with experimental data was found. For each set of data, the diffusion prefactor (D0) and the effective activation energy (QNi) can be directly evaluated by fitting to the Arrhenius relation DNi = D0 exp(−QNi/kBT). The results are summarized in Table 2. The values obtained from

Figure 3. Self-diffusion coefficients in fcc Ni computed using kinetic Monte Carlo plotted versus 1000/T. Single crystal experimental data [Bakker (1968),40 Maier (1976),41 Feller (1976),42 Wazzan (1965),43 Ivantsov (1966),44 and McEwan (1953)45] and values obtained from the Einstein equation13,46 are included for comparison. DFT-LDA and experimentally predicted vacancy concentration in Ni plotted as a function of 1000/T are shown in the inset. Values labeled by KMC′ and Einstein Eq.′ are calculated based on experimental vacancy concentration (see text for more details).

Table 2. Computed and Experimentally Measured [Bakker,40 Feller,42 Ivantsov,44 McEwan,45 Maier,41 and Wazzan43] Diffusion Prefactor D0 and Effective Activation Energy QNi for Self-Diffusion in Pure Ni D0 (m2/s) KMC KMC′ Bakker (1968) Feller (1976) Ivantsov (1966) McEwan (1953) Maier (1976) Wazzan (1965) Einstein Eq Einstein Eq′

diffusion coefficients are plotted on a logarithmic scale as a function of 1000/T. The fact that the data lie on a straight line shows that the Ni diffusivity is an Arrhenius function of the temperature. For reference, we calculated the self-diffusion with the Einstein equation13,46 using DFT-LDA calculated data: DNi = fa2Cvac,realΓ, where f, a, and Cvac,real are the jump correlation factor (= 0.78 for pure fcc), lattice constant of Ni, and the real vacancy concentration in Ni while the term Γ is the rate of Ni-vacancy exchange defined as Γ = kBT/h* exp(ΔGbarrier/ kBT) with ΔGbarrier defined as the free energy diffusion barrier. The results (labeled Einstein Eq.′) are displayed in Figure 3. The KMC prediction is found to be in excellent agreement with the values obtained from such an approach. We now compare our KMC prediction with existing experimental diffusion data based on single crystal sample.40−45

2.46 3.24 1.17 2.28 2.37 2.30 1.29 1.37 2.28 2.96

× × × × × × × × × ×

10−5 10−5 10−4 10−4 10−4 10−4 10−4 10−4 10−5 10−5

QNi (eV)

T (K)

2.90 2.77 2.91 2.87 3.02 2.98 2.91 2.88 2.90 2.77

700−1700 700−1700 1250−1650 1100−1270 1170−1470 1420−1675 880−1195 750−925 700−1700 700−1700

KMC and the Einstein equation compare extremely well. Although the calculated diffusion coefficients agree very well with experiments, we found noticeable discrepancies in the values of D0 and QNi. One possible explanation for this is that each set of experimental self-diffusion coefficient data span 22225



Article

limited temperature range, causing errors in the fit of the Arrhenius plots. Experiments measurements over far wider temperature range with consistent methodology is clearly needed to reach consensus. Another source of error could be due to the experimental CVac,real we used for scaling which is limited to high temperature range. To obtain the values at lower temperatures, extrapolation was done by fitting it to an Arrheniuslike function. The scatter in the experimental data and the limited temperature range covered could give rise to errors in the fit. Spot checks were carried out in order to ascertain the sensitivity of our results with respect to the size of the simulation box. The diffusion coefficients for selected temperatures of T = 1000, 1200, and 1400 K were determined using a larger 12 × 12 × 12 supercell. For this larger supercell, the KMC simulation yielded values of DNi = 3.39 × 10−19, 7.16 × 10−17, and 3.30 × 10−15 m2/s, respectively. The values for the 10 × 10 × 10 supercell are DNi = 3.36 × 10−19, 7.24 × 10−17, and 3.33 × 10−15 m2/s, which compare very well with the above results. 3.3. Fe Diffusion in Ni. We used the same approach to examine the tracer diffusivity of Fe in Ni. For reference we first examined the case where Fe is present in dilute amount. A vacant site and about 0.1 atomic % of Fe were introduced in the 10 × 10 × 10 simulation cell which was then equilibrated before sampling was carried out. Tracer diffusion coefficients for temperatures in the range 700−1700 K were calculated. Since Fe is present in dilute amount, we set the vacancy concentration of the real system to that of pure Ni and the predicted diffusivities were scaled accordingly to make direct comparison with experiments. The calculated variation of DFe with inverse temperature is plotted in Figure 4. For comparison, the predicted values corresponding to DFT-LDA and experimental vacancy concentration are shown (labeled KMC and KMC′, respectively). Included in Figure 4 for comparison is the experimental data by Bakker et al. based on single-crystal measurements obtained for

temperatures ranging from 1000 to 1700 K.18 For reference, data obtained from polycrystalline samples are included.15−17 Similar to the single crystal case, very limited data points were reported with no data available for lower temperatures. DFe calculated with the five-frequency model version of the Einstein equation using DFT-LDA calculated energetics and vacancy concentration was reported by Krcmar et al.12 For comparison, these data are included in Figure 4. The calculated and experimental QFe and D0 values are tabulated in Table 3. Table 3. Calculated and Experimentally Measured Single Crystal [Bakker (1971)18] and Polycrystalline Sample [Shinyaev (1958), Shinyaev (1970), and Guiraldenq (1962)15−17] Diffusion Prefactor D0 and Effective Activation Energy QFe for Dilute Fe in Ni KMC KMC′ Bakker (1971) Shinyaev (1958) Shinyaev (1970) Guiraldenq (1962) Krcmar (2005)

D0 (m2/s)

QFe (eV)

T (K)

1.35 × 10−5 1.75 × 10−5 −4 1.0+0.5 −0.4 × 10 9.0 × 10−7 7.4 × 10−6 8.0 × 10−5 1.10 × 10−5

2.84 2.71 2.79 ± 0.6 2.2 2.54 2.65 2.85

700−1700 700−1700 1478−1669 1223−1523 1293−1536 1173−1423 1173−1573

Our KMC results corresponding to the theoretical vacancy concentration gave DFe, QFe, and D0 values that are in excellent agreement with the values of Krcmar et al., which was also calculated using theoretical vacancy concentration. With respect to the single-crystal experimental data of Bakker et al., the agreement is found to be reasonable. The predicted diffusivity is roughly an order of magnitude lower and the QFe and D0 values of 2.84 eV and 1.35 × 10−5 m2/s compare fairly well with the reported experimental values of 2.79 ± 0.6 eV and 1.0 × 10−4 m2/s. Similar to the pure Ni case, the matching of predicted diffusivity with experiments is better when we considered the experimental vacancy concentration. That is, the values predicted with this approach are only ∼3 times lower, with QFe and D0 values of 2.71 eV and 1.75 × 10−5 m2/s, respectively. Again, a potential source of this variance is that the experimental data covers limited temperature range (at high temperature), causing statistical errors in the fit. The presence of dilute Fe has a very small effect on the diffusivity of Ni. The predicted diffusion coefficients of Ni for this system are within ≤3% of the pure Ni values. Additionally, the QNi (2.77 eV) and D0 (1.75 × 10−5 m2/s) values are essentially similar to the pure Ni case. Our results provide theoretical confirmation to the experimental finding of Bakker et al. that the presence of Fe impurity has a slight influence on Ni selfdiffusion.18 Comparison of predicted QFe and QNi indicates that the former is slightly lower (2.71 vs 2.77 eV). It should be noted that the diffusivity of Fe is also strongly influenced by the availability of first-neighbor vacancy. Inspection of Table1 (1) indicates that the values of ε(1) FeVac and εNiVac are comparable. This indicates that it is equally favorable from energetic standpoint for Fe and Ni to be sitting next to a vacant site. This plus the fact that the DFT calculated barrier for Fe hop is lower by ∼0.1 eV (Figure 1b) could explain why the effective activation energy of Fe is noticeably smaller. Our calculations were then extended to examine the case of nondilute Fe diffusion in Ni. Four compositions ranging from 1 up to 30% (i.e., Ni containing x = 0.01, 0.1, 0.2, and 0.3 atomic fraction of Fe which are alternatively referred to as Ni0.99Fe0.01,

Figure 4. Diffusion coefficients of dilute Fe in fcc Ni computed using kinetic Monte Carlo plotted versus 1000/T. Experimental data based on single crystal [Bakker (1971)18] and polycrystalline samples [Shinyaev (1958),15 Shinyaev (1970),16 and Guiraldenq (1962)17] including the calculated data based on the five-frequency model version of the Einstein equation [Krcmar (2005)]12 are included for comparison. Values labeled by KMC′ are calculated based on experimental vacancy concentration (see text for more details). 22226



Article

Ni0.9Fe0.1, Ni0.8Fe0.2, and Ni0.7Fe0.3) and temperatures in the interval 700−1700 K were chosen. Fe atoms were randomly distributed in the lattice sites of the supercell, with each configuration equilibrated using KMC as a preparatory step. The vacancy formation energy in the binary alloy which characterizes the vacancy concentration can be obtained within a pairwise approximation using the analytical expression51

has important consequence for the mobility of Fe in the presence of other Fe atoms. The interactions between Fe atoms are relatively less favorable and thus the diffusion of Fe occurs with energy penalty. That is, the generally slightly lower diffusivity observed for Fe with increasing Fe originates from reduced probability for Fe-vacancy exchange to occur.

Vac Vac Vac E NiFe = x NiE Ni + x FeE Fe + x Nix Fe ∑ zivi

4. CONCLUSIONS An atomistic kinetic Monte Carlo framework was developed in order to examine vacancy mediated diffusion process in fcc Ni−Fe alloy. It includes an efficient neighbor interaction based model for efficient determination of local environment dependent activation energies. The model was parametrized using relevant atom−atom interaction, atom−vacancy pair interaction, and Ni and Fe diffusion barriers determined from firstprinciples density functional theory calculations. For validation purposes, the framework was used to simulate diffusion in pure Ni and in Ni with dilute concentration of Fe. The KMC determined temperature dependent self-diffusion coefficients were compared to those determined from experiments. The calculations yielded values that are within a factor of approximately 10 with respect to experiments. The results indicate that the agreement improves if experimentally known vacancy concentration is utilized in the model. The tracer diffusion coefficient of Fe slightly decreases with increasing concentration of Fe. More insight about this trend can be gained from the analysis of the neighbor interactions between atoms. It arises from the slightly stronger Ni−Fe and Ni−Ni interactions compared to its Fe−Fe counterpart. That is the less energetically favorable Fe−Fe neighbor interactions result in reduced probability for Fe-vacancy exchange to occur as Fe concentration increases.

i Vac where xNi (xFe) is the atomic fraction of Ni (Fe), EVac Ni (EFe ) is the vacancy formation energy of Ni (Fe), zi is the coordination number for the ith nearest neighbors (z1 = 12, z2 = 6), and vi is the ordering energy for the ith neighbor defined as vi = (i) (i) ε(i) NiNi + εFeFe − 2εNiFe. Based on this expression, the computed vacancy formation energy for the various alloys considered here are found to be within ∼3% of the vacancy formation energy of pure Ni. Thus, the vacancy concentration of the real system was set to that in pure Ni. Shown in Figure 5 is the plot of the temperature dependence of DFe for the various alloys. Generally, we find a small change

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected]; Tel +1 (412) 386-4113 (D.R.A.). Notes

The authors declare no competing financial interest.

■

Figure 5. KMC predicted diffusion coefficients of Fe in fcc Ni at various concentration x plotted versus 1000/T. For reference, the composition dependent values of QFe are shown in the inset.

ACKNOWLEDGMENTS This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of author(s) expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

in the Fe diffusivity with increasing Fe content. Compared with the dilute Fe case, the values for Ni0.99Fe0.01 are reduced by a factor of up to 0.16 and 0.07 for T = 700−900 K and T = 900− 1700 K. For Ni0.9Fe0.1, Ni0.8Fe0.2, and Ni0.7Fe0.3, the reduction factor is 0.23−0.77 for T = 700−900 K and up to 0.22 for T = 900−1700 K. This effect is also reflected in the composition dependence of predicted QFe which is shown in the inset of Figure 5. It is found that increase in the Fe content generally induces small incremental increase in the effective activation energy of Fe. This demonstrates the influence of pair interactions and the above-mentioned behavior can be explained by a comparison (i) (i) of the calculated values of ε(i) NiNi, εNFe, and εFeFe in Table 1. The first-neighbor attraction between Ni and Fe (ε(1) NFe= −0.62 eV) is larger than between mutual Fe (ε(1) FeFe = −0.55 eV) while the (2) second-neighbor ones are equal ε(2) NFe = εFeFe = −0.28 eV). (1) Interactions between Ni atoms (εNiNi = −0.59 and ε(2) NiNi = −0.30 eV) are also more favorable than between Fe atoms. This

■

REFERENCES

(1) Pollock, T. M.; Tin, S. Nickel-Based Superalloys for Advanced Turbine Engines: Chemistry, Microstructure and Properties. J. Propul. Power 2006, 22, 361−374.

22227



Article

(2) Sims, S. T.; Stollof, N. S.; Hagel, W. C. SuperAlloys II; Wiley: New York, 1986. (3) Durand-Charre, M. The Microstructure of Superalloys; Gordon and Breach: Amsterdam, 1997. (4) Zhong, Z. H.; Gu, Y. F.; Yuan, Y.; Shi, Z. A New Wrought Ni-Fe Based SuperAlloy for Advanced Ultra-Supercritical Power Plant Applications Beyond 700 °C. Mater. Lett. 2013, 109, 38−41. (5) Heauer, A. H.; Hovis, D. V.; Smialek, J. L.; Gleeson, B. Alumina Scale Formation: A New Perspective. J. Am. Ceram. Soc. 2011, 94, S146−s153. (6) Brady, M. P.; Gleeson, B.; Wright, I. G. Alloy Design Strategies for Promoting Oxide Scale Formation. JOM 2000, 52, 16−21. (7) Fichthorn, K. A.; Scheffler, M. Island Nucleation in Thin-Film Epitaxy: A First-Principles Investigations. Phys. Rev. Lett. 2000, 84, 5371−5374. (8) Bogicevic, A.; Ovesson, S.; Hyldgaard, P.; Lundqvist, B. I.; Brune, H.; Jennsion, D. R. Nature, Strength and Consequences of Indirect Adsorbate Interactions on Metals. Phys. Rev. Lett. 2000, 85, 1910− 1913. (9) Negulyaev, N. N.; Stepanyuk, V. S.; Hergert, W.; Fangohr, H.; Bruno, P. Effect of the Long-Range Interactions on the Atomic SelfAssembly on Metal Surfaces. Surf. Sci. 2006, 600, L58−L61. (10) Sendner, C.; Sakong, S.; Gross, A. Kinetic Monte Carlo Simulations of the Partial Oxidation of Methanol on Oxygen-Covered Surface. Surf. Sci. 2006, 600, 3258−3265. (11) Alfonso, D. R. Influence of Sulfur Poisoning on CO Adsorption on Pd(100). Top. Catal. 2012, 55, 267−279. (12) Krcmar, M.; Fu, C. L.; Janotti, A.; Reed, R. C. Diffusion Rates of 3d Transition Metal Solutes in Nickel by First-Principles Calculations. Acta Mater. 2005, 53, 2369−2376. (13) LeClaire, A. D. Self Diffusion in Dilute Alloys. J. Nucl. Mater. 1978, 69−7, 70−96. (14) Swoboda, B.; Van der Ven, A.; Morgan, D. Assessing Concentration Dependence of FCC Metal Alloy Diffusion Coefficients Using Kinetic Monte Carlo. J. Phase Equilib. Diffus. 2010, 31, 250− 259. (15) Shinyaev, A. Y. Diffusion in Nickel Based Limited Solid Solutions. Fiz. Metal. Metalloved. Akad. Nauk SSSR 1958, 6, 68. (16) Shinyaev, A. Y. Diffusion in Alloys of the Ni-Fe System. Izv. Akad. Nauk SSSR, Met. 1969, 4, 182−186. (17) Guiraldenq, P. Metallographie - Influence des Impuretes et de la Stucture Fritee sur les Coefficients de Diffusion en Volume et aux Joints des Grains du Fer dans le Nickel. C. R. Acad. Sci. 1962, 254, 1994−1996. (18) Bakker, H.; Backus, J.; Waals, F. A Curvature in the Arrhenius Plot for the Diffusion of Iron in Single Crystals of Nickel in the Temperature Range from 1200 to 1400 °C. Phys. Status Solidi B 1971, 45, 633−638. (19) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. New Algorithm for Monte Carlo Simulation of Ising Spin Systems. J. Comput. Phys. 1975, 17, 10−18. (20) Voter, A. F. Introduction to the Kinetic Monte Carlo Method. In Radiation Effects in Solids; Springer: Dordrecht, Netherlands, 2007; Vol. 235, pp 1−23. (21) Battaile, C. C. The Kinetic Monte Carlo Method: Foundation, Implementation and Application. Comput. Methods Appl. Mech. Eng. 2008, 197, 3386−3398. (22) Van der Ven, A.; Ceder, G.; Asta, M.; Tepesch, P. D. FirstPrinciples Theory of Ionic Diffusion with Nondilute Carriers. Phys. Rev. B 2001, 64, 184307. (23) Van der Ven, A.; Thomas, J. C.; Xu, Q.; Swoboda, B.; Morgan, D. Nondilute Diffusion from First Principles: Li Diffusion in LixTiS2. Phys. Rev. B 2008, 78, 104306. (24) Michalak, W. D.; Miller, J. B.; Alfonso, D. R.; Gellman, A. J. Uptake, Transport and Release of Hydrogen from Pd(100). Surf. Sci. 2012, 606, 146−155. (25) Simonovic, D.; Ande, K. Diffusion of Carbon in BCC Fe in the Presence of Si. Phys. Rev. B 2010, 81, 054116−1−054116−7.

(26) Lee, E.; Prinz, F. B.; Cai, W. Enhancing Ionic Conductivity of Bulk Single-Crystal Ytria-Stabilized Zirconia by Tailoring Dopant Distribution. Phys. Rev. B 2011, 83, 052301−1−052301−4. (27) Boreskov, G. K. Heterogeneous Catalysis; Nova Science Publishers Inc.: New York, 2003. (28) Soisson, F.; Fu, C. C. Cu-Precipitation in Alpha-Fe from Atomistic Simulations: Vacancy-Trapping Effects and Cu-Cluster Mobility. Phys. Rev. B 2007, B76, 214102. (29) Kresse, G.; Furthmüller, J. Effeciency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (30) Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, A. T.; Joannopoulus, J. D. Iterative Minimization Techniques for Ab Initio Total-Energy Calculations: Molecular Dynamics and Conjugate Gradients. Rev. Mod. Phys. 1992, 64, 1045−1097. (31) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758. (32) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (33) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (34) Murnaghan, F. D. The Compressibility of Media under Extreme Pressures. Proc. Natl. Acad. Sci. U. S. A. 1944, 30, 244−247. (35) Kittel, C. Introduction to Solid State Physics, 6th ed.; Wiley: New York, 1986. (36) Janotti, A.; Krcmar, M.; Fu, C. L.; Reed, R. C. Solute Diffusion in Metals: Larger Atoms Can Move Faster. Phys. Rev. Lett. 2004, 92, 085901-1−085901-4. (37) Terentyev, D. Molecular Dynamics Study of Oxygen Transport and Thermal Properties of Mixed Oxide Fuel Cells. Comput. Mater. Sci. 2007, 40, 319−326. (38) Shewmon, G. Diffusion in Solids; McGraw-Hill: New York, 1963. (39) Mehrer, H. Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion Controlled Process; Springer-Verlag: Berlin, 2007. (40) Bakker, H. A Curvature in the ln D Versus 1/T Plot for Self Diffusion in Nickel at Temperatures from 980 to 1400 C. Phys. Status Solidi B 1968, 28, 569−576. (41) Maier, K.; Mehrer, H.; Lessmann, E.; Schule, W. Self-Diffusion in Ni at Low Temperatures. Phys. Status Solidi B 1976, 78, 689−698. (42) Feller, M.; Grundler, M.; Hefmeier, H. Self-Diffusion in HighPurity Ni Single Crystals and Ni-0.02 At-Percent Au Single Crystal. Z. Metallkd. 1976, 67, 533−537. (43) Wazzan, A. R. Lattice and Grain-Boundary Diffusion in Nickel. J. Appl. Phys. 1965, 36, 3596−3599. (44) Ivantsov, I. G. Self-Diffusion in Monocrystalline Nickel. Fiz. Met. Metalloved. 1966, 22, 77. (45) McEwan, J. R. Diffusion of Ni 63 in Iron, Cobalt, Nickel and Two Iron-Nickel Alloys. Can. J. Chem. 1959, 37, 1629−1636. (46) Peterson, N. L. Self Diffusion in Pure Metals. J. Nucl. Mater. 1978, 69−7, 3−37. (47) Scholz, H. P. Messungen Der Absoluten Leerstellenkonzentration in Nickel Und Geordneten Intermetallischen Nickel-legierungen Mit Einem Differentialdilatometer, University of Gottingen, Goettingen, Germany, 2001. (48) Carling, K.; Wahnstrom, G.; Mattsson, T. R.; Mattsson, A. E.; Sanberg, N.; Grimvall, G. Vacancies in Metals: From First-Principles Calculations to Experimental Data. Phys. Rev. Lett. 2000, 85, 3862− 3865. (49) Mattsson, T. R.; Mattsson, A. E. Calculating the Vacancy Formation Energy in Metals; Pt, Pd and Mo. Phys. Rev. B 2002, 66, 214110. (50) Delczeg, L.; Delczeg-Czirjak, E. K.; Johansson, B.; Vitos, L. Assessing Common Density Functional Approximations for Ab Initio Description of Monovacancies in Metals. Phys. Rev. B 2009, 80, 205121. (51) Nastar, M.; Soisson, F. Atomistic Modeling of Phase Transformations: Point Defect Concentrations and the Time Scale Problems. Phys. Rev. B 2012, 86, 220102-1−220102-5.

22228


Simulation of Diffusion in FCC NiFe Binary Alloys Using Kinetic Monte

Recommend Documents