1.16.3.4.1 Background

It is now well accepted that fracture of ferritic steels in the temperature range where it propagates by cleav­age originates in microcracks (mostly in precipitates) ahead of macrocracks, which could be precracks in a test specimen or surface cracks in structures. To explain the extremely high cleavage fracture tough­ness (>20 MPa Vm) compared with that («1 MPa Vm) calculated from surface energy alone assuming pure cleavage of Fe matrix, it was originally postulated by Orowan41 that fracture in ferritic steels could occur due to cleavage originated in microcracks situated ahead of the main crack. Later, it was found by experiments that these microcracks originate in precipitates,42,43 and that the propagation of these microcracks into the matrix was assumed to be the controlling step in the fracture of ferritic steels. Another observation, though less well established, is that the cleavage stress at fracture on these microcracks is invariant with

44,45

temperature.

Ritchie, Knott, and Rice (RKR)46 used the HRR solutions (Hutchinson, Rice, and Rosengren) and finite element analysis (FEA) to simulate the plastic zone, and used a critical tensile stress achieved over a characteristic distance ahead of the crack as the fail­ure criterion. This distance is essentially a fitting parameter, and RKR46 used a value equal to or twice the average the grain diameter. The model successfully predicts the lower-shelf fracture tough­ness, but fails to predict the upturn near the transition temperature. Statistical models were introduced to predict the brittle-ductile transition (BDT) in steels starting with Curry and Knott,43 most notable among them were by Beremin47 and Wallin et a/.48 In both of these models, FEA solutions of crack-tip plasticity were used to obtain the stress fields ahead of the crack. In the Beremin model,47 the maximum princi­pal stress is calculated for each volume element in the plastic zone and a probability of failure is assigned. The total probability of failure is then obtained by summing over the entire plastic zone. Wallin et a/.48 extended the modeling to the transition region by considering variation of the effective surface energy (gs + gp) with temperature, where gs is the true surface energy and gp the plastic work done during propaga­tion. This eventually led to the master curve (MC) hypothesis, which predicts that the BDT of all ferritic steels follows a universal curve.48,49 Even though the MC is used to check the reliability of structures under irradiation,50 a clear understanding of the physical basis of this methodology is still lacking.51 Odette and He52 explained the MC using a micro­scopic fracture stress varying with temperature. Most experimental findings53,54 indicate that the fracture stress is not sensitive to temperature, and more care­ful experiments and simulations may be required to resolve this issue.

Discrete dislocation simulations of crack tips were successful in predicting the BDT of simple single crystalline materials.55,56 The advantage of this approach over the continuum methods is that funda­mental material properties such as dislocation veloc­ity and their mutual interactions can be treated dynamically. By these simulations, it has been found that the dislocation mobility plays a significant role in determining the transition temperature.56 However, the variation of dislocation mobility alone cannot explain the BDT behavior. An earlier attempt to model the BDT of complex materials like steels57 predicted the lower-shelf fracture toughness, substan­tiating Orowan’s postulate41 of high fracture tough­ness measured at low temperature. However, the model failed to predict the sharp increase of fracture toughness around the transition temperature region. Here, we present a discrete dislocation simulation in which crack-tip blunting is accounted for the first time. The effects of blunting are incorporated in the simulation using elastic stress fields of blunted cracks. As the crack tip is blunted due to dislocation emission, and the position of the ‘virtual sharp crack tip’ retreats from the blunted tip thereby reducing the field at the microcrack further in addition to the contribution from emitted dislocations. The critical particle is assumed to be at a fixed distance from the blunted macrocrack tip.

Because of these difficulties, simulations of diffusive phase transformation kinetics are commonly based on various broken-bond models, in the framework of rigid lattice approximations.5,76 The total energy of the system is considered to be a sum of constant pair interaction energies, for example, e^B between A and

B atoms located on nth nn sites. Interactions between atoms and point defects can also be used to provide a better description of their formation energies and interactions with solute atoms, and other defects.

Various approximations are used to compute the migration barriers: a common one77 is writing the saddle-point energy of the system as the mean energy between the initial energy E and final energy Ef, plus a constant contribution Q(which can depend on the jumping atom, A or B). The migration barrier for an A-V exchange is then:

DEAVg = El——El + Qa [22]

where Ef — E corresponds to the balance of bonds destroyed and created during the exchange.

Another solution is to explicitly consider the interaction energy eAV of the jumping atom A with the system, when it is at the saddle point:

ifA7 = eAV — E *A? — E eVj [23]

i, n j, n

eAV itself can be written as a sum of interactions between A and the neighbors of the saddle point.12,78,79

Both approximations are easily extended to inter­stitial diffusion mechanisms, and their parameters can be fitted to experimental data and/or ab initio calculations. The first one has the drawback of imposing a linear dependence between the barrier and the difference between the initial and final ener­gies, which is not justified and has been found to be unfulfilled in the very few cases where it has been checked72 (with empirical potentials). The second one should better take into account the effect of the local configuration and, according to the theory of activated processes, does not impose a dependence of the barrier on the final state. However, a model of pair interactions on a rigid lattice does not give a very precise description of the energetic landscape in a metallic solid solution, so, the choice of approxima­tions [22] or [23] may not be crucial. Taking into account many-body interactions (fitted to ab initio calculations, using cluster expansion methods) could improve the description of migration barriers, but would significantly increase the simulation time.80,81

In this paragraph, we try to give some indication of what can be done with an ab initio code and how it is done in practice. The calculation starts with the posi­tioning of atoms of given types in a calculation cell of a certain shape. That would be all if the calculations were truly ab initio. Unfortunately, a few more pieces of information should be passed to the code; the most important ones are described in the final section. The first section introduces the basic outputs of the code, and the second one deals with the possible cell sizes and the associated CPU times.

1.08.2.3.1 Output

We describe in this section the output of ab initio calculations in general terms. The possible applica­tions in the nuclear materials field are given below. The basic output of a standard ab initio calculation is the complete description of the electronic ground state for the considered atomic configuration. From this, one can extract electronic as well as energetic information.

On the electronic side, one has access to the elec­tronic density of states, which will indicate whether the material is metallic, semiconducting, or insulating (or at least what the code predicts it to be), its possi­ble magnetic structure, and so on. Additional calcula­tions are able to provide additional information on the electronic excitation spectra: optical absorption, X-ray spectra, and so on.

On the energetic side, the main output is the total energy of the system for the given atomic configura­tion. Most codes are also able to calculate the forces acting on the ions as well as the stress tensor acting on the cell. Knowing these forces and stress, it is possible to chain ground-state calculations to perform various calculations:

• Atomic relaxations to the local minimum for the

atomic positions.

• From the relaxed positions (where forces are zero), one can calculate second derivatives of the energy to deduce, among other things, the phonon spec­trum. This can be done either directly, by the so — called frozen phonon approach, or by first-order perturbation theory (if such feature is implemen­ted in the code). In this last case, the third-order derivative of the energy (Raman spectrum, phonon lifetimes) can also be computed.

• Starting from two relaxed configurations close in space, one can calculate the energetic path in space joining these two configurations, thus allowing the calculation of saddle points.

• The integration of the forces in a Molecular Dynamics scheme leads to so-called ab initio molecular dynamics (see Chapter 1.09, Molecu­lar Dynamics). Car-Parrinello molecular dynam­ics18 calculations, which pertain to this class of calculations, introduce fictitious dynamics on the electrons to solve the minimization problem on the electrons simultaneously with the real ion dynamics.

Starting from the perfect crystal at equilibrium lattice constant a0, we can assign initial velocities to the atoms and perform MD simulations. In the simplest simulation, no thermostat is introduced to regulate the temperature, and no barostat is introduced to regulate the stress. The simulation then corresponds to the NVE ensemble, where the number of particles N, the cell volume V (as well as shape), and total energy E are conserved. This simulation is usu­ally performed as a benchmark to ensure that the numerical integrator is implemented correctly and that the time step is small enough.

The instantaneous temperature Tinst is defined in terms of the instantaneous kinetic energy K through the relation K = (3N/2)kBTinst, where kB is Boltzmann’s constant. Therefore, the velocity can be initialized by assigning random numbers to each component of every atom and scaling them so that Tinst matches the desired temperature. In practice, Tinst is usually set to twice the desired temperature for MD simulations of solids, because approximately half of the kinetic energy flows to the potential energy as the solids reach thermal equilibrium. We also need to subtract appropriate constants from the x, y z components of the initial velocities to make sure the center-of-mass linear momentum of the entire cell is zero. When the solid contains surfaces and is free to rotate (e. g., a nanoparticle or a nanowire), care must be taken to ensure that the center-of-mass angular momentum is also zero.

Figure 7(a) plots the instantaneous temperature as a function of time, for an MD simulation starting with a perfect crystal and Tinst = 600 K, using the Velocity Verlet integrator13 with a time step of A t = 1 fs. After 1 ps, the temperature of the simulation cell is equilibrated around 300 K. Due to the finite time step At, the total energy E, which should be a conserved quantity in Hamiltonian dynamics, fluctuates during the MD simulation. In this simula­tion, the total energy fluctuation is <2 x 10-eV per atom, after equilibrium has been reached (t> 1 ps). There is also zero long-term drift of the total energy. This is an advantage of symplectic integrators11, and also indicates that the time step is small enough.

The stress of the simulation cell can be computed by averaging the Virial stress for time between 1 and 10 ps. A hydrostatic pressure P = — (sxx + + ozz)/3 =

1.33 ± 0.01GPa is obtained. The compressive stress develops because the crystal is constrained at the zero-temperature lattice constant. A convenient way to find the equilibrium lattice constant at finite tem­perature is to introduce a barostat to adjust the vol­ume of the simulation cell. It is also convenient to introduce a thermostat to regulate the temperature of the simulation cell. When both the barostat and thermostat are applied, the simulation corresponds to the NPT ensemble.

The Nose-Hoover thermostat11,33,34 is widely used for MD simulations in NVT and NPT ensem­bles. However, care must be taken when applying it to perfect crystals at medium-to-low temperatures, in which the interaction between solid atoms is close to harmonic. In this case, the Nose-Hoover thermo­stat has difficulty in correctly sampling the equi­librium distribution in phase space, as indicated by periodic oscillation of the instantaneous temperature.

2	2
03 0
0 20 40 60 t (Ps)
0
80	100
(b)

Figure 7 (a) Instantaneous temperature Tinst and Virial pressure p as functions of time in an NVE simulation with initial temperature at 600 K. (b) Tinst and P in a series of NVT at T = 300 K, where the simulation cell length L is adjusted according to the averaged value of P.

The Nose-Hoover chain35 method has been devel­oped to address this problem.

The Parrinello-Rahman19 method is a widely used barostat for MD simulations. However, periodic oscillations in box size are usually observed during equilibration of solids. This oscillation can take a very long time to die out, requiring an unreasonably long time to reach equilibrium (after which meaningful data can be collected). A viscous damping term is usually added to the box degree of freedom to accel­erate the speed of equilibration. Here, we avoid the problem by performing a series of NVT simulations, each one lasting for 1 ps using the Nose-Hoover chain method with Velocity Verlet integrator and At = 1 fs. Before starting each new simulation, the simulation box is subjected to an additional hydro­static elastic strain of e = (P)/B0, where (P) is the average Virial pressure of the previous simulation, where B = 2000GPa is an empirical parameter.

The instantaneous temperature and Virial pressure during 100 of these NVT simulations are plotted in Figure 7(b). The instantaneous temperature fluctuates near the desired temperature (300 K) nearly from the beginning of the simulation. The Virial pressure is well relaxed to zero at t = 20 ps. The average box size from 50 to 100 ps is L = 16.5625 A, which is larger than the initial value of 16.5290 A. This means that the normal strain caused by thermal expansion at 300 K is exx = 0.00203. Hence, the coefficient of thermal expan­sion is estimated to be a = exx/T = 6.8 x 10_6K_1.35

MD simulations have been employed to investigate displacement cascade evolution in a wide range of materials. The literature is sufficiently broad that any list of references will be necessarily incomplete; Malerba,41 Stoller,43 and others52-70 provide only a representative sample. Additional references will be given below as specific topics are discussed. The recent review by Malerba41 provides a good summary of the research that has been done on iron. These MD investigations of displacement cascades have estab­lished several consistent trends in primary damage formation in a number of materials. These trends include (1) the total number of stable point defects produced follows a power-law dependence on the cas­cade energy over a broad energy range, (2) the ratio of MD stable displacements divided by the number obtained from the NRT model decreases with energy until subcascade formation becomes prominent, (3) the in-cascade clustering fraction of the surviving defects increases with cascade energy, and (4) the effect of lattice temperature on the MD results is rather weak. Two additional observations have been made regarding in-cascade clustering in iron, although the fidelity of these statements depends on the interatomic potential employed. First, the interstitial clusters have a complex, three-dimensional (3D) morphology, with both sessile and glissile configura­tions. Mobile interstitial clusters appear to glide with a low activation energy similar to that of the mono­interstitial (~0.1—0.2 eV).71 Second, the fraction of the vacancies contained in clusters is much lower than the interstitial clustering fraction. Each of these points will be discussed further below.

The influence of the interatomic potential on cascade damage production has been investigated by several researchers.42,72-74 Such comparisons generally show only minor quantitative differences between results obtained with interatomic potentials of the same general type, although the differences in cluster­ing behavior are more significant with some potentials. Variants of embedded atom or Finnis-Sinclair type potential functions (see Chapter 1.10, Interatomic Potential Development) have most often been used.

However, more substantial differences are sometimes observed that are difficult to correlate with any known aspect ofthe potentials. The analysis recently reported by Malerba41 is one example. In this case, it appears that the formation of replacement collision sequences (RCS) (discussed in Section 1.11.4.1) was very sensi­tive to the range over which the equilibrium part of the potential was joined to the more repulsive pair potential that controls short-range interactions. This changed the effective cascade energy density and thereby the number of stable defects produced.

Therefore, in order to provide a self-consistent database for illustrating cascade damage production over a range of temperatures and energies and to provide examples of secondary variables that can influence this production, the results presented in this chapter will focus on MD simulations in iron using a single interatomic potential. This

potential was originally developed by Finnis and Sinclair21 and later modified for cascade simulations by Calder and Bacon.58 The calculations were carried out using a modified version of the MOLDY code written by Finnis.75 The computing time with this code is almost linearly proportional to the number of atoms in the simulation. Simulations were carried out using periodic, Parrinello—Rahman boundary condi­tions at constant pressure.76 As no thermostat was applied to the boundaries, the average temperature of the simulation cell was increased as the kinetic energy of the PKA was dissipated. The impact of this heating appears to be modest based on the observed effects of irradiation temperature discussed below, and on the results observed in the work of Gao and coworkers.77 A brief comparison of the iron cas­cade results with those obtained in other metals will be presented in Section 1.11.5.

The primary variables studied in these cascade simulations is the cascade energy, EMD, and the irra­diation temperature. The database of iron cascades includes cascade energies from near the displacement threshold (~100 eV) to a 200 keV, and temperatures in the range of 100—900 K. In all cases, the evolution of the cascade has been followed to completion and the final defect state determined. Typically this is reached after a few picoseconds for the low-energy cascades and up to ~15 ps for the highest energy cascades. Because of the variability in final defect production for similar initial conditions, several simu­lations were conducted at each energy to produce statistically meaningful average values. The para­meters of most interest from these studies are the number of surviving point defects, the fraction of these defects that are found in clusters, and the size distribution of the point defect clusters. The total number of point defects is a direct measure of the residual radiation damage and the potential for long — range mass transport and microstructural evolution. In-cascade defect clustering is important because it can promote microstructural evolution by eliminat­ing the cluster nucleation phase.

The parameters used in the following discussion to describe results of MD cascade simulations are the total number of surviving point defects and the fraction of the surviving defects contained in clus­ters. The number of surviving defects will be expressed as a fraction of the NRT displacements listed in Table 1, whereas the number of defects in clusters will be expressed as either a fraction of the NRT displacements or a fraction of the total surviv­ing MD defects. Alternate criteria were used to define a point defect cluster in this study. In the case of interstitial clusters, it was usually deter­mined by direct visualization of the defect struc­tures. The coordinated movement ofinterstitials in a given cluster can be clearly observed. Interstitials bound in a given cluster were typically within a second nearest-neighbor (NN) distance, although some were bound at third NN. The situation for vacancy clusters will be discussed further below, but vacancy clustering was assessed using first, sec­ond, third, and fourth NN distances as the criteria. The vacancy clusters observed in iron tend to not exhibit a compact structure according to these defi­nitions. In order to analyze the statistical variation in the primary damage parameters, the mean value (M), the standard deviation about the mean (s), and the standard error of the mean (e) have been calculated for each set of cascades conducted at a given energy and temperature. The standard error of the mean is calculated as e = ff/я0 5, where n is the number of cascade simulations completed.78 The standard error of the mean provides a measure of how well the sample mean represents the actual mean. For example, a 90% confidence limit on the mean is obtained from 1.86e for a sample size of nine.79 These statistical quantities are summarized in Table 2 for a representative subset of the iron cascade database.

Table 2 Statistical analysis of primary damage parameters derived from MD cascade simulations

Energy (keV) Temperature (K) Number of cascades Surviving MD Clustered interstitials (mean /

displacements standard deviation / standard error)

(mean / standard deviation / standard error)

Table 2 Continued

Energy (keV) Temperature (K) Number of cascades Surviving MD Clustered interstitials (mean /

displacements standard deviation / standard error)

(mean / standard deviation / standard error)

The configuration, thermal stability and mobility of vacancy, and SIA clusters are of importance for the kinetics of damage accumulation and are different in the fcc and bcc metals. In the fcc metals, vacancy clusters are in the form of either dislocation loops or SFTs, depending on the stacking-fault energy, and the fraction of clustered vacancies, ev, is close to that for the SIAs, ei. In the bcc metals, nascent vacancy clusters usually form loosely correlated 3D configurations, and ev is much smaller than ei. Gen­erally, vacancy clusters are considered to be immo­bile and thermally unstable above the temperature corresponding to the recovery stage V.

In contrast to vacancy clusters, the SIA clusters are mainly in the form of a 2D bundle of crowdions or small dislocation loops. They are thermally stable and highly mobile, migrating 1D in the close-packed crystallographic directions.45 The ability of SIA clus­ters to move 1D before being trapped or absorbed by a dislocation, void, etc. leads to entirely different reac­tion kinetics as compared with that for 3D migrating defects, and hence may result in a qualitatively differ­ent damage accumulation than that in the framework of the FP3DM (see Section 1.13.6).

It should be noted that MD simulations provide maximum evidence for the high mobility of small SIA clusters. Numerous experimental data, which also sup­port this statement, are discussed in this chapter, how­ever, indirectly. One such fact is that most of the loops formed during ion irradiations of a thin metallic foil have Burgers vectors lying in the plane of the foil.54 It should also be noted that recent in situ experiments55-58 provide interesting information on the behavior of interstitial loops (>1 nm diameter, that is, large enough to be observable by transmission electron microscope, TEM). The loops exhibit relatively low mobility, which is strongly influenced by the purity of materials. This is not in contradiction with the simulation data. The observed loops have a large cross-section for interaction with impurity atoms, other crystal imperfections and other loops: all such interactions would slow down or even immobilize interstitial loops. Small SIA clusters produced in cascades consist typically of approximately ten SIAs and have, thus, much smaller cross-sections and consequently a longer mean-free path (MFP). The influence of impurities may, however, be strong on both the mobility of SIA clusters and, consequently, void swelling is yet to be included in the theory.

The emission or dissociation rate is usually the sum of the binding energy of the emitted particle and its migration energy. As in the case of migration energy, the binding energies can be obtained using either experimental studies, ab initio calculations, or MD.

As stated previously, three kinds of KMC techni­ques (AKMC, OKMC, and EKMC) have been used so far to model microstructural evolutions during radia­tion damage. In atomistic KMC, the evolution of a complex microstructure is modeled at the atomic scale, taking into account elementary atomic mech­anisms. In the case of diffusion, the elementary mech­anisms leading to possible state changes are the diffusive jumps of mobile point defect species, includ­ing point defect clusters. Typically, vacancies and SIAs can jump from one lattice site to another lattice site (in general first nearest neighbor sites). If foreign interstitial atoms such as C atoms or He atoms are included in the model as in Hin et a/.,70,71 they lie on an interstitial sublattice and jump on this sublattice.

In OKMC, the microstructure consists of objects which are the intrinsic defects (vacancies, SIAs, dis­locations, grain boundaries) and their clusters (‘pure’ clusters, such as voids, SIA clusters, He or C clus­ters), as well as mixed clusters such as clusters con­taining both He atoms, solute/impurity atoms, and
interstitials, or vacancies. These objects are located at known (and traced) positions in a simulation volume on a lattice as in LAKIMOCA or a known spatial position as in BIGMAC and migrate according to their migration barriers.

In the EKMC approach,72,73 the microstructure also consists of objects. The crystal lattice is ignored and objects’ coordinates can change continuously. The only events considered are those which lead to a change in the defect population, namely clustering of objects, emission of mobile species, elimination of objects on fixed sinks (surface, dislocation), or the recombination between vacancy and interstitial defect species. The migration of an object in its own right is considered an event only if it ends up with a reaction that changes the defect population. In this case, the migration step and the reaction are processed as a single event; otherwise, the migration is performed only once at the end of the EKMC time interval At. In contrast to the RTA, in which all rates are lumped into one total rate to obtain the time increment, in an EKMC scheme the time delays of all possible events are calculated separately and sorted by increasing order in a list. The event corresponding to the shortest delay, ts, is processed first, and the remaining list of delay times for other events is modified accordingly by eliminating the delay time associated with the particle that just disappeared, adding delay times for a new mobile object, etc.

To illustrate the power of KMC for modeling radi­ation effects in structural materials and nuclear fuels, this chapter next considers two examples, namely the use of AKMC simulations to predict the coupled evo­lution of vacancy clusters and copper precipitates dur­ing low dose rate neutron irradiation of Fe-Cu alloys and the use of an OKMC model to predict the trans­port and diffusional release of fission product, silver, in tri-isotropic (TRISO) nuclear fuel. These two examples will provide more details about the possible implementations of AKMC and OKMC models.

Another way to view solid solutions and nonstoichio­metry is as a function of defects in the ideal lattice.

This has been of particular interest for oxide fuels as they are seen to govern dissolution of cations and nonstoichiometry in oxygen behavior and as a result, transport properties. Defect concentrations are inherent in the CEF, as vacancies and interstitials on the oxygen lattice for fluorite-structure actinide systems are treated as constituents linked to cations (see Section 1.17.4.3). A more explicit treatment of oxide systems with point defects has been applied to a wide range of materials such as high-temperature oxide superconductors, TiO2, and ionic conducting membranes, among others. For oxide fuels, point defects have been described thermochemically by a number of investigators starting as early as 1965 with more recent treatments in fuels by Nakamura and Fujino,36 Stan eta/.,37 and Nerikar eta/.38 Oxygen site defects, which dominate in the fluorite-structure fuels, are of course driven by the multiple possible valence states of the actinides, most notably uranium, which can exhibit U+3, U+4, U+5, and U+6.

A simple example of the point defect treatment can be seen in Stan et a/.37 They optimized defect concentrations from the defect reactions described in the Kroger-Vinck notation

oo = + Vo

VO2 + 2U U = Ot + 2U’V

A dilute defect concentration was assumed such that there were no interactions between defects and thus no excess energy terms. The results of the fit to literature data are seen in Figure 7(a), where the stoichiometry of the fluorite-structure hyperstoi­chiometric urania is plotted as a function of defect concentration xa. The relationships were also used to compute oxygen potentials as a function of stoichi­ometry and are plotted in Figure 7(b) illustrating relatively good agreement with values computed by Nakamura and Fujino.36

AKMC simulations can be used to follow the atomic configuration of a finite-sized system, starting from a given initial condition, by performing successive

Figure 9 Comparison of Cr segregation profiles as a function of irradiation dose in FeNi12Cr19 at T = 635 K. Left cell represents typical experimental results of Busby et al. ,36 right cell is mean-field predicted result starting from the experimental profile observed just before irradiation. Reproduced from Nastar, M. Philos. Mag. 2005, 85, 641-647.

point defect jumps.16’17’126’127 Then, the migration barriers are exactly computed (in the framework of the used diffusion model) for each atomic configura­tion using equations such as [22] or [23], without any mean-field averaging. The jump to be performed can be chosen with a residence time algorithm,128,129 which can also easily integrate creation and annihila­tion events.16

Correlation effects between successive point defect jumps, as well as thermal fluctuations, are naturally taken into account in AKMC simulations; this pro­vides a good description of diffusion properties and of nucleation events (the latter being especially impor­tant for the modeling of RIP). The downside is that they are more time consuming than mean-field models, especially when correlation or trapping effects are significant. These advantages and drawbacks explain why AKMC is especially useful to model the first stages of segregation and precipitation kinetics.

AKMC simulations were first used to study RIS and RIP in model systems,16,17 to highlight the con­trol of segregation by point defect migration mecha­nism, and to test the results of classical diffusion equations. These studies show that it is possible — in favorable situations, where correlation and trapping effects are not too strong — to simulate microstruc­ture evolution with realistic dose rates and point defect concentrations, up to doses of typically 1 dpa.

In simple cases, AKMC simulations can validate the predictions of continuous models, on the basis of simple diffusion equations17,16: an example is given in Figure 10 for an ideal solid solution, that is, when RIS can be studied without any clustering or ordering tendency.16 In this simulation, the diffusion of

A atoms by vacancy jumps is more rapid than that of B atoms, and one observes an enrichment of B atoms at the sinks due to the IK effect (A and B atoms diffuse at the same speed by interstitial jumps; those jumps therefore do not contribute to the segre­gation). One can notice the nonmonotonous shape of the concentration profile, which corresponds to the prediction of Okamato and Rehn29 when the partial diffusion coefficients are dBV/dAV < dBI/du. In the more complicated case of a regular solution, Rottler et a/.17 have shown that the RIS profiles of AKMC simulations can be reproduced by a continu­ous model using a self-consistent formulation, which gives the dependence of the partial diffusion coeffi­cients with the local composition.

In alloys with a clustering tendency, AKMC simu­lations have been used to study microstructure evo­lution when RIS and precipitation interact, either in under — or supersaturated alloys. The evolution of the precipitate distribution can be quite complicated as the kinetics of nucleation, growth, and coarsening depend on both the local solute concentrations and the point defect concentrations (which control solute diffusion), concentrations that evolve abruptly in the vicinity of the sinks.16 A case of homogeneous RIP, due to a mechanism similar to the one proposed by Cauvin and Martin52 (see Section 1.18.2.4), has been simulated with a simplified interstitial diffusion model.130131

The application of AKMC to real alloys has just been introduced. Copper segregation and precipitation in a-iron has been especially studied, because of its role in the hardening of nuclear reactor pressure vessel steels.118’126’127’132’133 These studies are based on rigid

10-4

10-5

10-6

10-7

10-8

10-9

10-10

10-11

10-12

10-13

Figure 10 (a) Evolution of the B concentration profile in an A10B90 ideal solid solution under irradiation at 500 K and 10—3dpas—1 when dAV > dBV and dAi = dBi; (b) Point defect concentration profiles in the steady state. Reproduced from Soisson, F. J. Nucl. Mater. 2006, 349, 235-250.

lattice approximations, using parameters fitted to DFT calculations. The point defect formation energies are found to be much smaller in copper-rich coherent clusters than in the iron matrix,79,118 and there is a strong attraction between vacancies and copper atoms in iron, up to the second-nearest neighbor posi — tions.70,125 AKMC simulations show that in dilute Fe-Cu alloys, the Lcuv = — [Lcucu + LcuFe] is positive at low temperatures, because of the diffusion of Cu-V pairs. At higher temperatures, Cu-v pairs dissociate and LcuV becomes negative.118,119 The resulting segre­gation behaviors have been simulated, with homoge­neous formation of Frenkel pairs (i. e., conditions corresponding to electron irradiations). Only vacancy fluxes are found to contribute to RIS; copper concen­tration increases near the sinks at low temperatures and decreases at high temperatures.118 In highly supersatu­rated alloys, RIS does not significantly modify the evo­lution of the precipitate distribution, except for the acceleration proportional to the point defect supersat­uration. Figure 11 illustrates a simulation of RIS in a very dilute Fe-Cu alloy, performed with the para­meters of Soisson and Fu118,125; the segregation of copper at low temperatures produces the preferen­tial formation of small copper-rich clusters near the sinks, which could correspond to the beginning of a
heterogeneous precipitation. However, simulations are limited to very short irradiation doses because of the trapping of defects as soon as the first clusters are formed. This makes it difficult to draw conclusions from these studies.

Coprecipitation of copper clusters with other solutes (Mn, Ni, and Si) has been modeled by Vincent et a/.126,127 and Wirth and Odette133 under irradiation at very high radiation fluxes and with formation of point defects in displacement cascades. AKMC simulations display the formation of vacancy clusters surrounded by copper atoms, which could result both from the Cu-V attraction (a purely thermodynamic factor) and from the dragging of Cu by vacancies (effect of kinetic coupling). The forma­tion of Mn-rich clusters is favored by the positive coupling between fluxes of self-interstitials and Mn (DFT calculations show that the formation of mixed Fe-Mn dumbbells is energetically favored).

1.18.4 Conclusion

We started this review with a summary of the exper­imental activity on RIS. Intensive experimental work has been devoted to austenitic steels and its model fcc alloys (Ni-Si, Ni-Cr, and Ni-Fe-Cr) and, more recently, to ferritic steels. A strong variation of RIS with irradiation flux and dose, temperature, composition, and the grain boundary type was observed. One study tried to take advantage of the sensitivity of RIS to composition to inhibit Cr depletion at grain bound­aries. A small addition of large-sized impurities, such as Zr and Hf, was shown to inhibit RIS up to a few dpa in both austenitic and ferritic steels. On the other hand, the Cr depletion at grain boundaries was observed to be delayed when the grain boundary was enriched in Cr before irradiation. A ‘W-shaped’ transitory profile could be maintained until a few dpa before the grain boundary was depleted in Cr. The mechanisms involved in these recent experiments are still not well understood, although RIS model development was closely related to the experimental study.

The main RIS mechanisms had been understood even before RIS was observed. From the first models, diffusion enhancement and point defect driving forces were accounted for. The kinetic equations are based on general Fick’s laws. While in dilute alloys one knows how to deduce such equations from atomic jump frequencies, in concentrated alloys more empir­ical methods are used. In particular, the definition of the partial diffusion coefficients of the Fick’s laws in terms of the phenomenological L-coefficients of TIP has been lost over the years. This can be explained by the lack of diffusion data and diffusion theory to deter­mine the L-coefficients from atomic jump frequencies. For years, the available diffusion data consisted mainly of tracer diffusion coefficients, and the RIS models employed empirical Manning relations, which approximated partial diffusion coefficients based on tracer diffusion coefficients. However, recent improve­ments of the mean-field diffusion theories, including short-range order effects for both vacancy and intersti­tial diffusion mechanism, are such that we can expect the development of more rigorous RIS models for concentrated alloys. It now seems possible to overcome the artificial dichotomy between dilute and concen­trated RIS models and develop a unified RIS model with, for example, the prediction of the whole kinetic coupling induced by an impurity addition in a con­centrated alloy. Meanwhile, first-principles methods relying on the DFT have improved so fast in the last decades that they are able to provide us with activation energies of both vacancy and interstitial jump frequencies as a function of local environment. Therefore, it now seems easier to calculate the L — coefficients and their associated partial diffusion coef­ficients from first-principles calculations rather than estimating them from diffusion experiments.

An alternative approach to continuous diffusion equations is the development of atomistic-scale simu­lations, such as mean-field equations or Monte Carlo simulations, which are quite appropriate to study the nanoscale RIS phenomenon. Although the mean-field approach did not reproduce the whole flux coupling due to the neglect ofcorrelation effects, it predicted the main trends of RIS in austenitic steels, with respect to temperature and composition, and was useful to under­stand the interplay between thermodynamics and kinetics during the formation of an oscillating transi­tory profile. Monte Carlo simulations are now able to embrace the full complexity of RIS phenomena, including vacancy and split interstitial diffusion mechanisms, the whole flux coupling, the resulting segregation, and eventual nucleation at grain bound­aries. But these simulations become heavy time­consuming when correlation effects are important.

roS

1 + oS + ko(1 — S) 3

_______ rp3________

1 + oS + ko(1 — S) r(1 + oS)

1 + oS + ko(1 — S) 3r(1 + oS)

1 + oS + ko(1 — S) 3r(1 + oS)

1 + oS + ko(1 — S) r(1 + o)S

[2] + oS + ko(1 — S)

ro(1 — S)S

1 + oS + ko(1 — S)

6r2u m (1 — S)

[6] Chiang, Y. M.; Birnie, D.; Kingery, W. D. Physical Ceramics: Principles for Ceramic Science and Engineering; MIT Press: Cambridge, 1997.

[7] Kittel, C. Introduction to Solid State Physics; Wiley:

New York, 1996.

[8] keV, the exponent is 0.75. This is marginally lower than the value in eqn [2], possibly because the 20 keV data were used in the current fitting. An exponent of 1.03 was found in the range above

20 keV, which is dominated by subcascade formation. Only in the highest energy range do the MD results

approach the linear energy dependence predicted by the NRT model. The range of plus or minus

one standard error is barely detectable around the

data points, indicating that the change in slope is statistically significant.

The data from Figure 9(a) are replotted in

Figure 9(b) where the number of surviving displace­ments is divided by the NRT displacements at each energy. The rapid decrease in this MD defect survival

Dvexp(—Evcl/kB T) (k( + Zdpd) — £gGv

[16]

[17]

where Evcl is an effective binding energy of vacancies with the vacancy clusters and x[[d м are the mean size of the vacancy and SIA glissile clusters (see eqn [8]).

It should be noted that the SDF of sessile intersti­tial clusters, the sink strength of which is described by eqn [132], is limited by the maximum size of the clusters produced in displacement cascades (see Figure 2 in Singh eta/.22). This is because the clusters

[20] mm

Figure 12 Three-dimensional view for the formation of

dislocation channels on glide planes which have low

stacking fault tetrahedron (SFT) density, for irradiated

copper at a dose of 0.01 dpa. For clarity of visualization, the

apparent SFT density has been reduced by a factor of 100.

Note that other dislocation segments are inactive as a result

of the high density of surrounding clusters.

As differences in dose rates can confound direct comparison between neutron and ion irradiations, it is important to assess their impact. A simple method for examining the tradeoff between dose and temper­ature in comparing irradiation effects from different particle types is found in the invariance requirements. For a given change in dose rate, we would like to know what change in dose (at the same temperature) is required to cause the same number of defects to be absorbed at sinks. Alternatively, for a given change in dose rate, we would like to know what change in temperature (at the same dose) is required to cause the same number of defects to be absorbed at sinks. The number of defects per unit volume, NR, that have recombined up to time t, is given by Mansur14

t
Ci Cv dt	[12]
Nr
0
where kSj is the strength of sink j and Cj is the sink concentration. The ratio of vacancy loss to interstitial loss is
Nsv Nsi
[14]
Rs
where j = v or i. The quantity Ns is important in describing the microstructural development involving total point defect flux to sinks (e. g., RIS), while Rs is the relevant quantity for the growth of defect aggregates such as voids that require partitioning of point defects to allow growth. In the steady-state recombination dominant regime, for Ns to be invariant at a fixed dose, the follow­ing relationship between ‘dose rate (K) and temperature (Ti)’ must hold:
Figure 9 shows plots of the relationship between the ratio of dose rates and the temperature difference required to maintain the same point defect absorption at sinks (a), and the swelling invariance (b). The invariance requirements can be used to prescribe an ion irradiation temperature-dose rate combination that simulates neutron radiation. We take the example of irradiation of stainless steel under typical BWR core irradiation conditions of ^4.5 x 10-8 dpa s-1 at 288 °C. If we were to conduct a proton irradiation with a characteristic dose rate of 7.0 x 10-6dpas-1, then using eqn [15] with a vacancy formation energy of 1.9 eV and a vacancy migration
£M t
[15]
where Evm is the vacancy migration energy. In the steady-state recombination dominant regime, for Rs to be
Figure 9 Temperature shift from the reference 200 °C required at constant dose in order to maintain (a) the same point defect absorption at sinks, and (b) swelling invariance, as a function of dose rate, normalized to initial dose rate. Results are shown for three different vacancy migration energies and a vacancy formation energy of 1.5 eV. Adapted from Mansur, L. K. J. Nucl. Mater. 1993, 206, 306-323; Was, G. S. Radiation Materials Science: Metals and Alloys; Springer: Berlin, 2007.

energy of 1.3 eV, the experiment will be invariant in Ns with the BWR core irradiation (e. g., RIS) at a proton irradiation temperature of 400 °C. Similarly, using eqn [16], a proton irradiation temperature of 300 °C will result in an invariant Rs (e. g., swelling or loop growth). For a Ni2+ ion irradiation at a dose rate of 10~3 dpas, the respective temperatures are 675 °C (Ns invariant) and 340 °C (Rs invariant). In other words, the temperature ‘shift’ due to the higher dose rate is dependent on the microstructure feature of interest. Also, with increasing difference in dose rate, the AT between neutron and ion irradiation increases substantially. The nominal irradiation tem­peratures selected for proton irradiation, 360 °C and for Ni2+ irradiation, 500 °C represent compromises between the extremes for invariant Ns and Rs.