Quasi-atomistic modeling of the microstructure evolution in binary alloys and its application to the FeCr case

Acta Materialia 95 (2015) 324–334

Contents lists available at ScienceDirect

Acta Materialia journal homepage: www.elsevier.com/locate/actamat

Quasi-atomistic modeling of the microstructure evolution in binary alloys and its application to the FeCr case Ignacio Dopico a,⇑, Pedro Castrillo b, Ignacio Martin-Bragado a a b

Atomistic Materials Modeling Group, IMDEA Materials Institute, Eric Kandel 2, Getafe, 28906 Madrid, Spain UCAM, Universidad Católica de Murcia, Campus de los Jerónimos, Guadalupe, 30107 Murcia, Spain

a r t i c l e

i n f o

Article history: Received 1 April 2015 Revised 25 May 2015 Accepted 25 May 2015

Keywords: Atomistic modeling Phase transformation kinetics Monte Carlo techniques Fe–Cr Spinodal decomposition

a b s t r a c t In this work, we present a comprehensive quasi-atomistic Object Kinetic Monte Carlo (OKMC) model for diffusion-mediated decomposition in binary alloys, which is applied to the particular case of phase nucleation and spinodal decomposition in the iron–chromium system. The model describes atomistically the defects driving diffusion, while following the evolution of alloy concentrations by tracking the number of alloy atoms in the elements of an uniform mesh. Input parameters are defect diffusivities, tracer diffusivity ratios, and mixing energies, and they have been calibrated according to reported experiments and ab-initio calculations. Simulations based on this model are able to reproduce both phase nucleation in the metastable composition region and spontaneous phase decomposition and coarsening within the spinodal composition region. The convergence into the correct thermodynamics has been shown by comparing the simulation results to theoretical predictions, while the time evolution has been validated with experimental data for different alloy compositions. The simulation approach has proven to be suitable for extended annealing times and for domain sizes up to hundreds of nanometers. Ó 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

1. Introduction The understanding of the microstructure evolution of metals with time is critical to assess their structural properties as they age [1,2]. Of particular interest are the order–disorder transitions in binary alloys where, depending on temperature and composition, the alloy constituents might be homogeneously dispersed or form precipitates [3] that, in some cases, can even lead to full spinodal decomposition [4]. The kinetics of this evolution is extremely sensitive to the presence of radiation damage [5], which is especially important in the alloys used for nuclear applications [6]. Modeling and simulation of these systems is of extreme interest when experimental validation of material properties in real operation conditions is very costly or, even worse, undoable for still unexisting technologies [7]. However, predictive simulation is quite challenging, due to the involvement of atomistic mechanisms that have to be accounted at much higher length scales [8]. The FeCr system is a relevant example of what was mentioned above as, from a fundamental point of view, FeCr alloys show a rich physical phenomenology [9,10] and, from a practical perspective, ⇑ Corresponding author. Tel.: +34 91 549 34 22x1036; fax: +34 91 550 30 47. E-mail addresses: [email protected] (I. Dopico), [email protected] (P. Castrillo), [email protected] (I. Martin-Bragado). http://dx.doi.org/10.1016/j.actamat.2015.05.040 1359-6454/Ó 2015 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

FeCr-based steels exhibit excellent features for nuclear applications [11]. Under certain conditions, FeCr alloys decompose into an iron-rich phase (a) and a chromium-rich phase (a0 ) [12], with detrimental consequences for their structural properties [13]. This decomposition is driven by point-defect diffusion and it is an order–disorder transformation involving the a and a0 phases, both with the same body-centered cubic lattice (bcc) and very similar lattice parameters. The kinetics of a  a0 decomposition, usually very slow at moderate temperatures without irradiation, is accelerated by the presence of radiation-induced point defects, giving rise to Radiation Induce Segregation (RIS) [14]. In an effort to build a whole multiscale simulation platform, computational techniques going from ab-initio to continuum have been successfully used to test and understand material properties under different scenarios [7]. From these many techniques, kinetic Monte Carlo (KMC) methods, based either on event (EKMC) [2], atomistic (AKMC) [9,15] or object (OKMC) [16–18] algorithms, have been proposed [19]. Several studies have successfully shown the possibility to reproduce the Fe-Cr phase separation using AKMC [9,20,15] but, given the improved performance in computation time and memory consumption, quasi-atomistic models are desirable to overcome the limitations of AKMC models. Following this concept, one approach is the coarse-grained KMC simulation of interdiffusion proposed in Ref. [18], which demonstrates the

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

great potential of such models, keeping the links to atomistic mechanisms and comparing successfully to AKMC simulations, while avoiding a detailed treatment of the microscopic configuration of the alloy. A different approach is the defect-based quasi-atomistic paradigm, which allows to scale up simulation sizes and times and, hence, a favorable comparison to experimental results while still keeping an atomistic description of defects. Such an atomistic description of defects in KMC methods is critical for an accurate modeling of radiation-induced effects in iron [21] and FeCr steels, and needed for a comprehensive model including both damage evolution and Cr segregation and diffusion in an integrative way. In this work, a non-lattice OKMC model for decomposition between coherent phases in binary crystal alloys is presented. The model has been particularized to reproduce the thermodynamics and kinetics of the FeCr phase separation within a quasi-atomistic framework, and it has been implemented and tested on the MMonCa simulator tool [22]. The structure of the article is as follows: Section 2 describes the physics of the model and establishes a link between model parameters and measurable quantities obtained from experiments or ab-initio calculations. The relevant technical details related to the implementation are depicted in Section 3. Section 4 reports an illustrative set of simulation results for FeCr alloys, which are compared to experiments. Conclusions are wrapped up in Section 5.

2. Model Atom movements in crystalline metals are mainly driven by defect migration [23,24], with vacancies (V) being dominant near equilibrium conditions, but with self-interstitials (I) playing also a role under irradiation conditions [24]. Regarding alloys, atom diffusion may modify local concentrations and, consequently, intrinsic material properties. This diffusion is affected by the energetics of resulting microscopic configurations and it can give rise to either intermixing or phase decomposition. Thus, physically-based models for the structural evolution of crystal alloys must be built on a correct description of both point-defect diffusion and alloy formation energies. The model presented here inherits the ideas on defect diffusion developed in a previous work for semiconductor alloys [25], and incorporates the missing description of alloy formation energies in order to reproduce the phase diagram and decomposition behavior in binary crystal alloys and, in particular, in the interesting case of the FeCr system. This model is based on an atomistic treatment of defects, keeping track of their coordinates, while following the alloy evolution by counting the number of alloy atoms (but not the position) in every local region (cell). These cells are generated by splitting the whole simulation domain using a simple uniform mesh. The lattice is discarded everywhere in order to scale up the range of both simulation times and sizes, and the only position coordinates that are kept are those of defects. We do not distinguish between different self-interstitial or different vacancy types, lumping all of them into an effective self-interstitial or vacancy configuration [26]. This assumption has been previously done for iron [2], and it is extended here for iron-chromium, assuming that defect properties may depend on the alloy composition but not on the local atomic configuration. Fermi-level dependencies are negligible in metals and, therefore, no charge effects on diffusion [27] are considered in this work. Within the Transition State Theory [28], event rates are assigned to the different mechanisms of moving defects as r ¼ r0 expðE=kB TÞ, which are parametrized with the rate prefactor r 0 and the activation energy E, with kB T being the thermal energy. Once the complete list of rates is known, a regular KMC direct

325

method [29] is performed to follow the time evolution of the system. 2.1. Homogeneous alloys In homogeneous crystal materials, native point-defects K (K ¼ V or I) jump randomly with a thermally-activated migration frequency, mK , which is assigned to be:



mK ¼ m0K exp 

 EmK ; kB T

ð1Þ

with m0K and EmK being the migration prefactor and energy, respectively. The diffusivity, DK , is linked to the migration frequency by:

  k2 E mK ¼ D0K exp  mK ; 6 kB T

DK ¼

ð2Þ

where D0K is the diffusivity prefactor and k is the jump length, which is chosen to be on the order of the lattice parameter of the crystal. In equilibrium conditions, point-defect concentrations Ceq K are given by:

  EfK eq ; C eq ¼ C exp  K 0K kB T

ð3Þ

where C eq 0K is the equilibrium concentration prefactor, and EfK is the formation energy. During their movement, defects induce atom diffusion. In FeCr alloys, migration frequencies of atoms A (with A ¼ Fe or Cr) driven by K-defects, mKA , are related to defect migration frequencies as:

f

K

mK C K ¼ mKFe C Fe þ mKCr C Cr ;

ð4Þ

where C K is the defect concentration, C A is the atom concentration, and f

K

is the correlation factor for K-driven atom jumps, with

K

f 6 1. The K-components of tracer diffusivities, DKA , are given by:

k2 K m : 6 A

DKA ¼

ð5Þ

The different ability of K defects to move Cr or Fe atoms can be expressed by the ratio between their tracer diffusivity components, fK :

fK ¼

DKCr DKFe

ð6Þ

:

The input parameters of our defect-based model are taken to be eq K those related to DK C eq K (i.e.: C 0K ; EfK ; D0K ; EmK ) and to f , which is also assumed to follow an Arrhenius dependence, with prefactor fK0 and activation energy EKf . These parameters have been chosen to fit the experimental tracer diffusivities for pure materials (pure iron or pure chromium) near equilibrium [30,24]. Using Eqs. (4) and (6), DKFe and DKCr in Fe1x Crx near equilibrium conditions can be rewritten as: K

DKFe ¼

DK C eq K ; 1  x þ f x C tot

DKCr ¼

f fK DK C eq K ; 1  x þ fK x C tot

f

K

ð7aÞ

K

ð7bÞ

where x is the chromium fraction and C tot is the total atom concentration. In metals near equilibrium, vacancy mechanism dominates eq over self-interstitial mechanism (DV C eq V  DI C I ) and, therefore, self-interstitial contribution can be neglected. However, self-interstitials are considered in the formalism in order to be used in the future for the interesting case of irradiation conditions.

326

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

Table 1 Model parameters for vacancy-driven atom diffusion in pure iron and pure chromium materials near equilibrium conditions. C eq 0V Iron Chromium

EfV [eV]

D0V ½cm2 =s

EmV [eV]

fV0

½cm3  5e27 4e23

2.5 2

1e2 1.2e3

0.68 0.6

170 1

EfV

C tot ½cm3 

[eV] 0.32 0

8.77e22 8.77e22

Table 1 displays the input parameter values for the model, whereas experimental and modeled tracer diffusivities are shown in Fig. 1. In pure chromium, the diffusivities of Fe and Cr atoms, DVFe and DVCr , have been chosen to be equal due to the uncertainty of available experimental data and, thus, fV has been taken equal to the unit. The value of f

V

has been chosen to be the correlation V

factor for a bulk body-centered cubic lattice (f ¼ 0:727) [24]. For intermediate composition Fe1x Crx alloys, energies and entropies have been assumed to follow a linear dependence with x [25], implying an exponential dependence for prefactors.1 Defect to defect interaction occurs when a migrating defect jumps into the proximity of another defect. In particular, the interaction between two vacancies would give rise to a V 2 pair whereas the interaction between a vacancy and a self-interstitial would produce defect recombination. However, Ceq K is so low at usual temperatures that defect interactions are very improbable near equilibrium conditions and their effects are expected to be negligible. In contrast, these interactions would play a key role in the evolution of materials with a high concentration of radiation-induced defects. 2.2. Inhomogeneous alloys As previously mentioned, in our approach an inhomogeneous alloy is seen as divided in nanometer-sized homogeneous cells. Moreover, following a quasi-atomistic non-lattice scheme, we do not record the individual positions of lattice atoms (Fe or Cr) of the FeCr alloy in each cell i but just their integer numbers (nFei and nCri ). This allows to know the local concentration of Fe and Cr atoms keeping the discrete nature of matter, but without having to store in memory all the atomic positions. The local composition, characterized here by the relative Cr concentration, is given by:

nCrj Xj ¼ : nFej þ nCrj

ð8Þ

K For jumps inside a cell i, defect properties (mK i ; C eq K i , and fi ) depend

not only on its own local composition X i , but also on the chemical species of the atoms in adjacent cells that statistically would be first or second neighbors of atoms of cell i. In this way, the effective composition for defect properties, xi , is calculated as:

X xi ¼ wij X j ;

ð9Þ

j

where j denotes the cells adjacent to i including itself, wij is the weight of cell j on the effective composition of cell i, and X j is the unweighted local composition of cell j. Following Ref. [15], second neighbor contributions have been taken to weight half than first neighbor contributions. The expressions used to calculate the weights wij are explicitly derived in Appendix A. Defect jump frequencies crossing the borders between cell i and cell j, (mK i!j , from cell i to cell j, and mK j!i , from j to i) are taken to be related as [25]: 1

1x eq Linear dependence of entropies with x implies C V0 ðxÞ ¼ C eq C 0V ð1Þ x and 0V ð0Þ

mV0 ðxÞ ¼ m0V ð0Þ1x m0V ð1Þ x .

Fig. 1. (a) Self-diffusivity (solid line) and Cr diffusivity (dashed line) in iron used within this work, together with experimental reported values (dots) [30]. (b) Selfdiffusivity (dashed line) and Fe diffusivity (solid line) in chromium used within this work and experimental reported values (dots) [24].

eq mK i!j C eq K i ¼ mK j!i C K j :

ð10Þ

eq Thus, for DK i C eq K i > DK j C K j and using a step-like energy profile, these

frequencies are related to the bulk jump frequencies,

DK j C eq Kj DK i C eq Ki

mK i!j ¼

mK i ;

mK i , as [25]: ð11aÞ

mK j!i ¼ mK j :

ð11bÞ

A vacancy moving from cell i to cell j will displace an atom from j to i. In FeCr alloys, the ratio between vacancy-induced Cr and Fe fluxes from j to i (J VCrj!i and J VFej!i ) will depend on the relative concentrations in cell j and on the relative ability of a V to move a Cr or a Fe atom according to:

J VCrj!i J VFej!i

¼ fVj!i

Xj ; 1  Xj

ð12Þ

where:

fVj!i ¼

mVCrj!i : mVFej!i

ð13Þ

Consequently, the proportion of Cr atoms moved by vacancies from cell j to cell i, which we denote as X Vj!i , is given by:

X Vj!i ¼

J VCrj!i J VFej!i

þ

J VCrj!i

¼

fVj!i X j 1  X j þ fVj!i X j

:

ð14Þ

Here, fVj!i will depend not only on the bulk tracer ratios fV of Eq. (6), but also on the difference between the corresponding energy modifications. This is expressed as the ratio:

327

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

fVj!i

0 V 1 V qffiffiffiffiffiffiffiffiffiffi dEFei!j  dECri!j j!i j!i V V A; ¼ fi fj exp @ 2kB T

ð15Þ

V

i!j where dEAj!i stands for the total energy modification of the system

associated to a movement of an A atom from j to i due to a V-jump qffiffiffiffiffiffiffiffiffiffi from i to j, and fVi fVj is equal to the bulk value of fV for a composition x ¼

xi þxj . 2 V

V

The energy difference dEFei!j  dECri!j in Eq. (15) can be seen as j!i j!i the energy modification associated to an exchange of a Fe atom from cell j with a Cr atom from cell i. This energy difference is critical to model phase decomposition in FeCr alloys, as it controls diffusion-mediated segregation into Fe-rich (a) and a Cr-rich (a0 ) phases. Elastic energy contributions are assumed to be negligible in FeCr due to the small difference between the lattice parameters of a-iron and chromium (below 0.4%). This allows to write the total energy of an alloy cell with composition x as:

Etotal ¼ ðnL þ nI  nV ÞEbulk ðxÞ þ nI EfI ðxÞ þ nV EfV ðxÞ;

ð16Þ

where Ebulk is the internal bulk energy per atom in a defect-free Fe1x Crx crystal, and nL ; nI and nV are the number of lattice sites, self-interstitials and vacancies, respectively, in the cell. The internal bulk energy per atom, Ebulk can be written as:

Ebulk ¼ ð1  xÞEbulk ð0Þ þ xEbulk ð1Þ þ Emix ðxÞ;

ð17Þ

where Ebulk ð0Þ and Ebulk ð1Þ are the values for pure iron and pure chromium respectively, whereas Emix ðxÞ is the so called mixing energy and it accounts for the deviation from linear behavior. From Eqs. (16) and (17), and assuming first order approximation for nL  1, the difference between the energy modifications is evaluated to be (see Appendix B): V

V

dEFei!j  dECri!j ¼ j!i j!i

X ðwlj  wli ÞE0mix ðxl Þ;

ð18Þ

l

with E0mix ðxÞ being the first derivative of Emix respect to the effective composition x. A formally similar derivation can be obtained for self-interstitial assisted atom movements, but with defects and atoms moving in the same direction. This energy modification is also analyzed in Appendix B, and the result for atoms moved by I defects from j to i is analogous to the one of Eq. (18). In both cases (V an I), Eq. (15) results in:

fKj!i ¼

P  qffiffiffiffiffiffiffiffiffiffi 0 l ðwlj  wli ÞEmix ðxl Þ : fKi fKj exp 2kB T

ð19Þ

Emix has been postulated to depend also on temperature, T, as [15]:

Emix ðx; TÞ ¼



1

T

 Emix ðx; 0Þ;

H

Fig. 2. Mixing energy curves at 0 K calculated by different authors using ab-initio techniques [31–34], and the curve selected for this work.

ð20Þ

with Emix ðx; 0Þ being the mixing energy at 0 K and H being the unmixing critical temperature, which is related to vibrational entropy [15]. In that case, E0mix ðxÞ in Eq. (19) denotes the partial derivative of Emix ðx; TÞ with respect to x for the given T. Emix ðx; 0Þ has been calculated by different authors using ab-initio techniques [31–34], as shown in Fig. 2.

For an ideal binary alloy the mixing free energy at temperature T; Gmix , is given by [35]:

Gmix ¼ Emix ðx; TÞ þ kB T½x ln x þ ð1  xÞ lnð1  xÞ;

ð21Þ

where, in our case, the mixing energy Emix is equal to the mixing enthalpy since constant null pressure has been assumed. The last term of Eq. (21) corresponds to the contribution of the configurational entropy for an ideal binary alloy [15]. In the atomistic model described in Section 2.2, this term arises from atom statistics (see Eq. (12)) and, then, in contrast to Emix , configurational entropy has not to be included into the relative jump probabilities of individual atoms (Eqs. (15) and (19)). The mixing free energy can be related to the spinodal decomposition range and to the miscibility gap of the alloy. Therefore, Eq. (21) can be used as a link between Emix and the expected phase separation behavior. In particular, the alloy becomes unstable for G00mix < 0 [35], where G00mix is the second partial derivative of Gmix respect to x. Thus, spontaneous spinodal decomposition occurs for T < T sp , with:

 T sp ðxÞ ¼

1

H



kB xð1  xÞE00mix ðx; T ¼ 0Þ

1 ;

ð22Þ

where the temperature dependence of Emix ðx; TÞ formulated in Eq. (20) has been assumed. Phase coexistence between the equilibrium composition of the iron-rich phase, xa , and the equilibrium composition of the chromium-rich phase, xa0 , at a given T must satisfy that the first partial derivative of Gmix has to be equal for the two phases (G0mix ðxa Þ ¼ G0mix ðxa0 Þ) and the total free mixing energy of the system has to be minimized, being equivalent to minimize the average Gmix . For an average composition  x decomposing into xa and xa0 ; Gmix will be:

Gmix ¼

x  xa xa0  x Gmix ðxa Þ þ Gmix ðxa0 Þ: xa0  xa xa 0  x a

ð23Þ

Applying the above mentioned conditions to Eq. (23), it yields: 2.3. Comparison to thermodynamics The model described above can be readily used to simulate the concentration evolution of the alloy system. Even so, a comparison to thermodynamic predictions for a given Emix ðx; TÞ dependence is useful both as a guide for parameter calibration and as a test for model validation.

G0mix ðxa Þ ¼ G0mix ðxa0 Þ ¼

Gmix ðxa0 Þ  Gmix ðxa Þ ; xa0  xa

ð24Þ

which is equivalent to the common tangent method [35]. The points ðx; TÞ satisfying Eq. (24) constitute the curve that separates the stable region from the metastable region of the alloy and define the miscibility gap of the a þ a0 solid solution.

328

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

If Emix ðx; 0Þ and H are known, it is possible to infer the spinodal region, using Eq. (22), and the miscibility gap, using Eq. (24). Thus, a mixing energy curve Emix ðx; 0Þ (parametrized in Appendix C) and a value of H ¼ 1400 K have been selected in this work to be consistent to both the ab-initio calculated mixing energies [31–34] (see Fig. 2) and to the available information on phase separation in FeCr alloys [31–34,39]. In particular, a recent work reporting segregation of a Fe0:85 Cr0:15 sample at 475 °C [39], has been taken into account. The thermodynamic miscibility gap calculated from our parametrization of Emix is displayed in Fig. 3 together with reported reference results. Notice that the contribution of the a  a0 interface has been neglected within the thermodynamic considerations of this subsection. Hence, the miscibility gap concentrations calculated with Eq. (24) would correspond to separated bulk-like phases with no interface effects. In contrast, the simulation method of Section 2.2 accounts for the contribution of all cells, including those at the border between a and a0 domains. The smaller the phase domains, the higher the relative contribution of interface regions and, then, the lower the expected convergence between thermodynamic predictions and atomistic simulations.

Fig. 3. Miscibility gap for the a þ a0 FeCr solid solution. Dashed lines stand for some proposed miscibility gaps [10,36–38], whereas the solid line is the miscibility gap inferred from thermodynamics for the mixing energy proposed in this work.

K

3. Implementation The model has been implemented within a quasi-atomistic Object kinetic Monte Carlo framework into the MMonCa [22] Open Source Simulator.2 The non-lattice OKMC method divides the volume into a homogeneous mesh. Each point-defect is represented by its type (vacancy or self-interstitial) and position. Point-defects migrate performing jumps of length k in one of the six orthogonal directions of the space, with the frequency described in Section 2.1. k has been taken to be 0.3 nm, approximately equal to the lattice parameter of FeCr alloys. After a jump, interacting neighbor defects, if any, are searched within an radius equal to k, following the same scheme described in Ref. [22].3 Point-defect migration depends on the effective cell composition, xi , which is computed from Eq. (9). When a point defect migrates from cell i to cell j it faces an energy barrier due to the difeq ference between effective cell compositions [25]. If DK j C eq K j < DK i C K i , the jump is accepted with a probability PK i!j , where:

PK i!j ¼

DK j C eq Kj DK i C eq Ki

;

is always accepted (see Eq. (11b)). In local equilibrium conditions (i.e., when the number of K-defect jumps from i to j and from j to i are nearly the same), alloy interdiffusion can be atomistically described by a net atom exchange between cells induced by point-defect migration [25,18]. In FeCr alloys, the possible atom exchanges are: Fei $ Fej ; Fei $ Cr j ; Cr i $ Cr j and Cr i $ Fej . Exchanges of atoms of the same type leave the system unchanged and can be discarded in the simulations. Non-trivial exchanges have a probability of occurrence given by [25]:

2

ð26Þ

http://materials.imdea.org/MMonCa. Nevertheless, we have checked that defect interactions are very rare in the simulations presented in this work and they can be ignored with no detrimental consequences on the simulation results. 3

tion domain. The simulator time step is evaluated as the inverse of the total frequency of the executed events [22]. The other option is to consider periodic boundary conditions in all three axes and introduce one vacancy into the simulation domain. In this case, the real time of the simulation system can be approximated by generalization of the expression reported in Ref. [15]:

ð25Þ

which is the probabilistic counterpart of Eq. (11a). Conversely, if eq DK j C eq K j P DK i C K i , there is no additional energy barrier and the jump

  PKCri $Fej ¼ g K X Ki!j 1  X Kj!i ;

with g K ¼ f2Lk ; L being the cell length, and X Ki!j and X Kj!i being a generalization of Eq. (14) for K ¼ V; I. If an exchange is selected to occur, the atom counters nFe and nCr of the involved cells will be modified 1, keeping nFe þ nCr constant. This algorithm has been previously proven to be very efficient and suitable for SiGe alloys [25]. Periodic boundary conditions have been set in two of the three axes. For the other axes, two different options have been used: free surface plus mirror conditions or periodic conditions. In the first case, a free surface is set at the front side, and mirror conditions are imposed at the back side. The balance between point-defect emission and capture at the surface, together with the balance between cells described in Eq. (11), establishes near equilibrium point-defect concentrations (C K i ’ C eq K i ) in every cell i of the simula-

tP

tMC eq ; vol iCVi i

ð27Þ

where t MC refers to the Monte Carlo simulator time based in the executed migration events, and voli is the volume of each cell. As C eq V i is composition dependent and the compositions of the simulation cells are changing, time scaling has to be dynamically evaluated during the simulation. We have checked that the two approaches provide similar results and we have mainly used the second one. 4. Simulation results Fig. 4(a) validates our model with thermodynamic results by presenting the simulated miscibility gap for the FeCr system (color histogram) against the theoretical calculated miscibility gap (solid black line), previously displayed in Fig. 3. The figure has been built from a set of 40 simulations from 300 °C to 700 °C in 10 °C intervals. Each simulation uses a 16  16  16 nm3 domain with periodic conditions and begins with an in-depth Fe0:45 Cr0:55 / Fe0:55 Cr0:45 sinusoidal superlattice concentration profile, set to trigger and easily identify the alloy decomposition into a0 /a regions, as shown in Fig. 4(b) (squares). Once the simulation reach

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

Fig. 4. (a) Set of simulated composition histograms of the Cr fraction in FeCr alloys annealed at temperatures between 300 °C and 700 °C in 10 °C intervals (color histogram), compared to the miscibility gap Cr fractions predicted by thermodynamics for separated bulk phases (solid line). The spinodal decomposition onset is also shown (dashed line). Histograms have been computed by running simulations for alloys with an initial sinusoidal in-depth composition profile, with maximum and minimum amplitudes at Fe0:45 Cr0:55 and Fe0:55 Cr0:45 respectively, as displayed in (b) for 400 °C (squares). The final in-depth composition distribution, showing ironrich (a-phase), chromium-rich (a0 -phase), and interface regions, is also presented (color histogram). Finally, (c) displays the joint histogram of relative occurrence of chromium fraction composition at 400 °C. Thermodynamic miscibility gap values at 400 °C are marked by horizontal dashed lines in (b) and by vertical solid lines in (c), exhibiting a good agreement with the maxima of the simulated bimodal distribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

equilibrium, two nearly homogeneous Cr-rich and Fe-rich regions, separated by a heterogeneous interface, are formed. This can be seen in the in-depth color histogram displayed in Fig. 4(b) for a particular temperature (400 °C). Joint histograms are then built for each temperature, with the relative occurrence of chromium fractions as shown in Fig. 4(c), also for 400 °C, and represented in color scale at the corresponding temperature in Fig. 4(a). The theoretical miscibility gap chromium concentrations for 400 °C are indicated to allow comparison in Fig. 4(b)) and (c) as horizontal and vertical lines, respectively. The agreement of the phase diagram and bimodal histograms with theoretical values is remarkable and it validates our approach and implementation. Additionally, simulation results in Fig. 4 provide insights into the asymmetric concentration distribution of a and a0 phases, coming from the asymmetry of the two mixing free-energy valleys but also from the contribution of cells at the interface between a and a0 regions, not accounted by the thermodynamic calculation for bulk phases. For temperatures above 650 °C, no a  a0 phase separation

329

is expected, as it is confirmed by the presence of an undivided band at the top of Fig. 4(a), corresponding to a monomodal concentration distribution. Fig. 5 shows one of the main characteristics of spontaneous spinodal decomposition: the formation of an interwoven ‘‘vein-like’’ network of a and a0 phases starting from a homogeneous composition alloy [40]. This morphology has been produced by running a long term aging simulation at 470 °C of a Fe0:5 Cr0:5 alloy and displaying a section of 100  100 nm2 at four different times, with time evolution increasing left to right and top to bottom. Contrary to the simulations shown in Fig. 4, where the a and a0 phases had their seeds in the valleys and hills of the initial superlattice, in Fig. 5 the initial evolution towards a or a0 phases is related to the random concentration fluctuations of the initial sample, which evolves by coarsening and shape rounding, giving rise to the characteristic spinodal decomposition morphology. A phase diagram could be also obtained in this way, but due to the higher presence of interfaces and to the complex geometry generated, it is less efficient and more ambiguous than using a preset quasi-one dimensional structure, as in Fig. 4. It is worthy to mention that, in order to obtain the ‘‘vein-like’’ morphology of Fig. 5, it is necessary to include the contribution of neighboring cells on the energetics of a given cell. Within our approach, this is given by the weights wij of cells j on the effective composition xi (and, thus, on Emix ðxi Þ) of cell i (see Eq. (9)). If these intercell contributions were not taken into account, the concentration of each cell would still evolve spontaneously towards the miscibility gap values, but forming an artificial ‘‘chess-like’’ topology, and no connected a and a0 phase domains would be formed. In contrast, when these contributions are included, the system evolves toward configurations with connected bulk-like a and a0 domains, separated by a reduced amount of cells acting as an interface between phases. These cells have intermediate effective composition (coming from the weighted average of cells with high and low Cr content, see Eq. (9)) and, therefore, high mixing energy (see Fig. 2), resulting in an extra energy related to interfaces. This energy penalty is minimized through morphology evolution. In this way, the interface contribution appears in a simple but natural way in the model, without the need of an additional implementation. Fig. 6 shows the impact of mesh spacing L on the accuracy and morphology of simulated results. A set of simulations with different mesh spacing, L, and a fixed volume of 60  60  60 nm3 were performed to rule out artifacts related to mesh spacing. An initial homogeneous Fe0:5 Cr0:5 alloy is annealed a fixed time with different L values of 0.75 nm, 1 nm, 1.5 nm, and 2 nm (left axis). For each L, the intercell weights wij were recalculated according to Eqs. (A.1)–(A.3). In order to display the results with the same spatial resolution, the simulations with smaller cells were merged into larger cells. In particular, (2  2  2) and (3  3  3) 0.75 nm cells were merged into 1.5 and 2.25 nm cells respectively, and (2  2  2) 1 nm cells into 2 nm cells. The diagonal values represent the original simulations, the off-diagonal, the merged ones. It can be seen that, within this spacing range, the effect of changing mesh spacing in the simulation is quite similar to a modification of the output resolution, while still producing the same physics. A lower limit for L in the current algorithm is L ¼ 2a  0:6 nm, with a being the lattice parameter (see Appendix A), whereas the larger bound is imposed by the lowest size of the features to be observed. Larger mesh spacing offers higher computational efficiency, with the drawback of a lower resolution. As a compromise, a typical spacing of 1 nm has be chosen in our simulations. A similar analysis with analogous conclusions has been performed for simulations near the nucleation regime. In Fig. 7, the experimental results reported in Ref. [40] are compared to our model predictions. A Fe0:68 Cr0:32 sample was aged for

330

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

Fig. 5. Simulation of the evolution of an initially homogeneous Fe0:5 Cr0:5 alloy under a long term aging at 470 °C. Sections at different times are shown, with time increasing left to right and top to bottom.

and even the initial distribution of Cr is different, full agreement cannot be expected, but a good match in the evolution of the concentration and the maximum amplitudes is observed, showing agreement in the time evolution of the spinodal decomposition of the sample. Fig. 8(a) shows a 50  9  9 nm3 section of a Fe0:8 Cr0:2 sample aged at 500 °C for a maximum of 1067 h. Experimental measurements of Energy Compensated Tomographic Atom Probe (EC-TAP) from Ref. [12] (left) are compared to the simulated nanostructure evolution obtained with OKMC (right). Nucleation, growth, and coarsening of the precipitates can be tracked in the sequence, and qualitative agreement with our simulations is observed for both space and time evolutions, taking into account that in the EC-TAP images only the precipitates with a chromium fraction > 0:4 were shown. A quantitative comparison between experiment [12] and OKMC simulations is presented in Fig. 8(b)–(d), including also lattice AKMC results of Ref. [15]. A 60  60  60 nm3 domain with cell spacing of 1 nm was used in our OKMC simulations. An initial binomial Cr distribution composition centered in x ¼ 0:2 was assumed to describe a random Fe0:8 Cr0:2 alloy. In order to distinguish precipitate cells from matrix cells and interface cells in the analysis of annealing simulations, a cell has been considered to be part of an a0 precipitate if it is above 0.4 chromium fraction and at least five of its six next-neighbor cells are above a 0.4 chromium fraction. Analogously, a cell is considered to be part of the a matrix if it is below a 0.4 chromium fraction and at least five of its six next-neighbors are below a 0.4 chromium fraction. Otherwise the cell is considered to be part of an interface. Fig. 8(b) displays the mean radius of a0 precipitates. The radius R of each precipitate in the simulation is estimated assuming a spherical shape and, thus, taking

R¼

Fig. 6. Simulations of an initially homogeneous Fe0:5 Cr0:5 alloy under a long term aging at 500 °C using different cell spacing (vertical axis). Diagonal views correspond to as-simulated results, while off-diagonal views are produced by merging cells into larger ones (see text). In this way, each row shows same simulations displayed with different spatial resolution (horizontal axis), while columns compare simulations performed with different cell spacing but displayed with similar resolution.

669 h at 470 °C revealing fully spinodal decomposition (‘‘vein-like’’ morphology). In-depth composition profiles at 0, 50, 193 and 669 h are shown in the Figure. A domain size of 100  20  20 nm3 has been used in the simulations. Since the process is fully stochastic,

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 3ncell v ol ; 4p

ð28Þ

where ncell is the number of cells belonging to the precipitate, and vol is the cell volume. Since the cell size of the simulated volume was chosen to be 1 nm, precipitates with a radius lower than about 1 nm (reached in average in 150 h) are below the resolution of the simulator, represented in the figure by shaded areas, and should not be considered. Fig. 8(c) shows the number density evolution of a0 precipitates with time, and Fig. 8(d) represents the Cr fraction of both the precipitates and the matrix with time. Shaded areas represent, again, regions where the mean radius is below the resolution of our simulations and, hence, OKMC results are not reliable for the chosen cell size. Otherwise the agreement to experimental measurements is very good and comparable to AKMC simulations. Ref. [12] also reported the experimentally estimated fraction of percolated precipitates. A comparison with simulation results can be seen in Fig. 9. The top panel displays a sequence of simulation images where percolation can be traced back from early stages, and the subsequent growth produces a spherical precipitate. In our model, the driving force of this morphological evolution is the minimization of the mixing energy Emix , which is related to the reduction of the a–a0 interface. Finally, Fig. 10 compares precipitate nucleation and growth of a Fe0:85 Cr0:15 alloy annealed at 475 °C as observed by bright field transmission electron microscopy [39], to the results of OKMC simulations. Images correspond to annealing times of 500, 2000 and 5000 h, and a domain of 400  400  20 nm3 has been used in the simulations. A good agreement, both in morphology and kinetics, is observed. It should be noticed that, given the volume of material involved, simulations like the previous one are challenging for current AKMC models [15], but can be done in regular desktop computers within

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

331

Fig. 7. In-depth chromium concentration profile of a Fe0:68 Cr0:32 sample as quenched (top left) and after 50 h (top right), 193 h (bottom left) and 669 h (bottom right) aging at 470 °C. The experimental data [40] are plotted with solid black lines whereas OKMC simulation results are plotted with solid red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Compilation of OKMC simulation results for a Fe0:8 Cr0:2 sample aged at 50 °C, compared to experimental data [12] and lattice AKMC simulations [15]. (a) Evolution of chromium-rich a0 precipitates in a (50  9) nm2 section of the sample. Left column images correspond to Atom Probe Tomography [12]5 and right column images are the simulation results. (b) Mean radius evolution of the a0 precipitates. (c) Number density evolution of the a0 precipitates, together with a slope trend. (d) Average chromium concentration evolution of a0 precipitates and a matrix. Symbols are experimental data [12], dashed lines are simulation results of a lattice AKMC model [15], and solid lines are the results of the presented OKMC model. Shaded areas in (b)–(d) correspond to conditions in which the results of our OKMC approach are not reliable due to spacial resolution limitation.

our OKMC approach. Using our model, the computation time for the simulation in Fig. 10 was 19 days on a single core processor. Computation time for the simulation of Fig. 8 for a volume of 100  30  30 nm3 was 19 h on a single core, and only 30 min for the simulation of Fig. 7. To further quantify the performance improvement with respect to AKMC models, we have run a 9  9  50 nm3 simulation, comparable to the 32  32  176 (lattice units)3 simulation of Ref. [15]. The wall clock time was 37 min with our model, whereas it was 11.8 days with the AKMC approach4 for the same 1067 h annealing. Therefore, in this example, the performance improvement achieved with our method is about a 4

Martinez, E. Personal Communication. 05/05/2015.

factor 460. The trade off incurred for such efficiency is the limited capability of our OKMC model to track very small precipitates, with radius typically bellow 1 nm. Regarding other quasi-atomistic models, previous coarse grained KMC models are, so far, not reporting comparisons with meso-scale experiments, and their authors recognize limitations for the simulation of nucleation [18]. In contrast, our work shows nucleation and phase decomposition, and how to correlate them with the input parameters that define the system thermodynamics,

5 Reprinted from Journal of Nuclear Materials, vol. 384, S. Novy, P. Pareige, C. Pareige, Atomic scale analysis and phase separation understanding in a thermally aged Fe–20 at.%Cr alloy, pp. 96–102, Copyright Ó2008, with permission from Elsevier.

332

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

and diffusion of individual defects, which makes it promising for a future OKMC description of Radiation Induced Segregation [21]. Further work will be carried to extend the model beyond the local equilibrium assumption used here. 5. Conclusions

Fig. 9. Evolution of the fraction of percolated chromium precipitates in Fe0:8 Cr0:2 at 500 °C. Symbols: experimental data [12]. Line: simulation results. Top panel: time sequence of a percolation process as seen in the simulations.

In this work, an efficient atomistic model for diffusion and phase separation in binary alloys has been reported. The model uses a quasi-atomistic approach that, while still keeping track of the atomistic defects in the material, follows the evolution of the alloy composition using mesh elements. Input parameters are related to point-defect diffusivities, tracer-atom diffusivities, and mixing energy of the alloy. The model has been applied to the FeCr system. It has been implemented into MMonCa [22], a modular Open Source Kinetic Monte Carlo simulator. The model has been validated against theoretical predictions for bulk-like phases, providing the correct thermodynamics. Spatial resolution of the method, related to mesh spacing, has been analyzed. Simulations have been compared to microscopy experiments for FeCr alloys with different Cr concentrations and annealing times. A good agreement between simulations and experiments is obtained in the morphology and kinetics of the system, both in the spinodal decomposition regime and in the precipitate nucleation and growth regime. A computational efficiency improvement of more than two orders of magnitude has been obtained over state-of-the-art AKMC methods. These results demonstrate that significant domain sizes and annealing times can be accurately simulated using the present OKMC model, without the practical limitation of current AKMC models. This capability, added to the atomistic description of defects, makes the model promising to be extended to the simulation of radiation-induced segregation, while providing a formalism that can be easily applied to other binary alloys. Acknowledgments I.D. and I.M.-B. acknowledge funding of the project MASTIC (PCIG09-GA-2011-293783) by the Marie Curie Actions Grant FP7-PEOPLE-2011-CIG program, I.M.-B. also acknowledges funding by the Spanish Ministry of Economy and Competitiveness (MINECO) under Ramón y Cajal fellowship RYC-2012-10639, and partial funding from the FerroGENESYS project (MAT201347460-C5-3-P) funded by the same Ministry. P.C. acknowledges funding by UCAM under PMAFI/04/14 project. This work contributes to the Joint Programme of Nuclear Materials (EERA JPNM) of the European Energy Research Alliance. We thank E. Martinez for providing us the data on AKMC simulation times. Appendix A. Contribution of adjacent cells to the effective composition Let us consider a body-centered cubic (bcc) lattice divided in cubic regions (cells), each cell formed by N 3 unit-cells and, thus,

Fig. 10. Bright field transmission electron micrographs [39] (left) and simulated sections (right) of Fe0:85 Cr0:15 alloy annealed at 475 °C at different aging times: (a) before aging, and after (b) 500, (c) 2000, and (d) 5000 h.6

even for FeCr alloys, which exhibit a phase diagram less symmetric than the FeCu case of Ref. [18]. Moreover, our model establishes a correlation between alloy segregation evolution and the position 6 Reprinted from Journal of Nuclear Materials, vol. 455, D. Chen, A. Kimura, W. Han, Correlation of Fe/Cr phase decomposition process and age-hardening in Fe–15 at.%Cr ferritic alloys, pp. 436–439, Copyright Ó2014, with permission from Elsevier.

containing 2N 3 atoms each, with N being a natural number. In a bcc lattice, every atom has 8 first-neighbors and 6 secondneighbors. Hence, there are

1 2

 8  2N 3 ¼ 8N 3 first-neighbor

(bond-like) interactions and 12  6  2N 3 ¼ 6N 3 second-neighbor interactions per cell. It can be assumed that the properties of a cell depend mainly on the chemical nature of the atoms of their firstand second-neighbor interactions, and that the two atoms of the interaction contribute with the same weight [15]. Although for large enough cells (for example, N > 3) most of the interactions occur with atoms of the same cell, some of them take place with atoms of adjacent cells.

333

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

In the case of FeCr alloys and assuming a random atom distribution within each cell, the average fraction of Cr atoms in an interaction between two lattice sites of cells i and j would be ðX i þ X j Þ=2, with X i ; X j being the Cr fraction in cells i and j. Therefore, the properties of a cell i and, thus, its ‘‘effective composition’’, xi , would depend on the Cr fraction of adjacent cells, X j , as P P xi ¼ j wij X j (Eq. (9)), with j wij ¼ 1. The weights wij can be derived from the number of interactions between atoms of i and j cells. 2

It is possible to show that each cell has 6ð2N  1Þ first-neighbor interactions with adjacent cells sharing a cube face (‘‘face cells’’), 6ð2N  1Þ with cells sharing one cube edge (‘‘edge cells’’), and 2 with cells sharing one cube corner (‘‘corner cells’’). Consequently, first-neighbor contributions of face, edge, and cor-

Table 2 Initial and final Cr fraction of cells i and j (X iini ; X jini ; X ifin ; X jfin ) for Fe or Cr atom movements from j to i, associated to vacancy movement from i to j. The modifications dX i and dX j , expressed in terms of X iini and X jini , are also shown. X iini

X jini

X ifin

X jfin

dX i

dX j

V i!j Fej!i

nCri

ini

nCrj

ini

nCri

nCrj

nL

nL

X i ini nL

Xj ini nL 1

V i!j Crj!i

nCri

ini

nCrj

ini

nCri

1X i ini nL

1þX j ini nL 1

nL 1

nL 1

nL

ini

ini

ini

nL 1 þ1

nL

nCrj

ini

1

nL 1

ner cells are 3ð2N  1Þ2 ; 3ð2N  1Þ, and 1, respectively, over a total

cell. This can be computed using Eq. (9) and Appendix A. According to that, the modification of the effective composition of a cell l; dxl , can be calculated as:

of 8N 3 . Regarding second-neighbor interactions (and assuming

dxl ¼

2

3

N P 2), face cells have a contribution of 6N over a total of 6N , whereas there are no second-neighbors interactions with edge or corner cells. We denote by p the weight of a second-neighbor interaction relative to a first-neighbor interaction, with 0 6 p 6 1. Therefore, if cell j is one of the 6 face cells of i, the value of wij is taken to be:

wij ¼ wface ¼

1 3ð2N  1Þ2 þ 6pN 2 r 4 þ 2p  4r þ r 2 ¼ ; 6 4 4 þ 3p 8N3 þ 6pN3

ðA:1Þ

where r ¼ N1 ¼ aL, with a being the lattice parameter and L the cube size. For j corresponding to an edge cell, wij is taken to be:

wij ¼ wedge ¼

1 3ð2N  1Þ r2 2  r ¼ ; 3 3 12 8N þ 6pN 8 4 þ 3p

ðA:2Þ

which is actually the average weight of the 12 edge cells. Analogously, for corner cells, wij is taken to be the average of the 8 corner cells:

wij ¼ wcorner

1 1 r3 1 ¼ ¼ : 3 3 8 8N þ 6pN 16 4 þ 3p

ðA:3Þ

The weight of a cell on itself can be computed using that P wii ¼ 1  j–i wij . Expressions (A.1)–(A.3) have been derived for natural values of N, with N P 2. In our non-lattice approach, these expressions are extrapolated to arbitrary cell sizes L (with L P 2a), and r is calculated as aL with a ¼ 0:287 nm (equal to the lattice parameter of bcc iron). A value of p ¼ 0:5 has been used, following Ref. [15]. For a typical case of L ¼ 1 nm, the weight values obtained 2

X wlm dX m ¼ wli dX i þ wlj dX j ;

where it has been used that dX m ¼ 0 for m – i; j. Therefore, according to Table 2, the values of dxl for Fe or Cr movements   V V dxl i!j and dxl i!j are: Fej!i Crj!i

V i!j

dxl

Fej!i V i!j

dxl

Cr j!i

¼ wli

¼ wli

X iini X jini þ wlj ; nL nL  1

1  X iini 1  X jini  wlj : nL nL  1

The actual Cr fraction in a region (cell) i of a binary FeCr alloy is given by:

Xi ¼

nCri nCri ¼ ; nFei þ nCri nL  nV i þ nIi

ðB:1Þ

with nFei ; nCri ; nV i ; nIi ; nL being the number of Fe atoms, Cr atoms, vacancies, self-interstitials, and lattice sites, respectively, in the cell. The same size (implying the same nL ) is assumed for all cells. Let us consider the movement of an atom (either Fe or Cr) from cell j to cell i induced by a vacancy jumping from i to j in a FeCr alloy. The initial and final values of X i and X j are related to the initial nCr values according to Table 2. The modifications of X i and X j (denoted as dX i and dX j ) expressed in terms of the initial values are also indicated in the table. The effective composition of a cell accounts for the average Cr fraction of the atoms that are neighbors of the atoms of the given

ðB:3bÞ

difference can be evaluated from Eqs. (16) and (17) as the difference of the total energy of the final states:

    V i!j V i!j V V i!j  E dEFei!j  dE ¼ E x þ dx x þ dx j j j j fV fV Cr ini ini j!i j!i Fej!i Cr j!i    

X  V i!j V i!j nL  dlj Emix xlini þ dxl  Emix xlini þ dxl : þ Fej!i Cr j!i l ðB:4Þ Using Eqs. (B.3a) and (B.3b), and taking first order approximation in Emix ðxÞ and EfV ðxÞ for nL  1, Eq. (B.4) becomes: V

V

dEFei!j  dECri!j  j!i j!i

X

 wlj  wli E0mix ðxl Þ:

ðB:5Þ

l

with Eqs. (A.1)–(A.3) are: wface ¼ 4:8  10 ; wedge ¼ 2:7  10 ;

Appendix B. Difference between total energy modifications

ðB:3aÞ

We analyze now the difference between the total energy modifications associated to the above-considered movements  V  V dEFei!j . As the initial state is the same both cases, this  dECri!j j!i j!i

3

wcorner ¼ 2  104 , and wii ¼ 0:678.

ðB:2Þ

m

E0mix ðxl Þ denotes the derivative of Emix ðxÞ with respect to the composition x evaluated at the effective composition of cell l or, if the temperature dependence is considered (as in Eq. (21)), the partial derivative of Emix ðx; TÞ with respect to x, evaluated at xl for the given temperature T. For self-interstitial driven diffusion, an atom movement from cell j to cell i is induced by a self-interstitial going from j to i. Analogously to the V-driven case, the effective composition modifications are found to be:

dxl

dxl

V i!j Fej!i V i!j Cr j!i

¼ wli

¼ wli

X iini Xj þ wlj ini ; nL þ 1 nL

1  X iini 1  X jini  wlj : nL þ 1 nL

ðB:6aÞ

ðB:6bÞ

and the difference between energy modifications results to be: I

I

dEFej!ij!i  dECrj!ij!i 

X

 wlj  wli E0mix ðxl Þ;

l

which is analogous to Eq. (B.5).

ðB:7Þ

334

I. Dopico et al. / Acta Materialia 95 (2015) 324–334

Table 3 Coefficient values for Eq. (C.1). x0

x1

a3

a2

a1

a0

b2

b1

b0

0.09

0.20

210

170

56

7.9

0.15

0.535

0.05

Appendix C. Composition dependence of mixing energy In this work, the mixing energy of Fe1x Crx alloys at T ¼ 0 K has been taken to be a two-region polynomial function of the Cr fraction x as:

( Emix ðx; 0Þ ¼

xðx  x0 Þða3 x3 þ a2 x2 þ a1 x þ a0 Þ if x < x1 ð1  xÞðb2 x2 þ b1 x þ b0 Þ

if x P x1

ðC:1Þ

In order to get T sp ðxÞ continuum and derivable (see Eq. (22)), Emix ðx; 0Þ has been imposed to be continuum up to the third derivative. The values of the coefficients are listed in Table 3. References [1] Committee A.I.H. ASM Metals Handbook, Metallography and Microstructures, vol. 9, ASM International, 1985. [2] C.C. Fu, J.D. Torre, F. Willaime, J.L. Bocquet, A. Barbu, Multiscale modelling of defect kinetics in irradiated iron, Nat. Mater. 4 (1) (2005) 68–74. [3] S.M. Allen, J.W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall. 27 (6) (1979) 1085–1095. [4] J.W. Cahn, On spinodal decomposition, Acta Metall. 9 (9) (1961) 795–801. [5] M. Bachhav, G.R. Odette, E.A. Marquis, a0 precipitation in neutron-irradiated FeG (Cr alloys, Scr. Mater. 74 (2014) 48–51. [6] R.L. Klueh, D.R. Harries, High-Chromium Ferritic and Martensitic Steels for Nuclear Applications, ASTM, West Conshohocken, PA, 2001. [7] L. Malerba, A. Caro, J. Wallenius, Multiscale modelling of radiation damage and phase transformations: The challenge of FeCr alloys, J. Nucl. Mater. 382 (2) (2008) 112–125. [8] T. Suzudo, Y. Nagai, D. Schwen, A. Caro, Hardening in thermally-aged Fe–Cr binary alloys: Statistical parameters of atomistic configuration, Acta Mater. 89 (2015) 116–122. [9] G. Bonny, D. Terentyev, L. Malerba, D.V. Neck, Early stages of a–a0 phase separation in Fe–Cr alloys: an atomistic study, Phys. Rev. B 79 (10) (2009) 104207. [10] W. Xiong, P. Hedström, M. Selleby, J. Odqvist, M. Thuvander, Q. Chen, An improved thermodynamic modeling of the Fe–Cr system down to zero kelvin coupled with key experiments, Calphad 35 (3) (2011) 355–366. [11] F. Garner, M. Toloczko, B. Sencer, Comparison of swelling and irradiation creep behavior of fcc-austenitic and bcc-ferritic/martensitic alloys at high neutron exposure, J. Nucl. Mater. 276 (1) (2000) 123–142. [12] S. Novy, P. Pareige, C. Pareige, Atomic scale analysis and phase separation understanding in a thermally aged FeG (20 at.% Cr alloy, J. Nucl. Mater. 384 (2) (2009) 96–102. [13] O. Soriano-Vargas, E.O. Avila-Davila, V.M. Lopez-Hirata, N. Cayetano-Castro, J.L. Gonzalez-Velazquez, Effect of spinodal decomposition on the mechanical behavior of Fe–Cr alloys, Mater. Sci. Eng. A 527 (12) (2010) 2910–2914. [14] Z. Jiao, G. Was, Novel features of radiation-induced segregation and radiationinduced precipitation in austenitic stainless steels, Acta Mater. 59 (3) (2011) 1220–1238.

[15] E. Martínez, O. Senninger, C.C. Fu, F. Soisson, Decomposition kinetics of Fe–Cr solid solutions during thermal aging, Phys. Rev. B 86 (2012) 224109. [16] C. Ortiz, M. Caturla, C.C. Fu, F. Willaime, He diffusion in irradiated a-Fe: An abinitio-based rate theory model, Phys. Rev. B 75 (10) (2007) 100102. [17] C. Björkas, K. Nordlund, M. Caturla, Influence of the picosecond defect distribution on damage accumulation in irradiated a-Fe, Phys. Rev. B 85 (2) (2012) 024105. [18] T. Garnier, M. Nastar, Coarse-grained kinetic Monte Carlo simulation of diffusion in alloys, Phys. Rev. B 88 (13) (2013) 134207. [19] C. Becquart, A. Barbu, J. Bocquet, M. Caturla, C. Domain, C.C. Fu, et al., Modeling the long-term evolution of the primary damage in ferritic alloys using coarsegrained methods, J. Nucl. Mater. 406 (1) (2010) 39–54. [20] G. Bonny, R. Pasianot, L. Malerba, A. Caro, P. Olsson, M. Lavrentiev, Numerical prediction of thermodynamic properties of iron–chromium alloys using semiempirical cohesive models: the state of the art, J. Nucl. Mater. 385 (2) (2009) 268–277. [21] D. Terentyev, I. Martin-Bragado, Evolution of dislocation loops in iron under irradiation: the impact of carbon, Scr. Mater. 97 (2015) 5–8. [22] I. Martin-Bragado, A. Rivera, G. Valles, J. Gomez-Selles, M. Caturla, MMonCa: an object kinetic Monte Carlo simulator for damage irradiation evolution and defect diffusion, Comput. Phys. Commun. 184 (12) (2013) 2703–2710. [23] A. Seeger, The mechanisms of diffusion in metals and alloys, J. Less Common Met. 28 (2) (1972) 387–418. [24] H. Mehrer, N. Stolica, N. Stolwijk, Diffusion in Solids, Springer, 2007. [25] P. Castrillo, R. Pinacho, M. Jaraiz, J. Rubio, Physical modeling and implementation scheme of native defect diffusion and interdiffusion in SiGe heterostructures for atomistic process simulation, J. Appl. Phys. 109 (2011) 103502. [26] L.A. Marqués, L. Pelaz, P. Castrillo, J. Barbolla, Molecular dynamics study of the configurational and energetic properties of the silicon self-interstitial, Phys. Rev. B 71 (8) (2005) 085204. [27] I. Martin-Bragado, P. Castrillo, M. Jaraiz, R. Pinacho, J. Rubio, J. Barbolla, Physical atomistic kinetic Monte Carlo modeling of Fermi-level effects of species diffusing in silicon, Phys. Rev. B 72 (3) (2005) 035202. [28] G.H. Vineyard, Frequency factors and isotope effects in solid state rate processes, J. Phys. Chem. Solids 3 (1) (1957) 121–127. [29] D.T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys. 22 (4) (1976) 403– 434. [30] B. Jönsson, Assessment of the mobilities of Cr, Fe, and Ni in fcc Cr–Fe–Ni alloys, Z. Metallkd. 86 (10) (1995) 686–692. [31] P. Olsson, C. Domain, J. Wallenius, Ab initio study of Cr interactions with point defects in bcc Fe, Phys. Rev. B 75 (2007) 014110. [32] P. Olsson, I. Abrikosov, L. Vitos, J. Wallenius, Ab initio formation energies of FeG (Cr alloys, J. Nucl. Mater. 321 (1) (2003) 84–90. [33] M. Levesque, E. Martínez, C-C. Fu, M. Nastar, F. Soisson, Simple concentrationdependent pair interaction model for large-scale simulations of Fe–Cr alloys, Phys. Rev. B 84 (18) (2011) 184205. [34] G. Bonny, R. Pasianot, D. Terentyev, L. Malerba, Iron chromium potential to model high-chromium ferritic alloys, Philos. Mag. 91 (12) (2011) 1724–1746. [35] J.M. Philibert, Atom Movements-Diffusion and Mass Transport in Solids, EDP Sciences, 2012. [36] G. Bonny, D. Terentyev, L. Malerba, New contribution to the thermodynamics of Fe–Cr alloys as base for ferritic steels, J. Phase Equilib. Diffus. 31 (5) (2010) 439–444. [37] J.O. Andersson, B. Sundman, Thermodynamic properties of the Cr Fe system, Calphad 11 (1) (1987) 83–92. [38] H. Kuwano, Mossbauer effect study on the miscibility gap of the iron– chromium binary system, Transa. Jpn. Inst. Met. 26 (7) (1985) 473–481. [39] D. Chen, A. Kimura, W. Han, Correlation of Fe/Cr phase decomposition process and age-hardening in Fe–15Cr ferritic alloys, J. Nucl. Mater. 455 (1) (2014) 436–439. [40] S. Brenner, M. Miller, W. Soffa, Spinodal decomposition of iron-32 at.% chromium at 470 C, Scr. Metall. 16 (7) (1982) 831–836.