1.18.4.2.1 Dilute alloy models

RIS measurements in dilute alloys are less numerous than in concentrated alloys because the required grain boundary concentrations are usually smaller. However, some of the first RIS observations concerned the seg­regation of a dilute element, Si, in austenitic steels. In this specific case, observations were easy because the RIS of Si was accompanied by precipitation of Ni3Si.

The first mechanism that was proposed to explain the observed solute segregation was the diffusion of solute-point defect complexes towards sinks.105 Since then, more rigorous models that rely on the linear response theory have been established and applied to the RIS description of Mn and P in nickel108 and P in ferritic steels.87,107 Although good precision of the microscopic parameters was still missing, the formu­lation of the kinetic equations was general enough to be used almost without modification.105 Recent ab initio calculations not only provided accurate atomic jump frequencies of P in Fe,7,70 but they also called into question the jump interstitial diffusion mechanism that had to be considered.7 Indeed, the octahedral and the (110) mixed dumbbell configura­tions have almost the same stability and migration enthalpies. The resulting effective diffusion energy estimated by the transport model was found to be smaller than the self-interstitial atom migration enthalpy, confirming the classical statement that a solute atom with a negative size effect tends to segre­gate at the grain boundary. However, as emphasized in Meslin et a/.,7 the current interpretation of the interstitial contribution to RIS in terms of size effects is certainly oversimplified. A very large ab initio value of 1.05 eV for the binding energy between a mixed dumbbell and a substitutional P atom may lead to a large activation energy for P interstitials and a drastic reduction of P segregation predictions.7 To consider this new blocking configuration with two P atoms, a concentrated alloy diffusion model including short — range-order effects is required.

It is interesting to note that the same solute seems to have a positive coupling with the vacancy also (although the calculation was not as precise as for the interstitials as it was based on an empirical poten — tial).117 In the same way, recent ab initio calculations showed that a Cu solute is also expected to be dragged by vacancy at low temperatures in Fe.11 ,119

In the process ofsetting up an ion irradiation experi­ment, a number of parameters that involve beam

characteristics (energy, current/dose) and beam — target interaction must be considered. ASTM E 521 provides standard practice for neutron radiation damage simulation by charged-particle irradiation49 and ASTM E 693 provides standard practice for characterizing neutron exposures in iron and low alloy steels in units of dpa.9 One of the most important considerations is the depth of penetration. Figure 35 shows the range versus particle energy for protons, helium ions, and nickel ions in stainless steel as calculated by SRIM.50 The difference in penetration depth between light and heavy ions is over an order of magnitude in this energy range. Figure 36 shows how several other parameters describing the target

behavior during proton irradiation vary with energy, dose rate, the time to reach 1 dpa, deposited energy, and the maximum permissible beam current (which will determine the dose rate and total dose), given a temperature limitation of 360 °C. With increasing energy, the dose rate at the surface decreases because of the drop in the elastic scattering cross-section (Figure 36(a)). Consequently, the time to reach a target dose level, and hence the length of an irradia­tion, increases rapidly (Figure 36(b)). Energy depo­sition scales linearly with the beam energy, raising the burden of removing the added heat in order to control the temperature of the irradiated region (Figure 36(c)). The need to remove the heat due to higher energies will limit the beam current at a specific target temperature (Figure 36(d)), and a limit on the beam current (or dose rate) will result in a longer irradiation to achieve the specified dose. Figure 37 summarizes how competing features of an irradiation vary with beam energy, creating trade­offs in the beam parameters. For example, while greater depth is generally favored in order to increase the volume of irradiated material, the higher energy required leads to lower dose rates near the surface and higher residual radioactivity. For proton irradia­tion, the optimum energy range, achieved by balanc­ing these factors, lies between 2 and 5 MeV as shown by the shaded region.

1.10.1 Basics

Materials relevant to reactors are held together by electrons. An interatomic potential expresses energy in terms of atomic positions. the electronic and mag­netic degrees of freedom are integrated out.

Most empirical potentials are derived on the basis of some approximation to quantum mechanical energies. If they are subsequently used in MD of solids, then what is actually used are forces: the derivative of the energy. For near-harmonic solids, it is actually the second derivative of the energy that governs the behavior.

A dilemma: Does the primary term covering ener­getics also dominate the second derivative of the energy? To take an extreme example, an equation which calculates the energy of a solid exactly for all configurations to 0.1% is: E = mc2 : most of the energy is in the rest mass of the atoms. But this is patently useless for calculating condensed matter properties. We encounter the same problem in a less extreme form in metals: should we concentrate the energy gained in delocalizing the electrons to form the metal, or is treating perturbations around the metallic state more useful? In general, the issue is ‘What is the reference state.’ Most potentials implic­itly assume that the free atom is the reference state.

As discussed elsewhere,59,63,65 in-cascade vacancy clustering in iron is quite low (~10% of NRT) when a NN criterion for clustering is applied. This was identified as one of the differences between iron and copper in the comparison of these two materials reported by Phythian and coworkers.59 However, when the coordinates of the surviving vacancies in 10, 20, and 40keV cascades were analyzed, clear spatial correlations were observed. Peaks in the distributions of vacancy-vacancy separation dis­tances were obtained for the second and fourth NN locations.64 These radial distributions are shown in Figure 19. Similar results were obtained from the analysis of the vacancy distributions in higher energy cascades at 100 and 600 K. The peak observed for vacancies in second NN locations is consistent with the di-vacancy binding energy being greater for second NN (0.22 eV) than for first NN (0.09 eV).90 The reason for the peak at fourth NN is presumably related to this also since two vacancies that are second NN to a given vacancy would be fourth NN. In addi­tion, work discussed by Djurabekova and coworkers91 indicates that there is a small binding energy between two vacancies at the fourth NN distance.

Sonic
shock wave "V
(nondestructive)
separation ^

Transonic shock-= Spaghetti

wave limit

(d)

Figure 18 Schematic representation of cascade development leading to the formation of interstitial and vacancy clusters formation. Reproduced from Calder, A. F.; Bacon, D. J.; Barashev, A. V.; Osetsky, Yu. N. Phil. Mag. 2010, 90, 863-884.

An example of a locally vacancy-rich region in a 50keV, 100 K cascade is shown in Figure 20, where the region around a collection of 14 vacancies has

vacancy within the fourth NN spacing of 1.66a0, where a0 is the iron lattice parameter. The ‘cluster’ is shown in two views: a 3D perspective view and an orthographic projection (—x) in Figure 20. Such an arrangement of vacancies is similar to some of the vacancy clusters observed by Sato and coworkers in field ion microscope images of irradiated tungsten.9 Since the time period of the MD simulations is too short to allow vacancies to jump (<100 ps), the possi­bility that these closely correlated vacancies might collapse into clusters over somewhat longer times has

been investigated using MC simulations. The vacancy coordinates at the end ofthe MD simulations were extracted and used as the starting configuration in MC cascade annealing simulations. The expecta­tion of vacancy clustering was confirmed in the MC simulations, where many of the isolated vacancies had clustered within 70 ms.90,

The energy and temperature dependence of in-cascade vacancy clustering as a fraction of the NRT displacements is shown in Figure 21 for cas­cade energies of 10-50 keV. Results are shown for clustering criteria of first, second, third, and fourth NN. A comparison of Figure 21 and Figure 12 demonstrates that in-cascade vacancy clustering in iron remains lower than that of interstitials even when the fourth NN criterion is used. This is consis­tent with the experimentally observed difficulty of forming visible vacancy clusters in iron as discussed by Phythian and coworkers,59 and the fact that only relatively small vacancy clusters are found in positron annihilation studies of irradiated ferritic alloys.94 However, it should be pointed out that work with more recently developed iron potentials finds less dif­ference between vacancy and interstitial clustering.74 The cascade energy dependence of vacancy cluster­ing is similar to that of interstitials; there is essentially zero clustering at the lowest energies but it rapidly increases with cascade energy and is relatively inde­pendent of energy above ^10 keV. However, vacancy clustering decreases as the temperature increases,

which is consistent with vacancy clusters being ther­mally unstable.

Fractional vacancy cluster size distributions are shown in Figure 22, for which the fourth NN clus­tering criterion has been used. Figure 22(a) illus­trates that the vacancy cluster size distribution shifts to larger sizes as the cascade energy increases from 10 to 50keV. This is similar to the change shown for interstitial clusters in Figure 13(a). There is a corresponding reduction in the fraction of single vacancies. However, as mentioned above, the effect of cascade temperature shown in Figure 22(b) and 22(c) is the opposite of that observed for interstitials. The magnitude of the temperature effect on the vacancy cluster size distributions also appears to be weaker than in the case of interstitial clusters. The fraction of single vacancies increases and the size dis­tribution shifts to smaller sizes as the temperature increases from 100 to 900 K for the 10keV cascades, and from 100 to 600 K for 20 keV cascades. Similar to the case of interstitial clusters, the effect of tempera­ture seems to be greater at 20 keV than at 10 keV.

The main idea of the grouping methods for numerical evaluation of the ME is to replace a group of equations described by the ME with an ‘averaged’ equation. Such a procedure was proposed by Kiritani75 for describing the evolution ofvacancy loops during aging of quenched metals. Koiwa84 was the first to examine the Kiritani method by com­paring numerical results with the results of an analytical solution for a simple problem. Serious disagreement was found between the numerical and analytical results, raising strong doubts regard­ing the applicability of the method. The main objec­tion to the method75 in Koiwa84 is the assumption used by Kiritani75 that the SDF within a group does not depend on the size of clusters. However, Koiwa did not provide an explanation of where the inaccuracy comes from. The Validity of the Kiritani method was examined thoroughly by Golubov et a/.85 The general conclusion of the analysis is that the grouping method proposed by Kiritani is not accurate. The origin of the error is the approxi­mation that the SDF within a group is constant as was predicted by Koiwa.84 Thus, the disagree­ment found in Koiwa84 is fundamental and cannot be circumvented. Because it is important for under­standing the accuracy of the other methods sug­gested for numerical calculations of cluster evolution,
the analysis performed in Golubov et a/.85 is briefly highlighted below.

It follows from eqn [18] that the total number of clusters, N(t) = ^2T=2 f (x, t) and total number of defects in the clusters, S(t) = ^Ц=2 xf (x, t), are described by the following equations:

f = J 0.t) N

@ S 1

— = J(1, t)+ J(x, t) [37]

x=1

where the generation term in eqn [18] is dropped for simplicity. Equations [36] and [37] are the con­servation laws which can be satisfied when one uses a numerical evaluation of the ME. When a group method is used, the conservation laws can be satisfied for reactions taking place within each group.69 How­ever, this is not possible within the approximation used by Kiritani75 because a single constant can be used to satisfy only one of the eqns [36] and [37]. To resolve the issue, Kiritani75 used an ad hoc modifica­tion of the flux J(x,); therefore, the final set of equations for the density of clusters within a group, Fj, are as follows

dF, dt
Jj ]	[38]
It is worth noting that this comparison also sheds light on the relative accuracy of other numerical solutions of the F—P equation such as in Bondarenko and Konobeev,76 Ghoniem and Sharafat,77 Stoller and Odette,78 and Hardouin Duparc et a/.79 Equations [36] and [37] provide a way of getting a simple but still reasonably correct grouping method for numerical integration of the ME. Indeed, the two conservation laws, eqns [36] and [37], require two parameters within a group at least. The simplest approximation of the SDF within a group of clusters (sizes from x,—1 to x, = x,—1 + Ax, — 1) can be achieved using a linear function
[39]
Jj
where Ax, is the width of the ‘j’ group. Equations [38] and [39] indeed satisfy both the conservation laws. However, they do not provide a correct description of cluster evolution described by the ME because the flux Jj in eqn [39] depends on the widths of groups and these widths have no physical meaning. An example of a comparison of the calcu­lation results obtained using the Kiritani method with the analytical and numerical calculations based on a more precise grouping method is pre­sented in Figure 3. Note that in the limiting case where the widths of group are equal, Ax, = Ax, +1, the flux Jj is equal to the original one, J(x, t). In this limiting case, eqns [38] and [39] correspond to those that can be obtained by a summation of the ME within a group and therefore they provide con­servation of the total number of clusters, N(t), only. This limiting case is probably the simplest way to demonstrate the inaccuracy of the Kiritani method.
f (x) = L0(x — (x), ) + L1 N where (x)j = x, — 1/2 (Ax, — 1) is the mean size of the group. Equations for L0, L] are as follows69
dI0 dt
[J (x,—1)—Jx (x,)]
[41]
di; dt —<xi—‘’+j <x’)—■J (x’)—0}|42]

This chapter has attempted to illustrate the power of KMC for modeling radiation effects in structural materials and nuclear fuels, following an introduction and review of the Monte Carlo technique. Monte Carlo modeling, first developed by Metropolis and coworkers9,10 during the Manhattan project does pro­vide a physically satisfying technique to simulate the stochastic evolution of defect evolution in materials science and in fact has been used to simulate irradia­tion effects on materials for four decades. There are three main types of KMC modeling used in irradia­tion effects, namely event Monte Carlo, object Monte Carlo, and atomistic Monte Carlo.

This chapter has focused on describing the atom­istic KMC and OKMC methods by providing two examples of successful KMC simulations to predict the coupled evolution of vacancy clusters and copper precipitates during low dose rate neutron irradiation of Fe-Cu alloys and the transport and diffusional release of the fission product, silver, in TRISO nu­clear fuel. These examples clearly demonstrate the power and ability of KMC models to capture the spatial correlations that can be an important compo­nent of microstructural evolution in nuclear materi­als. Yet, the further widespread application of KMC models will require algorithmic developments that can more readily treat the wide range of time scales inherent in microstructural evolution and yet effec­tively incorporate the rare-event dynamics in inte­grating system performance to realistic time and irradiation dose exposures.

Thus, the challenges that must be overcome in future nuclear materials modeling include:

• bridging the inherently multiscale time and length scales which control materials degradation in nu­clear environments;

• dealing with the complexity of multicomponent materials systems, including those in which the chemical composition is continually evolving due to nuclear fission and transmutation;

• discovering the unknown to prevent technical surprises;

• transcending ideal materials systems to engineer materials and components; and

• incorporating error assessments within each mod­eling scale and propagating the error through the scales to determine the appropriate confidence bounds on performance predictions.

Successful overcoming of these challenges will result in nuclear materials performance models that can predict the properties, performance, and lifetime of nuclear fuels, cladding, and components in a vari­ety of nuclear reactor types throughout the full life cycle, and provide the scientific basis for the compu­tational design of advanced new materials. While the current chapter is focused on the KMC modeling methodology, it is important to note the challenges of predictive materials models of irradiation effects. High performance computing at the petascale, exas — cale, and beyond is a necessary and indeed critical tool in resolving these challenges, yet it is important to realize that exascale computing on its own will not be sufficient. This is best recognized from a simple example considering the computational degrees of freedom in a MD simulation. Assuming that reliable, multicomponent interatomic potentials existed for the nuclear fuel rod and cladding in a nuclear power plant and that a constant time-step of 2 x 10~15s could sufficiently capture the physics of high-energy atomic collisions to conserve energy; then to simulate 1 day of evolution of 1 cm tall, 1 cm diameter fuel pellet clad and zirconium clad would require ^6 x 1022 atoms for ^4 x 1019 time — steps. For comparison, the LAMPPS MD code using classical force fields has been benchmarked with 40 billion atoms (4 x 1010) and 100 time-steps on 10 000 processors of the RedStorm at Sandia National Laboratory with a wall clock time of 980 s and on 64 000 processors of the BlueGene Light at Lawrence Livermore National Laboratory with a wall clock time of 585 s.103 Thus, even assuming optimistic scaling and parallelization, brute force atomistic simulation of the first full power day of a nuclear fuel pellet in a reactor by MD will remain well beyond the reach of high performance comput­ing capabilities for the next decade.

77. .

1.18.1 Introduction

Irradiation creates excess point defects in materials (vacancies and self-interstitial atoms), which can be eliminated by mutual recombination, clustering, or annihilation of preexisting defects in the micro­structure, such as surfaces, grain boundaries, or dis­locations. As a result, permanent irradiation sustains fluxes of point defects toward these point defect sinks and, in case of any preferential transport of one of the alloy components, leads to a local chemical redistribution. These radiation-induced segregation (RIS) phenomena are very common in alloys under irradiation and have important technological implica­tions. Specifically in the case of austenitic steels, because Cr depletion at the grain boundary is sus­pected to be responsible for irradiation-assisted stress corrosion, a large number of experiments have been conducted on the RIS dependence on alloy composi­tion, impurity additions, irradiation flux and time, irradiation particles (electrons, ions, or neutrons), annealing treatment before irradiation, and nature of grain boundaries.1-5

The first RIS models generally consisted of appli­cation of Fick’s laws to reproduce two specific effects of irradiation: diffusion enhancement due to the increase ofpoint defect concentration, and the driving forces associated with point defect concentration gra­dients. According to these models, RIS is controlled by kinetic coefficients D or L (defined below) relating atomic fluxes to gradients of concentration or chemi­cal potentials. It was shown that these coefficients are best defined in the framework of the thermodynamics of irreversible processes (TIPs) within the linear response theory. RIS models were then separated into two categories: models restricted to dilute alloys, and models developed for concentrated alloys.

From the beginning until now, the dilute alloy models have benefited from progress made in the diffusion theory.6 The explicit relations between the phenomenological coefficients L and the atomic jump frequencies have been established, at least for alloys with first nearest neighbor (nn) interactions. In principle, such relations allow the immediate use of ab initio atomic jump frequencies and lead to pre­dictive RIS models.7

While the progress of RIS models of dilute alloys is closely related to that of diffusion theory, most segregation models for concentrated alloys still use oversimplified diffusion models based on Manning’s relations.8 This is mainly because the jump sequences of the atoms are particularly complex in a multicomponent alloy on account of the multiple jump frequencies and correlation effects that are involved. Only very recently has an interstitial diffu­sion model been developed that could account for short-range order effects, including binding energies with point defects.9, Emphasis has so far been

placed on comparisons with experimental observations. The continuous RIS models have been modified to include the effect of vacancy trapping by a large-sized impurity or the nature and displacement of a specific grain boundary. Most of the diffusivity coefficients of Fick’s laws are adjusted on the basis of tracer diffusion data. Paradoxically, the first RIS models were more rigorous11 than the present ones in which thermody­namic activities, particularly some of the cross-terms, are oversimplified. In this review, we go back to the first models starting from the linear response theory, albeit slightly modified, to be able to reproduce the main characteristics of an irradiated alloy. It is then possible to rely on the diffusion theories developed for concen­trated alloys.

Then again, lattice rate kinetic techniques12-14 and atomic kinetic Monte Carlo (AKMC) methods15-17
have become efficient tools to simulate RIS. Thanks to a better knowledge of jump frequencies due to the recent developments of ab initio calculations, these simulations provide a fine description of the thermo­dynamics as well as the kinetics of a specific alloy. Moreover, information at the atomic scale is precious when RIS profiles exhibit oscillating behavior and spread over a few tens of nanometers.

Discoveries and typical observations of RIS are illustrated in the first section. In the second section, the formalism of TIP is used to write the alloy flux couplings. It is explained that fluxes can be estimated only partially from diffusion experiments and ther­modynamic data. An alternative approach is the cal­culation of fluxes from the atomic jump frequencies. The third section presents more specifically the con­tinuous RIS models separated into the dilute and concentrated alloy approaches. The last section intro­duces the atomic-scale simulation techniques.

Vacancies and interstitials, with the associated forma­tion energy driving their concentration and migra­tion energy driving their displacement in the solid; the sum of these two energies is the activation energy for diffusion at equilibrium. For such simple defects, it is possible to go beyond the 0 K energies and to access the free energies of formation and migrations by calculating the vibrational spectra in the presence of the defect in the stable position and at the saddle point (see Section 1.08.4.2.3).

1.08.3.2.1 Hetero-defects

In the nuclear context, such defects can be fission products in a fuel material, actinide atoms in a waste material, helium gases in structural materials, and so on; ab initio gives access to the solution energy of these impurities, which allows one to determine their most favored positions in the crystal: interstitial position, substitution for host atoms, and so on. The kinetic energies of migration of interstitial impurities are accessible as well as the kinetic barrier for the extraction of an impurity from a vacancy site.

In the approach described earlier, which considers low concentrations of solutes and defects, the number of independent configurations is rather small, and they can be easily taken into account in kinetics model. The situation is much more complex when considering Fe-Cr with Cr concentration in the range 10-20%. Nevertheless, first results have been obtained by considering the interaction of defects with one or two Cr atoms in the Fe matrix.81 These data could ideally be used to fit an improved empiri­cal potential, but the Fe-Cr system is rather difficult to model because of the strong interplay between magnetic and chemical interactions. This is also clearly one of the challenges in the field.

1.08.4.2.2
Point defects in hcp-Zr

Point defects in hcp-Zr have also been studied via DFT calculations. It was found in particular that the vacancy migration energy is lower by ~0.15 eV within the basal plane than out of the basal plane.82 The situation for the self-interstitial is quite complex, since among the known configurations, at least three configurations are found to have almost the same formation energy (within 0.1 eV): the octahedral (O), split dumbbell (S), and basal octahedral (BO) configurations.83,84