As shown in Section 1.13.5, in the framework of FP3DM, the swelling rate depends on the dislocation density and becomes small for a low dislocation den­sity, dS/df « Bdpd/k2 ! 0 at pd! 0 (see eqn [96]). Thus, it was a common belief that the swelling rate in well-annealed metals has to be low at small doses, that is, when the dislocation density increase can be neglected. Under neutron irradiation, the effect of dislocation bias on swelling is even smaller because of intracascade recombination: (dS/df)^^ =

(dS/dfED(1 — er) < (dS/df)Zf. It has been

found experimentally, however, that the void swelling rate in fully annealed pure copper irradiated with fission neutrons up to about 10—2dpa (see Singh and Foreman18) is of ~1% per dpa, which is similar to the maximum swelling rate found in materials at high enough irradiation doses. This observation was one of those that prompted the development of the PBM. The production bias term in eqn [138] allows the understanding of these observations. Indeed, at low doses of irradiation, the void size is small, and therefore, the void cross-section for the inter­action with the SIA glissile clusters is small (Krc2Nc/Lg ^ 1). As a result, the last term in the pro­duction bias term is negligible and thus the swelling rate is driven by the production bias:

 

dS df

 

[138] where f = GNRTt is the NRT irradiation dose. The first term in the brackets on the right-hand side of eqn [138] corresponds to the influence of the dislocation bias and the second one to the production bias. The factor (1 — er) describes intracascade recombination of defects, which is a function of the recoil energy and may reduce the rate of defect production by up to an order of magnitude that can be compared to the NRT value: (1 — er)! 0.1 at high PKA energy (see Section 1.13.3). As indicated by this equation, the swelling rate is a complicated function of dislocation density, dislocation bias factor, and the densities and sizes of voids and PD clusters. It also demonstrates the dependence of the swelling rate on the recoil energy, determined by eig, which increases with increasing PKA energy up to about 10-20 keV. The main predictions of the PBM are dis­cussed below.

 

1.13.6.2 Main Predictions of Production Bias Model As can be seen from eqn [138], the action and con­sequences of the two biases, the dislocation and production ones, is quite different. As shown in Section 1.13.5, the dislocation bias depends only slightly on the microstructure and predicts indefinite void growth. In contrast, the production bias can be positive or negative, depending on the microstruc­ture. The reason for this is in negative terms in eqn [138]. The first term decreases the action of the

sessile vacancy and SIA clusters, the swelling rate is given by dS/df ~ 1/2(1 — er)eg where the sink strength ratio, k//(k2 + Z/pd), is taken to be equal to 1/2, as frequently observed in experiments. Taking into account the magnitude of the cascade para­meters er and eg estimated in Golubov et a/.24 and neglecting the dislocation bias term in eqn [138], one may conclude that the maximum swelling rate under fast neutron irradiation may reach about 1% per dpa. As pointed out in Section 1.13.5, in the case of FP production, that is, in the FP3DM, the maximum swelling rate is also ~1% per dpa. This coincidence is one of the reasons why an illusion that the FP3DM model is capable of describing damage accumula­tion in structural and fuel materials in fission and future fusion reactors has survived despite the fact that nearly 20 years have passed since the PBM was introduced.

Note that the production bias provides a way to understand another experimental observation, namely, that the swelling rate in some materials decreases with increasing irradiation dose (see, e. g., Figure 5 in Golubov et a/.24). Such a decrease is predicted by eqn [138], as the negative term of the production bias, %r Nc/Lg, increases with an increase in the void size. As the first term in the 10-4

Figure 5 Experimentally measured133 and calculated24 levels of void swelling in pure copper after irradiation with 2.5 MeV electrons, 3MeV protons, and fission neutrons. The calculations were performed in the framework of the FP3DM for the electron irradiation and using the production bias model as formulated in Singh etal.22 for irradiations with protons and fission neutrons. From Golubov etal.24

production bias is proportional to the void radius and the second to the radius squared, the swelling rate may finally achieve saturation at a mean void radius equal to Rmax ~ 2яг/.19,30,35

Finally, the cascade production of the SIA clusters may strongly affect damage accumulation. As can be seen from eqn [132], the steady-state sink strength of the sessile SIA clusters is inversely proportional to the fraction of SIAs produced in cascades in the form of mobile SIA clusters, thus k2d ! 1 when eg! 0. This limiting case was considered by Singh and Foreman18 to test the validity of the original frame­work of the PBM,16,17 where all the SIA clusters produced by cascades were assumed to be immobile (hereafter this case of eg = 0 is called the Singh— Foreman catastrophe). If for some reasons this case is realized, void swelling and the damage accumulation in general would be suppressed for the density of SIA clusters, hence, their sink strength would reach a very high value by a relatively low irradiation dose, f ^ 1dpa, (see Singh and Foreman18). Thus, irradiation with high-energy particles, such as fast neutrons, provides a mechanism for suppressing damage accumulation, which may operate if the SIA clusters are immobilized. In alloys, the interaction with impurity atoms may provide such an immobili­zation. The so-called ‘incubation period’ of swelling observed in stainless steels under neutron irradia­tion for up to several tens of dpa (Garner32,33) might be due to the Singh—Foreman catastrophe. A possible scenario of this may be as follows: during the incubation period, the material is purified by RIS mainly on SIA clusters because of their high density. At high enough doses, that is, after the incu­bation period, the material becomes clean enough to provide the recovery of the mobility of small SIA clusters created in cascades that triggered on the production bias mechanism. As a result, the high number density of SIA clusters decreases via the absorption of the excess of vacancies, restoring con­ditions for damage accumulation.

From numerous careful experimental studies con­ducted on the BDT behavior of steels, it is now established that precracks (macrocracks) blunt sub­stantially before the fracture of the specimen occurs at the transition region. However, the examination of the fracture surface reveals that cracks propagate predominantly by cleavage.44 Several cracked brit­tle particles are found to be present in the broken samples,45,53 and the measured microscopic fracture stress (at the microcracks) is found to be a few orders of magnitude higher than that of the pure Griffith value.53,54 All these observations are considered in our model as follows:

1. We implemented the blunting of macrocracks by using the elastic crack-tip stress field for blunted cracks. As dislocations are emitted, the crack blunts and the radius of curvature increases. The notional crack tip, which is taken as reference for calculation ofimage stresses, retreats away from the blunted tip.

2. A microcrack is placed in the field of a macrocrack and the failure criterion used in the calculation of the cleavage crack propagation from this microcrack.

3. We consider the emission of dislocations and subsequent shielding from the microcrack tip (a detailed study of this and the observed constancy for microscopic fracture stress is reported in an earlier study ).

The geometry of the model used for simulation is shown in Figure 13. A semi-infinite crack (macro­crack) with a finite microcrack situated ahead of it on its crack plane is loaded starting from K = 0. Disloca­tion sources are assumed to exist at a distance x0 from the tip, and are situated on the slip planes passing through the crack tip. During loading, dislocations are emitted from source positions (x0) when the resolved shear stress reaches a value of 2 т, The resolved shear stresses are obtained using expressions based on deri­vations for a semi-infinite crack58 and a finite crack59 for the respective cases. The emitted dislocations move along the slip plane away from the crack tip, and the stress at the source increases until another dislocation is emitted. The emitted dislocations move with a velocity based on the following expression:

Vx,. = jj-^ (Itx,.)mAe(-Ea’iT> [24]

V

У

у

The values for the parameters were obtained by fit­ting the data of screw dislocation velocities in iron.59 The value of m has a linear dependence on tempera­ture T: m = 400/T + 1.2, A = 3.14 x 10-4, and Ea = 0.316 eV. The first term restricts the motion of dis­locations below the friction stress value (t,), making sure that v = 0 for |tx, |< t,, and hence, most of the dislocations are in near-equilibrium positions at any given time. When the dislocations are in their equi­librium positions, the temperature and strain-rate dependence of measured fracture toughness (KF), plastic zone size (df), crack-tip opening displacement, etc. are determined only by the temperature and strain-rate dependence of the friction stress (t,). The friction stress used is chosen to be equal to the shear yield stress ay/2 when the Tresca yield criterion is assumed. Thus, the temperature dependence of frac­ture toughness is obtained by inputting the corresponding friction stress value for each tempera­ture. Simulations were done for temperatures range from -180 to -60 °C with corresponding yield stress values (ot,,) from 910 to 620 MPa.

The arrays of emitted dislocations form the plastic zones of the crack. The crack may also get blunted due to dislocation emission. In our case, the effects of blunting will be negligible for microcracks since the number of dislocations emitted is only up to 102. However, the effects of blunting will be significant in the case of macrocracks, because the number of dislocations emitted is of the order 105; here the blunting effects are accounted for by using the elastic crack-tip field for a blunted crack.60 The plastic zone developed at the macrocrack modifies the field ahead

of it so that it is the same as an elastic-plastic material with hardening.61 The microcrack placed in this field experiences a tensile stress and is assumed to propa­gate, leading to fracture when it reaches a critical value F (estimated on the basis of similar dislocation simulation of finite crack emitting dislocations along the slip planes). Computer simulations are performed in two stages. First, the microcrack is loaded to failure and the microscopic fracture stress is estimated for specific crystallographic orientations and crack sizes. The obtained microscopic fracture stress (sf) values are then used as the fracture criterion in the macro­crack simulation.

Given an atomic jump frequency model, transitions of the alloy configurations are described by a master equation. With one point defect in the system, a Monte Carlo simulation produces a trajectory of the alloy in the configuration space. The L-coefficients can be obtained from a Monte Carlo simulation. Measurements are performed on an equilibrated system using the generalized Einstein formula of Allnatt.82 This numerical approach has proved its efficiency; however, the achievement of a predictive model using this method is limited to short ranges of composition and temperature. Simulations become rapidly unworkable when, for example, binding ener­gies between interstitials and neighboring atoms are significant.9,83 The simulated trajectory can be trapped in some configurations due to the correlation effects of the diffusion mechanism; after a jump, an atom has a finite probability to exchange again with the same point defect and cancel its first jump. The escape probability of the point defect from an atom decreases with the binding energy between the two species.

In the limited case of dilute alloys in which a few point defect jump frequencies are involved, it is pos­sible to consider all the vacancy paths and deduce analytical expressions for the L-coefficients. On the other hand, diffusion models of concentrated alloys lead to approximate expressions of the transport coefficients.

The calculation time of ab initio calculations varies — to first order — as the cube of the number of atoms or equivalently of electrons (the famous N3 dependence) in the cell. If a fine к-point sampling is needed, this dependence is reduced to N2 as the number of к points decreases in inverse proportion with the size of the cell. On the other hand, the number of self­consistent cycles needed to reach convergence tends to increase with N. Anyway, the variation of calcula­tion time with the size of the cell is huge and thus strongly limits the number of atoms and also the cell size that can be considered. On one hand, calculations on the unit cell of simple crystalline materials (with a small number of atoms per unit cell) are fast and can easily be performed on a common laptop. On the other hand, when larger simulation cells are needed, the calculations quickly become more demanding. The present upper limit in the number of atoms that can be considered is of the order of a few hundreds. The exact limit of course depends on the code and also on the number of electrons per atoms and other technicalities (number of basis functions, к points, available computer power, etc.), so it is not possible to state it precisely. Considering such large cells leads anyway to very heavy calculations in which the use of parallel versions of the codes is almost mandatory. Various parallelization schemes are possible: on k points, fast Fourier transform, bands, spins; the parallelization schemes actually available depend on the code.

The situation gets even worse when one notes that a relaxation roughly involves at least ten ground-state calculations, a saddle point calculation needs about ten complete relaxations, and that each molecular dynamics simulation time step (of about 1 fs) needs a complete ground-state calculation. Overall, one can understand that the CPU time needed to complete an ab initio study (which most of the time involves vari­ous starting geometry) may amount up to hundreds of thousands or millions of CPU hours.

Dislocations are line defects in crystals, and their motion is the carrier for plastic deformation of
crystals under most conditions (T< Tn/2).36,37 The defects produced by irradiation (such as vacancy and interstitial complexes) interact with dislocations, and this interaction is responsible for the change in the mechanical properties by irradiation (such as embrittlement).38 MD simulations of dislocation interaction with other defects are discussed in detail in Chapter 1.12, Atomic-Level Level Dislocation Dynamics in Irradiated Metals. Here, we describe a more basic case study on the mobility of an edge dislocation in Ta. In Section 1.09.6.2.1, we describe the method of computing its Peierls stress, which is the critical stress to move the dislocation at zero temperature. In Section 1.09.6.2.2, we describe how to compute the mobility of this dislocation at finite temperature by MD.

The evolution of displacement cascades is similar at all energies, with the development of a highly ener­getic, disordered core region during the initial, colli — sional phase of the cascade. Vacancies and interstitials are created in equal numbers, and the number of point defects increases sharply until a peak value is reached. Depending on the cascade energy, this occurs at a time in the range of 0.1-1 ps. This evolu­tion is illustrated in Figure 5 for a range of cascade energies, where the number of vacancies is shown as a function of the cascade time. Many vacancy — interstitial pairs are in quite close proximity at the time of peak disorder. An essentially athermal process of in-cascade recombination of these close pairs takes place as they lose their kinetic energy. This leads to a reduction in the number of defects until a quasi­steady-state value is reached after about 5-10 ps. As interstitials in iron are mobile even at 100 K, further short-term recombination occurs between some vacancy-interstitial pairs that were initially separated by only a few atomic jump distances. Finally, a stage is reached where the remaining point defects are sufficiently well separated that further recombina­tion is unlikely on the time scale (a few hundred picoseconds) accessible by MD. Note that the number of stable Frenkel pair is actually somewhat lower than the value shown in Figure 5 because the values obtained using the effective sphere identification
procedure were not corrected to account for the interstitial structure discussed above.

A mechanism known as RCS may help explain some aspects of cascade structure.24,41 An RCS can be visualized as an extended defect along a close — packed row of atoms. When the first atom is pushed off its site, it dissipates some energy and pushes a second atom into a third, and so on. When the last atom in this chain is unable to displace another, it is left in an interstitial site with the original vacancy several atomic jumps away. Thus, RCSs provide a mechanism of mass transport that can efficiently sep­arate vacancies from interstitials. The explanation is consistent with the observed tendency for the final cascade state to be characterized by a vacancy-rich central region that is surrounded by a region rich in interstitial-type defects. However, although RCSs are observed, particularly in low-energy cascades, they do not appear to be prominent enough to explain the defect separation observed in higher energy cascades.58 Visualization of cascade dynamics indicates that the separation occurs by a more collective motion of multiple atoms, and recent work by Calder and coworkers has identified a shockwave-induced mech­anism that leads to the formation of large interstitial clusters at the cascade periphery.80 This mechanism will be discussed further in Section 1.11.4.3.1. Coherent displacement events involving many atoms have also been reported by Nordlund and coworkers.81

Defect production tends to be dominated by a series of simple binary collisions at low PKA ener­gies, while the more collective, cascade-like behav­ior dominates at higher energies. The structure of typical 1 and 20 keV cascades is shown in Figure 6, where parts (a) and (b) show the peak damage state and (c) and (d) show the final defect configurations. The MD cells contained 54 000 and 432 000 atoms for the 1 and 20 keV simulations, respectively. Only the vacant lattice sites and interstitial atoms identi­fied by the effective sphere approach described above are shown. The separation of vacancies from interstitials can be seen in the final defect config­urations; it is more obvious in the 1 keV cascade because there are fewer defects present. In addition to isolated point defects, small interstitial clusters are also clearly observed in the 20 keV cascade debris in Figure 6(d). In-cascade clustering is dis­cussed further in Section 1.11.4.3.

The morphology of the 20 keV cascade in Figure 6(b) exhibits several lobes which are evidence of a phenomenon known as subcascade formation.82 At low energies, the PKA energy tends to be dis­sipated in a small volume and the cascades appear as compact, sphere-like entities as illustrated by the
1 keV cascade in Figure 6(a). However, at higher energies, some channeling82,83 of recoil atoms may occur. This is a result of the atom being scattered into a relatively open lattice direction, which may permit it to travel some distance while losing rela­tively little energy in low-angle scattering events. The channeling is typically terminated in a high — angle collision in which a significant fraction of the recoil atom’s energy is transmitted to the next generation knock-on atom. When significant sub­cascade formation occurs, the region between high-angle collisions can be relatively defect-free as the cascade develops. This evolution is clearly shown in Figure 7 for a 40 keV cascade, where the branching due to high-angle collisions is observed on a time scale of a few hundreds of femto seconds. One practical implication of subcascade formation is that very high-energy cascades break up into what looks like a group of lower energy cascades. An example of sub­cascade formation in a 100 keV cascade is shown in Figure 8 where the results of 5 and 10 keV cascades have been superimposed into the same block of atoms for comparison. The impact of subcascade formation on stable defect production will be discussed in the next section.

Figure 6 Structure of typical 1 keV (a, c) and 20 keV (b, d) cascades. Peak damage state is shown in (a and b) and the final stable defect configuration is shown in (c and d).

MD cascade simulations in iron at 100 K: peak damage

10keV

100 keV

5 keV

5keV-0.26ps 10keV-0.63ps 100keV-0.70 ps

Figure 8 Energy dependence of subcascade formation.

Crystal microstructure under irradiation consists of two qualitatively different defect types: mobile (single vacancies, SIAs, and SIA and vacancy clus­ters) and immobile (voids, SIA loops, dislocations, etc.). The concentration of mobile defects is very small (^ 10 10—10 6 per atom), whereas immobile defects may accumulate an unlimited number of PDs, gas atoms, etc. The mathematical description of these defects is, therefore, different. Equations for mobile defects describe their reactions with immo­bile defects and are often called the rate (or balance) equations. The description of immobile defects is more complicated because it must account for nucle- ation, growth, and coarsening processes.

1.13.4.1 Concept of Sink Strength

The mobile defects produced by irradiation are absorbed by immobile defects, such as voids, disloca­tions, dislocation loops, and GBs. Using a MFA, a crystal

can be treated as an absorbing medium. The absorption rate of this medium depends on the type of mobile defect, its concentration and type, and the size and spatial distribution of immobile defects. A parameter called ‘sink strength’ is introduced to describe the reac­tion cross-section and commonly designated as k(, k(, and k? cl (x) for vacancies, SIAs, and SIA clusters of size x (the number of SIAs in a cluster), respectively. The role of the power ‘2’ in these values is to avoid the use of square root for the MFPs of diffusing defects between production until absorption, which are correspond­ingly k-1, k-1, and kcKx). There are a number of publications devoted to the derivation of sink strengths.40,59-61 Here we give a simple but sufficient introduction to this subject.

1.13.4.2 Equations for Mobile Defects

For simplicity, we use the following assumptions:

• The PDs, single vacancies, and SIAs, migrate 3D.

• SIA clusters are glissile and migrate 1D.

• All vacancy clusters, including divacancies, are immobile.

• The reactions between mobile PDs and clusters are negligible.

• Immobile defects are distributed randomly over the volume.

Then, the balance equations for concentrations of

mobile vacancies, Cv, SIAs, Ci, and SIA clusters,

Cg. l(x), are as follows

dCv = Gnrt (1 — er)(1 — ev) + G;h

— k? dycy — mRDiCiCv N

Gnrt(1 — er)(1 — ei)-k? DiCi

mRDi Ci Cv

G«l(x)-k.’,(x)DidC‘l.

x 2. 3. ••• xmax [ 12 ]

where G* is the rate of thermal emission of vacancies from all immobile defects (dislocations, GBs, voids, etc.); Dv, Di, and Dicl(x) are the diffusion coefficients of vacancies, single SIAs, and SIA clusters, respec­tively; and mR is the recombination coefficient of PDs. Since the dependence of the cluster diffusivity, Dicl(x), and sink strengths, k? d(x), on size x is
rather weak,45,46 the mean-size approximation for the SIA clusters may be used, where all clusters are assumed to be of the size <xf) . In this case, the set of eqn [12] is reduced to the following single equation

^ =<x?)-1Gnrt(1 — er)eg — k? lAdC? d [13] where eqn [9] is used for the cluster generation rate. To solve eqns [10]—[13], one needs the sink strengths k?, k2, and k2d, the rates of vacancy emission from various immobile defects to calculate G*, and the recombination constant, mR. The reaction kinetics of 3D diffusing PDs is presented in Section 1.13.5, while that of 1D diffusing SIA clusters in Section

1.13.5. In the following section, we consider equa­tions governing the evolution of immobile defects, which together with the equations above describe damage accumulation in solids both under irradiation and during aging.

1.13.4.3 Equations for Immobile Defects

The immobile defects are those that preexist such as dislocations and GBs and those formed during irradi­ation: voids, vacancy — and SIA-type dislocation loops, SFTs, and second phase precipitates. Usually, the defects formed under irradiation nucleate, grow, and coarsen, so that their size changes during irradiation. Hence, the description of their evolution with time, t, should include equations for the size distribution func­tion (SDF), f (X, t), where X is the cluster size.

1.13.4.3.1 Size distribution function

The measured SDF is usually represented as a func­tion of defect size, for example, radius, X = R : f (R, t). In calculations, it is more convenient to use x-space, X = x, where x is the number of defects in a cluster: f (x, t). The radius of a defect, R, is connected with the number of PDs, x, it contains as:

4p 3 _

—R3 = xO 3

pR2 b = xO [14]

for voids and loops, respectively, where O is the atomic volume and b is the loop Burgers vector. Cor­respondingly, the SDFs in R — and x-spaces are related to each other via a simple relationship. Indeed, ifsmall dx and dR correspond to the same cluster group, the number density of this cluster group defined by two functions f (x)dx and f (R)dR must be equal, f(x)dx = f(R)dR, which is just a differential form

of the equality of corresponding integrals for the total number density:

1

f (x)dx

x=2

The relationship between the two functions is, thus,

f (R]=f (x)dR [16]

For voids and dislocation loops

Note the difference in dimensionality: the units of f (x) are atom — (or m- ), while f (R) is in m — atom-1 (or m-4), as can be seen from eqn [15]. Also note that these two functions have quite different shapes, see Figure 1, where the SDF of voids obtained by Stoller et a/.62 by numerical integration of the master equation (ME) (see Sections 1.13.4.3.2 and 1.13.4.4.3) is plotted in both R — and x-spaces.

1.13.4.3.2 Master equation

The kinetic equation for the SDF (or the ME) in the case considered, when the cluster evolution is driven by the absorption of PDs, has the following form

@f @x’ f) = Gs(x) + J(x — 1, t)-J(x, t), x > 2 [18]

where Gs(x) is the rate of generation of the clusters by an external source, for example, by displacement

cascades, and J(x, t) is the flux of the clusters in the size-space (indexes ‘i’ and ‘v’ in eqn [18] are omitted). The flux J (x, t) is given by

J(x, t) = P(x, t)f (x, t) — Q(x + 1, t)f (x + 1, t) [19]

where P(x, t) and Q(x, t) are the rates of absorption and emission of PDs, respectively. The boundary conditions for eqn [18] are as follows

f (1) = C

f (x! 1) = 0 [20]

where C is the concentration of mobile PDs.

If any of the PD clusters are mobile, additional terms have to be added to the right-hand side of eqn [19] to account for their interaction with immobile defect which will involve an increment growth or shrinkage in the size-space by more that unity (see Section 1.13.6 and Singh eta/.22 for details).

The total rates of PD absorption (superscript!) and emission (t—) are given by

11 J! = P(x)f (x), Jto = Q (x )f (x) [21]

x= 2 x= 2

where the superscript arrows denote direction in the size-space. J! and are related to the sink strength

of the clusters, thus providing a link between equa­tions for mobile and immobile defects. For example, when voids with the SDF fc(x) and dislocations are only presented in the crystal and the primary damage is in the form of FPs, the balance equations are

dQ = gNRT — £r)

dt

— [uRACiCv + ZdrdDv(Cv — Cvo)]

— [Pc(1)fc(1, t) — Qvc(2)f (2, t)]

x=1

-E(pc(x)f (x’t)

x=1

— Qvc(x + 1)fc(x + 1, t)) [22]

GNRT (1 — er)

— [mRDiCiCv + ZfpdDiCi]

x=1

-^2 Aic(x + 1)fc(x +1’t) where рА and are the dislocation density and its efficiencies for absorbing PDs, mR, is the recom­bination constant (see Section 1.13.5); the last two terms in eqn [22] describe the absorption and

emission of vacancies by voids and the last term in eqn [23] describes the absorption of SIAs by voids. The balance equations for dislocation loops and sec­ondary phase precipitations can be written in a similar manner. Expressions for the rates P(x, t), Q(x, t), the dislocation capture efficiencies, Zdv, and mR are derived in Section 1.13.5.

Cu is of primary importance in the embrittlement of the neutron-irradiated RPV steels. It has been observed to separate into copper-rich precipitates within the ferrite matrix under irradiation. As its role was discovered more than 40 years ago,74-76 Cu precipitation in a-Fe has been studied extensively under irradiation as well as under thermal aging using atom probe tomography, small angle neutron scattering, and high resolution transmission electron microscopy. Numerical simulation techniques such as rate theory or Monte Carlo methods have also been used to investigate this problem, and we next describe one possible approach to modeling microstructure evolution in these materials.

The approach combines an MD database of pri­mary damage production with two separate KMC simulation techniques that follow the isolated and clustered SIA diffusion away from a cascade, and the subsequent vacancy and solute atom evolution, as discussed in more detail in Odette and Wirth,2 Monasterio eta/.,34 and Wirth eta/.77 Separation of the vacancy and SIA cluster diffusional time scales natu­rally leads to the nearly independent evolution of these two populations, at least for the relatively low dose rates that characterize RPV embrittlement.1,78-81

The relatively short time (~100 ns at 290 °C) evolu­tion of the cascade is modeled using OKMC with the BIGMAC code.77 This model uses the positions of vacancy and SIA defects produced in cascades obtained from an MD database provided by Stoller and coworkers82,83 and allows for additional SIA/ vacancy recombination within the cascade volume and the migration of SIA and SIA clusters away from the cascade to annihilate at system sinks. The duration of this OKMC is too short for significant vacancy migra­tion and hence the SIA/SIA clusters are the only diffus­ing defects. These OKMC simulations, which are described in detail elsewhere,77 thus provide a database of initially ‘aged’ cascades for longer time AKMC cas­cade aging and damage accumulation simulations.

The AKMC model simulates cascade aging and damage evolution in dilute Fe-Cu alloys by following vacancy — nearest neighbor atom exchanges on a bcc lattice, beginning from the spatial vacancy popula­tion produced in ‘aged’ high-energy displacement cascades and obtained from the OKMC to the ulti­mate annihilation of vacancies at the simulation cell boundary, and including the introduction of new cas­cade damage and fluxes of mobile point defects. The potential energy of the local vacancy, Cu-Fe, envi­ronment determines the relative vacancy jump prob­ability to each of the eight possible nearest neighbors in the bcc lattice, following the approach described in eqn [4]. The unrelaxed Fe-Cu vacancy lattice ener­getics are described using Finnis-Sinclair N-body type potentials. The iron and copper potentials are from Finnis and Sinclair84 and Ackland eta/.,85 respec­tively; and the iron-copper potential was developed by fitting the dilute heat of the solution of copper in iron, the copper vacancy binding energy, and the iron-copper [110] interface energy, as described else­where.55 Within a vacancy cluster, each vacancy main­tains its identity as mentioned above, and while vacancy-vacancy exchanges are not allowed, the clus­ter can migrate through the collective motion of its constituent vacancies. The saddle point energy, which is Ea0 in eqn [4], is set to 0.9 eV, which is the activation energy for vacancy exchange in pure iron calculated with the Finnis-Sinclair Fe potential.84

The time (AtAKMC) of each AKMC sweep (or step) is determined by AtAKMC = (nPmax)~ , where Pmax is the highest total probability of the vacancy population and n is an effective attempt frequency. This is slightly different than the RTA, in which an event chosen at random sets the timescale as opposed to always using the largest probability as done here. In this work, n = 1014s_1 to account for the intrinsic vibrational frequency and entropic effects associated with vacancy formation and migration, as used in the previous AKMC model by Odette and Wirth.2 As mentioned, the possible exchange of every vacancy (i) to a nearest neighbor is determined by a Metropolis random number test15 of the relative vacancy jump probability (P/Pmax) during each Monte Carlo sweep. Thus at least one, and often multiple, vacancy jumps occur during each Monte Carlo sweep, which is dif­ferent from the RTA. Finally, as mentioned above, as the total probability associated with a vacancy jump depends on the local environment, the intrinsic time­scale (AtAKMc) changes as a function of the number and spatial distribution of the vacancy population, as well as the spatial arrangement of the Cu atoms in relation to the vacancies.

The AKMC boundary conditions remove (annihi­late) a vacancy upon contact, but incorporate the ability to introduce point defect fluxes through the simulation volume that result from displacement cas­cades in neighboring regions as well as additional displacement cascades within the simulated volume. The algorithms employed in the AKMC model are described in detail in Monasterio et al3 and the remainder of this section will provide highlights of select results.

The AKMC simulations are performed in a ran­domly distributed Fe-0.3% Cu alloy at an irradiation temperature of 290 °C and are started from the spa­tial distribution ofvacancies from an 20 keV displace­ment cascade. The rate of introducing new cascade damage is 1.13 x 10~5 cascades per second, with a cascade vacancy escape probability of 0.60 and a vacancy introduction rate of 1 x 10~4 vacancies per second, which corresponds to a damage rate of this simulation at ^10~ apas~ . Thus a new cascade (with recoil energy from 100 eV to 40 keV) occurs within the simulated volume (a cube of ^86 nm edge length) every 8.8 x 104s (~1 day), while an individ­ual vacancy diffuses into the simulation volume every 1 x 104s (~3h). AKMC simulations have also been performed to study the effect of varying the cascade introduction rate from 1.13 x 10~3 to 1.13 x 10~7 cascades per second (dpa rates from 1 x 10~9 to 1 x 10_ 3dpas~ ). The simulated conditions should be compared to those experienced by RPVs in light water reactors, namely from 8 x 10~12 to 8 x 10~n dpas-1, and to model alloys irradiated in test reac­tors, which are in the range of 10~9-10~1°dpas~1.

Figures 2 and 3 show representative snapshots of the vacancy and Cu solute atom distributions as a function of time and dose at 290 °C. Note, only the Cu atoms that are part of vacancy or Cu atom clusters are presented in the figure. The main simulation volume consists of 2 x 106 atoms (100a0 x 100a0 x 100a0) of which 6000 atoms are Cu (0.3 at.%). Figure 2 demonstrates the aging evolution of a single cascade (increasing time at fixed dose prior to introducing additional diffusing vacancies or new cascade), while Figure 3 demonstrates the overall evolution with increasing time and dose. The aging of the single cascade is representative of the average behavior observed, although the number and size distribution of vacancy-Cu clusters do vary consider­ably from cascade to cascade. Further, Figure 2 is representative of the results obtained with the previ­ous AKMC models of Odette and Wirth21 and Becquart and coworkers,24,26 which demonstrated the formation and subsequent dissolution of vacancy-Cu clusters. Figure 3 represents a significant extension of that previous work21,24,26 and demonstrates the for­mation of much larger Cu atom precipitate clusters that results from the longer term evolution due to multiple cascade damage in addition to radiation enhanced diffusion.

Figure 2(a) shows the initial vacancy configura­tion from an aged 20-keV cascade. Within 200 ps at 290 °C, the vacancies begin to diffuse and cluster, although no vacancies have yet reached the cell boundary to annihilate. Eleven of the initial vacancies remain isolated, while thirteen small vacancy clusters rapidly form within the initial cascade volume. The vacancy clusters range in size from two to six vacancies. At this stage, only two of the vacancy clus­ters are associated with copper atoms, a divacancy

Figure 2 Representative vacancy (red circles) and clustered Cu atom (blue circles) evolution in an Fe-0.3% Cu alloy during the aging of a single 20keV displacement cascade, at (a) initial (200 ns), (b) 2 ms, (c) 48 ms, (d) 55 ms, (e) 83 ms, and (f) 24.5 h.

cluster with one Cu atom and a tetravacancy cluster with two Cu atoms. From 200 ps to 2 ms, the vacancy cluster population evolves by the diffusion of isolated vacancies through and away from the cascade region, and the emission and absorption of isolated vacancies in vacancy clusters, in addition to the diffusion of the small di-, tri-, and tetravacancy clusters. Figure 2(b) shows the configuration about 2 ms after the cascade. By this time, 14 of the original vacancies have dif­fused to the cell boundary and annihilated, while 38 vacancies remain. The vacancy distribution includes six isolated vacancies and seven vacancy clusters, ranging in size from two divacancy clusters to a ten vacancy cluster. The number of nonisolated copper atoms has increased from 223 in the initial random distribution to 286 following the initial 2 ms of cas­cade aging.

The evolution from 2 to 48.8 ms involves the dif­fusion of isolated vacancies and di — and trivacancy clusters, along with the thermal emission of vacancies from the di — and trivacancy clusters. Over this time, 7 additional vacancies have diffused to the cell boundary and annihilated, and 20 additional Cu atoms have been incorporated into Cu or vacancy clusters. Figure 2(c) shows the vacancy and Cu clus­ter population at 48.8 ms, which now consists of three isolated vacancies and four vacancy clusters, includ­ing a 4V-1Cu cluster, a 6V-4Cu cluster, a 7V cluster, and an 11V—1Cu cluster. Figure 2(d) and 2(e) shows the vacancy and Cu cluster population at 54.8 and

82.8 ms, respectively. During this time, the total num­ber of vacancies has been further reduced from 31 to 21 of the original 52 vacancies, the vacancy cluster
population has been reduced to three vacancy clus­ters (a 4V—1Cu, 7V, and 9V—1Cu), and 30 additional Cu atoms have incorporated into clusters because of vacancy exchanges.

Over times longer than 100 ms, the 4V-1Cu atom cluster migrates a short distance on the order of 1 nm before shrinking by emitting vacancies, while the 7 V and 9V-1Cu cluster slowly evolve by local shape rearrangements which produces only limited local diffusion. Both the 7V and 9V-1Cu cluster are ther­modynamically unstable in dilute Fe alloys at 290 °C and ultimately will shrink over longer times. The vacancy and Cu atom evolution in the AKMC model is now governed by the relative rate ofvacancy cluster dissolution, as determined from the ‘pulsing’ algo­rithm, and the rate of new displacement damage and the diffusing supersaturated vacancy flux under irra­diation. Figure 2(f) shows the configuration about

8.8 x 104s (^24 h) after the initial 20keV cascade. Only 17 vacancies now exist in the cell, an isolated vacancy which entered the cell following escape from a 500 eV recoil introduced into a neighboring cell plus two vacancy clusters, consisting of 7V-1Cu and 9V—1Cu. Three hundred and forty-five Cu atoms (of the initial 6000) have been removed from the super­saturated solution following the initial 24 h of evolu­tion, mostly in the form ofdi — and tri-Cu atom clusters.

Figure 3(a) shows the configuration at about 0.1mdpa (0.097mdpa) and a time of 7.1 x 106s (^82 days). Ten vacancies exist in the simulation cell, consisting of eight isolated vacancies and one 2V
cluster, while 807 Cu atoms have been removed from solution in clusters, although the Cu cluster size distri­bution is clearly very fine. The majority of Cu clusters contain only two Cu atoms, while the largest cluster consists of only five Cu atoms. Figure 3(b) shows the configuration at a dose of 0.33 mdpa and time of

2.1 x 107s (245 days). Only one vacancy exists in the simulation volume, while 1210 Cu atoms are now part of clusters, including 12 clusters containing 5 or more Cu atoms. Figure 3(c) shows the evolution at 1 mdpa and 7.2 x 107s (^2.3 years). Again, only one vacancy exists in the simulation cell, while 1767 Cu atoms have been removed from supersaturated solution. A handful of well-formed spherical Cu clusters are visible, with the largest containing 13 Cu atoms. With increasing dose, the free Cu concentration in solution continues to decrease as Cu atoms join clusters and the average Cu cluster size grows. Figure 3(d) and 3(e) shows the clustered Cu atom population at about 2 and

4.4 mdpa, respectively. The growth of the Cu clusters is clearly evident when Figure 3(d) and 3(e) is com­pared. At a dose of 4.4 mdpa, 48 clusters contain more than 10 Cu atoms, and the largest cluster has 28 Cu atoms. The accumulated dose of 5.34 mdpa is shown in Figure 3(f). At this dose, more than one-third of the available Cu atoms have precipitated into clusters, the largest of which contains 42 Cu atoms, corresponding to a precipitate radius of ^0.5 nm.

Figure 4 shows the size distribution of Cu atom clusters at 5.34 mdpa, corresponding to the configu­ration shown in Figure 3(f). The vast majority of the

Cu clusters consist of di-, tri-, tetra-, and penta-Cu atom clusters. However, as shown in the inset of Figure 4 and as visible in Figure 3(f), a significant number of the Cu atom clusters contain more than five Cu atoms. Indeed, 29 clusters contain 15 or more Cu atoms (a number density of

1.2 x 1024m~3), which corresponds to a cluster con­taining a single atom with all first and second near­est neighbor Cu atoms and a radius of 0.29 nm. An additional 45 clusters contain at least nine Cu atoms (atom + all first nearest neighbors), while 9 clusters contain 23 or more atoms (number density of

3.8 x 1023 m~3). This AKMC simulation is currently continuing to reach higher doses. However, the ini­tial results are consistent with experimental obser­vations and show the formation of a high number density of Cu atom clusters, along with the contin­ual formation and dissolution of 3D vacancy-Cu clusters.

Figure 5 shows a comparison of varying the dose rate from 10~9 to 10~13 dpas-1. Each simulation was performed at a temperature of 290 °C and introduced additional vacancies into the simulation volume at the rate of 10~4s-1. The effect of increasing dose

. ,v4i. ‘• ".-і-. ‘:,

v> Л %

4′- *■

, ЛЛГ4*F., * . 4

‘І/ч

(b)

(c)

Figure 5 Comparison of the representative vacancy (red) and clustered Cu atom (blue) population at a dose of ~1.9 mdpa and 290 "Casa function of dose rate, at (a) 10~11 dpas~1, (b) 10~9dpas~ (c) 10^13dpas^1.

rate at an accumulated dose of 1.9 mdpa is especially pronounced when comparing Figure 5(c) (10~ 3dpa s-1) with Figure 5(a) (10~n dpas-1) and Figure 5(b) (10~9dpas-1). At the highest dose rate, a substan­tially higher number density of small 3D vacancy clusters is observed, which are often complexed with one or more Cu atoms. Vacancy cluster nucle — ation occurs during cascade aging (as described in Figure 2) and is largely independent of dose rate, but cluster growth is dictated by the cluster(s) thermal lifetime at 290 °C versus the arrival rate of additional vacancies, which is a strong function of the damage rate and vacancy supersaturation under irradiation. Thus, the higher dose rates produce a larger number of vacancies arriving at the vacancy cluster sinks, resulting in the noticeably larger number of growing vacancy clusters. Also, there is a corresponding decrease in the amount of Cu removed from the solution by vacancy diffusion. In contrast, the effect of decreasing dose rate is greatly accelerated Cu precipitation. Already at 1.9 mdpa, a number of large Cu atom clusters exist at a dose rate of 10- dpas~ , with the largest containing 35 Cu atoms, as shown in Figure 5(c). The increased Cu clustering caused by a decrease in dose rate results from a reduction in the number ofcascade vacancy clusters, which serve as vacancy sinks. Thus, a higher number of free or isolated vacancies are available to enhance Cu diffusion required for the clustering and precipi­tation of copper. While these flux effects are antici­pated and have been predicted in rate theory calculations performed by Odette and coworkers,78,79 the spatial dependences of cascade production and microstructural evolution, in addition to correlated diffusion and clustering processes involving multiple vacancies and atoms are more naturally modeled and visualized using the AKMC approach.

While the results just presented in Figures 2-5 have shown the formation of subnanometer Cu-vacancy clusters and larger growing Cu precipitate clusters that result from AKMC simulations, which only consider vacancy-mediated diffusion, Becquart and coworkers have shown that Cu atoms in tensile posi­tions can trap SIAs and therefore the Cu clustering behavior may also be influenced by interstitial — mediated transport. Ngayam Happy and coworkers63,86 have developed another AKMC model to model the behavior of FeCu under irradiation. In this model, diffusion takes place via both vacancy and self­interstitial atoms jumps on nearest neighbor sites. The migration energy of the moving species is also deter­mined using eqn [4], where the reference activation
energy Ea0 depends only on the type of the migrating species. Ea0 has been set equal to:

• the ab initio vacancy migration energy in pure Fe when a vacancy jumps towards an Fe atom (0.62 eV);

• the ab initio solute migration energy in pure Fe when a vacancy jumps towards a solute atom (0.54 eV for Cu); and

• the ab initio migration-600 rotation energy of the migrating atom in pure Fe when a dumbbell migrates (0.31 eV).

Ei and Ef are determined using pair interactions, according to the following equation:

E = e(<)(Sj — — Sk) + Edumb [5]

i=1,2 j<k

where i equals 1 or 2 and corresponds to first or second nearest neighbor interactions, respectively, and where j and k refer to the lattice sites and Sj (respectively Sk) is the species occupying site j (respectively k): Sj in {X, V} where X = Fe or Cu. A more detailed description of the model can be found in Ngayam Happy et al63

In this study, various Cu contents were simulated (0, 0.18, 0.8, and 1.34 at.%) at three different tem­peratures (300, 400, and 500 0C). Without going into too much detail here, one can state that these AKMC simulation results are qualitatively similar to those presented in Figures 2-5, which showed the formation of small, vacancy-solute clusters and copper enriched cluster/precipitate formation at 300 0C. Similarly, the effect of decreasing dose rate in high Cu content alloys was also found to accelerate Cu precipitation.

This model does show that the formation of the Cu clusters/precipitates during neutron irradiation takes place via two different mechanisms depending upon the Cu concentration. In a highly Cu supersaturated matrix, precipitation is accelerated by irradiation, whereas in the case of low Cu contents, Cu precipi­tates form by induced segregation on vacancy clusters.

The influence of temperature was investigated for an Fe-0.18wt% Cu alloy irradiated at a flux of

2.3 x 10-5dpas-1. At 400 and 500 0C, neither Cu precipitates nor Cu-vacancy clusters were formed, in agreement with the results of Xu et a/.87 At these temperatures, the model indicates that the vacancy clusters are not stable and induced segregation is thus hindered. Another interesting result obtained with this model is that the presence of Cu atoms in the matrix was found to decrease the point defect cluster sizes because of the strong interactions of Cu with both vacancies and SIAs.

The liquid phases in nuclear fuels are important to model so that the phase equilibria can be completely assessed through comparison of experimental and computed phase diagrams. The availability of solidus and liquidus information also provides necessary boundaries for modeling the solid-state behavior. Finally, safety analysis requirements with regard to the potential onset of melting will benefit from accu­rate representations of the complex liquids.

Ideal, regular/subregular, or Bragg-Williams formulations are not very successful in representing metal and especially oxide liquids where there are strong interactions between constituents. The CEF model is designed for fixed lattice sites, and thus it too will not handle liquids. The issues for these complex liquids involve the short-range ordering that generally occurs and its effect on the form of the Gibbs free energy expressions. One approach to dealing with the issue of these strong interactions is the modified associate species method.

The modified associate species technique for crystalline materials was discussed to an extent in Section 1.17.4.2. Its application to, for example, oxide melts has been more broadly covered recently by Besmann and Spear39 with much of the original development by Hastie and coworkers.40-43 The approach assumes that the liquid can be modeled by an ideal solution of end-member species together with intermediate species. The modified term refers to the fact that an ideal solution cannot represent a miscibility gap in the liquid as that requires repulsive (positive) interaction energy terms. Thus, when a mis­cibility gap needs to be included, interaction energies between appropriate associate species are added to the formulation.

In the associate species approach, the system standard Gibbs free energy is simply the sum of the constituent end-member and associate free energies, for example, A, B, and A2B, where inclusion of the A2B associate is found to reproduce the behavior well,

G° = XaGa + XbGb + Xa2bGA2b [16]

Consequently, ideal mixing among end members and associates generates the entropy contribution

Gid = RT (Xa ln Xa + Xb ln Xb + Xa2b ln Xa2b) [17]

Should a nonideal term providing positive interac­tion energies be needed to properly address a misci­bility gap, it would be added into the total Gibbs free energy as in eqn [10]. For example, for an interaction between A and A2B in the Redlich-Kister-Muggianu formulation the excess term is expressed as

Gex = XaXa2bELk, A:A2B(XA — Xa2b)* [18]

k

The associate species are typically selected from the stoichiometry of intermediate crystalline phases, but others as needed can be added to accurately reflect the phase equilibria even when no stable crys­talline phases of that stoichiometry exist. Gibbs free energies for these species can be derived from fits to the phase equilibria and other data following the CALPHAD method with first estimates gener­ated from crystalline phases of the same stoichiome­try or weighted sums of existing phases when no stoichiometric phase exists. The application of the method for the liquid phase in the Na2O-Al2O3 is seen in the computed phase diagram in Figure 8. For this system, the associate species required to represent the liquid were only Na2O, NaAlO2, (1/3) Na2Al4O7, and Al2O3. In nuclear fuel systems, Chevalier et al44 applied an associate species approach using the components O, U, and O2U, although it deviated from the associate species approach in using binary interaction parameters in a Redlich-Kister-Muggianu form. The computed

phase diagram showing agreement with the liquidus/ solidus data is seen in Figure 9.

The use of the modified associate species model with ternary and higher order systems can require the use of ternary or possibly quaternary associates. Another issue with the modified associate species approach is that in the case of a highly ordered solution which requires an overwhelming content of an associate compared to an end-member, the rela­tions do not follow what should be Raoult’s law for dilute solutions. At the other extreme, it is also appar­ent that in the case of essentially zero concentration of associates, the relationships do not default to an ideal solution as one would expect.

Ion irradiations have been critical to the development of both our fundamental and applied understanding of radiation effects. As discussed in Sections 1.07.2 and 1.07.3, it is the flexibility of such irradiations and our firm understanding of atomic collisions in solids that afford them their utility. Principally, ion irradia­tions have enabled focused studies on the isolated effects of primary recoil spectrum, defect displace­ment rate, and temperature. In addition, they have provided access to the fundamental properties of point defects, defect creation, and defect reactions. In this section, we highlight a few key experiments that illustrate the broad range of problems that can be addressed using ion irradiations. We concentrate our discussion on past ion irradiations studies that have provided key information required by modelers in their attempts to predict materials behavior in existing and future nuclear reactor environments, and particularly information that is not readily available from neutron irradiations. In addition, we include a few comparative studies between ion and neutron irradiations to illustrate, on one hand,
the good agreement that is possible, while on the other, the extreme caution that is necessary in extrapolating results of ion irradiations to long-term predictions of materials evolution in a nuclear environment.