The embedded atom and FS potentials fall into a general class of potentials of the form:

U, =£ F

i

with a many-body cohesive part and a two-body repulsion. Both f(rj) and V(r, y) are short ranged, so MD with these potentials is at worst only as costly as a simple pair potential (computer time is propor­tional to number of particles).

These models are sometimes referred to as glue potentials,16 the many-body F term being thought of as describing how strongly an atom is held by the electron ‘glue’ provided by its environment.

The pragmatic approach to fitting in all glue schemes is to regard the pair potential as repulsive at short-range with long-range Friedel oscillations. Compared with most pair potential approaches, this is unusual in that the repulsive term is longer ranged than the cohesive one.

1.10.7.1 Crystal Structure

According to the tight-binding theory on which the FS potentials are based, the relative stability of bcc and fcc is determined by moments above the second, which in turn relate to three center and higher hops. These third and higher moments effects are explicitly absent in second-moment models, and so by implica­tion, the correct physics of phase stability is not contained in them. There is no such clear result in the derivation of the EAM; however, since the forms are so similar, the same problem is implicit.

In glue models, energy is lowered by atoms having as many neighbors as possible; thus, fcc, hcp, and bcc crystal structures (and their alloy analogs) are nor­mally stable (see Table 1); bcc is normally stable in potentials when the attractive region is broad enough
to include 14 neighbors, fcc/hcp are stable for nar­rower attractive regions in which only the eight near­est bcc-neighbors contribute significantly to the bonding. Indeed, without second neighbor interac­tions, bcc is mechanically unstable to Bain-type shear distortion. The fcc-hcp energy difference is related to the stacking fault energy: it is common to see MD simulations with too small an hcp-fcc energy differ­ence producing unphysically many stacking faults and over widely separated partial dislocations.

Phase transitions are observed in some potentials. As free energy calculations are complicated and time consuming,17 it is impractical to use them directly in fitting — one would require the differen­tial of the free energy with respect to the potential parameters, and this could only be obtained numer­ically. Consequently, most potentials are only fitted to reproduce the zero temperature crystal structure, and high-temperature phase stability is unknown for the majority of published potentials. One counter­example is in metals such as Ti and Zr, where the bcc structure is mechanically unstable with respect to hcp, but becomes dynamically stabilized at high temperatures. Here, the transition temperature is directly related to a single analytic quantity: the energy difference between the phases. Although about half of this difference comes from electronic entropy,18 which suggests a temperature-dependent potential, phase transition calculations have been explicitly included in some recent fits.19 The case of iron is also anomalous, as the phase transition is related to changes in the magnetic structure.

Y. N. Osetsky

Oak Ridge National Laboratory, Oak Ridge, TN, USA D. J. Bacon

The University of Liverpool, Liverpool, UK Published by Elsevier Ltd.

1.12.1 Introduction 334

1.12.2 Radiation Effects on Mechanical Properties 334

1.12.2.1 Radiation-Induced Obstacles to Dislocation Glide 334

1.12.2.2 Effects on Mechanical Properties 335

1.12.3 Method 336

1.12.3.1 Why Atomic-Scale Modeling? 336

1.12.3.2 Atomic-Level Models for Dislocations 336

1.12.3.3 Input Parameters 338

1.12.3.4 Output Information 338

1.12.4 Results on Dislocation-Obstacles Interaction 339

1.12.4.1 Inclusion-Like Obstacles 339

1.12.4.1.1 Temperature T — 0K 339

1.12.4.1.2 Temperature T> 0K 342

1.12.4.2 Dislocation-Type Obstacles 345

1.12.4.2.1 Stacking fault tetrahedra 345

1.12.4.2.2 Dislocation loops 348

1.12.4.3 Microstructure Modifications due to Plastic Deformation 352

1.12.5 Concluding Remarks 353

References 355

1.12.1 Introduction

Structural materials in nuclear power plants suffer a significant degradation of their properties under the intensive flux of energetic atomic particles (see Chapter 1.03, Radiation-Induced Effects on Microstructure). This is due to the evolution of microstructures associated with the extremely high concentration of radiation-induced defects. The high supersaturation of lattice defects leads to microstruc­tures that are unique to irradiation conditions. Irra­diation with high energy neutrons or ions creates initial damage in the form of displacement cascades that produce high local supersaturations of point defects and their clusters (see Chapter 1.11, Pri­mary Radiation Damage Formation). Evolution of the primary damage under the high operating tem­perature (^600 K to >1000 K) leads to a micro­structure containing a high concentration of defect clusters, such as voids, dislocation loops (DLs), stack­ing fault tetrahedra (SFTs), gas-filled bubbles, and precipitates, and an increase in the total dislocation network density (see Chapter 1.13, Radiation Dam­age Theory; Chapter 1.14, Kinetic Monte Carlo Simulations of Irradiation Effects and Chapter 1.15, Phase Field Methods). These changes affect material properties, including mechanical ones, which are the subject of this chapter.

A general theory of radiation effects has not yet been developed, and currently the most promising way to predict materials behavior is based on multi­scale materials modeling (MMM). In this framework, phenomena are considered at the appropriate length and times scales using specific theoretical and/or modeling approaches, and the different scales are linked by parameters/mechanisms/rules to provide integrated information from a lower to a higher level.

Research on the mechanical properties of irradiated materials, a topic of crucial importance for engineering solutions, provides a good example of this. The lowest level treats individual atoms by first principles, ab initio methods, by solving Schrodinger’s equation for moving electrons and ions. Calculations based on electron den­sity functional theory (DFT) (see Chapter 1.08, Ab Initio Electronic Structure Calculations for Nuclear Materials) and its approximations, such as bond order potentials (BOPs), can consider a few hundred atoms over a very short time of femtoseconds to picoseconds. Delivery of the resulting information to higher level models can be achieved through effective interatomic potentials (IAPs) (see Chapter 1.10, Interatomic Potential Development), in which the adjustable parameters are fitted to the basic chemical and struc­tural properties obtained ab initio. IAPs are required for atomic-scale modeling methods such as molecular stat­ics (MS) and molecular dynamics (MD), which are used to simulate millions of atoms. Time spanning nanose­conds to microseconds can be simulated by MD if the number of atoms is not large (see Section 1.12.3.3). This level can provide properties of point and extended defects and interactions between them (see Chapter 1.09, Molecular Dynamics). For mechanical properties, important interactions are between moving dislocations, which are responsible for plasticity, and defects created by irradiation. Mechanisms and para­meters determined at this level can then inform dislo­cation dynamics (DD) models based on elasticity theory of the continuum (see Chapter 1.16, Disloca­tion Dynamics). DD models can simulate processes at the micrometer scale and mesh with the mechanical properties of larger volumes of material used in finite elements (FEs) methods, that is, realistic models for the design of core components. In this chapter, we consider direct interactions at the atomic scale between moving dislocations and obstacles to their motion.

The structure of the chapter is as follows. First, we summarize the main features of the irradiation microstructure of concern. Then we provide a short description of atomic-scale methods applied to dislo­cation modeling, bearing in mind the details pre­sented in Chapter 1.09, Molecular Dynamics. This is followed by a review of important results from simulations of the interaction between disloca­tions and obstacles. We then describe how dislocations modify microstructure in irradiated metals. Finally, we indicate some issues that will hopefully be resolved by atomic-scale modeling in the near future. Our main aim is to give the reader a general picture of the phenomena involved and encourage further research in this area. The following sources1-4 provide a more general and deeper understanding of disloca­tions and modeling of plasticity issues.

Mutual annihilation of PDs happens either by direct interaction between single vacancies and SIAs in the matrix or within a certain type of neutral sink which we call ‘saturable.’ The fluxes of vacancies and SIAs to them are equal. An example of such sinks is vacancy loops, which were considered in the frame­work of the BEK model29 and PBM,22 that is, in the case where the vacancy clusters are generated in cascades. The BEK model is not discussed further in the present work because it does not correspond to any realistic situation in solids under irradiation; vacancy clustering in cascades is always accompanied with the SIA clustering, which is accounted for in the framework of the PBM but not in the BEK model.

The balance equations in the case considered are as follows

G — mRDi C Cv — lN A C — Cv — zdpdDvCv = 0

g — mRDi c Cv — iNa ci — A Ci — zd rdDi ci = 0 [99]
where kN2 is the strength of neutral sinks. Note that absorption rate of both vacancies and SIAs in eqn [99] is described by the same quantity, kNDiQ, which reflects neutrality of this sink with respect to vacancies and SIAs.115,116

The defect concentrations and swelling rate are

G 1 1

+ ZvPd 1 + JR 1 + JN

d5 = B k2zvdrd _J_________________ L_

df d (12 + Zdrd)2 1 + JR 1 + JN

where

4PrG

(k2 + zd Pd + kN )2dv

|2

kN 12 + zd Pd

In the absence of an effect on the sink structure, mutual recombination reactions are important at low temperature, when the vacancy diffusion is slow, and for high defect production rates, when the vacancy concentration is sufficiently high to provide higher sink strength for SIAs than that of existing extended defects. This can be expressed mathematically by an inequality JR > 1 or more explicitly as a temperature boundary IbT < £m/ln[2Dv0(k^ + Z^pd + )2/^rG] where Dv0 is the preexponential factor in the vacancy diffusion coefficient and Evm is the effective activation energy for the vacancy migration. In practice, this situation is unlikely to occur because the radiation — induced sink strength rapidly increases at low tem­peratures. In this case recombination at sinks is of greater importance.

One of the important aspects that recombination reactions introduce to microstructural evolution is the appearance of a temperature dependence; at low temperatures, an increase of the swelling rate with increasing temperature is predicted, which is also observed experimentally in the fcc-type materials. The question of whether it was possible to explain the experimental reduction of swelling rate with decreasing temperature by recombination was addressed.29 It was found that the observed tempera­ture effect on swelling rate was much stronger than predicted by recombination alone.

The impact of neutral sinks on swelling rate is significant when they represent the dominant sink

in the system: kN ^ kc + Z^pA. The swelling rate in the case is given by

dS= B k2ZydPd

df d (k;: + ZdPd)(k2 + ZdPd + kN)

k;zdPd

~ d (k; + zdpd)kN [102]

Such a situation may occur, for example, at low enough temperature, when the thermal stability of vacancy loops and SFTs becomes high enough, lead­ing to their accumulation up to extremely high con­centrations. Another possibility is when a high density (about 1024m-3) of second phase particles exists, as in the oxide dispersion strengthened (ODS) steels.

models of plastic deformation are extensively used in engineering practice, their range of application is limited by the underlying database. The reliability of continuum plasticity descriptions is dependent on the accuracy and range of available experimental data. Under complex loading situations, however, the database is often hard to establish. Moreover, the lack of a characteristic length scale in continuum plasticity makes it difficult to predict the occurrence of critical localized deformation zones. In small volumes, or in situations where submicrometer reso­lution of material deformation is required (e. g., in thin films, micropillars, near dislocations, grain boundaries, etc.), the use of continuum plasticity models would be questionable.

Although homogenization methods have played a significant role in determining the elastic properties of new materials from their constituents (e. g., composite materials), the same methods have failed to describe plasticity. It is widely appreciated that plastic strain is fundamentally heterogenous, displaying high strains concentrated in small material volumes, with virtually undeformed regions in between. Experimental obser­vations consistently show that plastic deformation is heterogeneous at all length scales. Depending on the deformation mode, heterogeneous dislocation struc­tures appear with definitive wavelengths. A satisfactory

description of realistic dislocation patterning and strain localization has been rather elusive. Attempts aimed at this question have been based on statistical mechanics, reaction-diffusion dynamics, or the theory of phase transitions. Much of the effort has aimed at clarifying the fundamental origins of inhomogeneous plastic deformation. On the other hand, engineering descrip­tions of plasticity have relied on experimentally ver­ified constitutive equations.

Planar dislocation arrays are formed under mono­tonic stress deformation conditions, and are com­posed of parallel sets of dislocation dipoles. While persistent slip bands (PSBs) are found to be aligned in planes with their normal parallel to the direction of the critical resolved shear stress (CRSS), planar arrays are aligned in the perpendicular direction. Dislocation cell structures, on the other hand, are honeycomb configurations in which the walls have high dislocation density, while the cell interiors have low dislocation density. Cells can be formed under both monotonic and cyclic deformation conditions. However, dislocation cells under cyclic deformation tend to appear after many cycles. Direct experimental observations of these structures have been reported for many materials.

Two main approaches have been advanced to model the mechanical behavior in this meso length scale. The first is based on statistical mechanics meth­ods. In this approach, evolution equations for statisti­cal averages (and possibly for higher moments) are to be solved for a complete description of the deforma­tion problem. The main challenge in this regard is that, unlike the situation encountered in the develop­ment of the kinetic theory of gases, the topology of interacting dislocations within the system must be included. The second approach, commonly known as dislocation dynamics (DD), was initially motivated by the need to understand the origins of hetero­geneous plasticity and pattern formation. An early variant of this approach (the cellular automata) was first developed by Lepinoux and Kubin,1 and this was followed by the proposal of DD by Ghoniem and Amodeo.2-4

In these early efforts, dislocation ensembles were modeled as infinitely long and straight in an isotro­pic, infinite elastic medium. The method was further expanded by a number of researchers,5 with appli­cations demonstrating simplified features of defor­mation microstructure. DD has now become an important computer simulation tool for the descrip­tion of plastic deformation at the micro — and meso — scales (i. e., the size range of a fraction of a micrometer to tens of micrometers). The method is based on a hierarchy of approximations that enable the solution of relevant problems with today’s computational resources. In its early versions, the collective behav­ior of dislocation ensembles was determined by direct numerical simulations of the interactions between infinitely long, straight dislocations.5 Recently, sev­eral research groups have extended the DD method­ology to the more physical, yet considerably more complex case of three-dimensional (3D) simulation. The method can be traced back to the concepts of internal stress fields and configurational forces. The more recent development of 3D lattice DD by Kubin and coworkers6-8 has resulted in greater confidence in the ability of DD to simulate more complex defor­mation microstructures. More rigorous formulations of 3D DD have contributed to its rapid development and applications in many systems.9-15

Many experimental observations have shown that neutron irradiation of metals and alloys at temper­atures below recovery stage V causes a substantial increase in the upper yield stress (radiation hardening) and, beyond a certain dose level, induces a yield drop and plastic instability (see Chapter 1.03, Radiation — Induced Effects on Microstructure and Chapter 1.04, Effect of Radiation on Strength and Ductility of Metals and Alloys). Furthermore, the postdeformation microstructure of a specimen showing an upper yield point has demonstrated two significant features. First, the onset of plastic deformation is generally found to coincide with the formation of‘cleared’ channels, where practically all plastic deformation takes place. The sec­ond feature refers to the fact that the material volume in between cleared channels remains almost undeformed (i. e., no new dislocations are generated during deforma­tion). In other words, the initiation of plastic defor­mation in these irradiated materials occurs in a very localized fashion. This specific type ofplastic flow local­ization is considered to be one of many possible plastic instabilities in both irradiated and unirradiated materials.16,17

A theory of radiation hardening was proposed by Seeger18 in terms of dislocation interaction with radiation-induced obstacles (referred to as depleted zones). Subsequently, Foreman19 performed computer simulations of loop hardening, in which the elastic interaction between the dislocation and loops was neglected. The model is based on Orowan’s mecha­nism, which assumes that the obstacles are indestructi­ble. In this view of matrix hardening, the stress necessary to overcome localized interaction barriers leads to the increase in the yield strength, while long-range elastic interactions are completely ignored. These short-range dislocation-barrier interaction models lead to an estimate of the increase of the CRSS of the form: At = a mb/l, where a is a numerical constant represent­ing obstacle strength, m the shear modulus, b the mag­nitude of the Burgers vector, and l is the average interobstacle distance. Kroupa,20 on the other hand, viewed hardening as a result of the long-range elastic interaction between slip dislocations and prismatic loops. In their model of friction hardening, the force necessary to move a rigid, straight dislocation on its glide plane past a prismatic loop was estimated. In these two classes of hardening models, all dislocation sources are assumed to be simultaneously activated at the yield point, and plastic deformation is assumed to be homogeneous throughout the material volume. Thus, they do not address the physics of plastic flow localization. Singh et al. introduced the concept of ‘stand-off distance’ for the decoration of dislocations with small self-interstitial atom (SIA) loops, and pro­posed the ‘cascade-induced source hardening’ (CISH) model, in analogy with Cottrell atmospheres.

The central question of the formation of dislocation decoration was treated analytically by Trinkaus etal23 Subsequent detailed elasticity calculations showed that glissile defect clusters that approach dislocation cores within a stand-off distance are absorbed, while clusters can accumulate just outside this distance. Trinkaus etal. concluded that dislocation decoration is a consequence of defect cluster mobility and trapping in the stress field of grown-in dislocations. The CISH model was used to calculate the stress necessary to free decorated dislocations from the atmosphere ofloops around them (upper yield, followed by yield drop), so that these freed dislocations can act as dislocation sources. The CRSS increase in irradiated Cu was shown by Singh et al21 to be given by: At ~ 0.1m(b/l)(d/y)2, where m is the shear modulus, and b, l, d, andy are the Burgers vector, interdefect distance, defect diameter, and stand­off distance, respectively. Assuming that y ~ d, and l~ 35b, they estimated At/m ~ 5 x 10~3. The phe­nomenon of yield drop was proposed to result from the unpinning of grown-in dislocations, decorated with small clusters or loops of SIAs. Neutron irradiation of pure copper at 320 K leads to an increase in the upper yield stress and causes a prominent yield drop and the initiation of plastic flow localization.

In this chapter, we first present a description of the mathematical and computational foundations of DD. We then assess the physical mechanisms, that are responsible for the initiation of plastic insta­bility in irradiated face-centered cubic (fcc) metals by detailed numerical simulations of the interaction between dislocations and radiation-induced defect clusters. Two distinct problems that are believed to cause the onset of plastic instability are addressed in the present study. First, we aim to determine the mechanisms of dislocation unlocking from defect cluster atmospheres as a result of the long-range elastic interaction between dislocations and sessile prismatic interstitial clusters situated just outside the stand-off distance. The main new feature of this analysis is that dislocation deformation is explicitly considered during its interaction with SIA clusters. Second, we investigate the mechanisms of structural softening in flow channels as a consequence of dis­location interaction with stacking fault tetrahedra (SFTs). Based on these numerical simulations, a new mechanism of channel formation is proposed, and the magnitude of radiation hardening is also computed. As an important application of DD simu­lations, we then present a model that describes the shift in the ductile-to-brittle transition temperature (DBTT) as a result of neutron irradiation, and com­pare this physically based model to experimental data. Finally, we discuss the future outlook for DD simulations and the role they may play in under­standing radiation effects on the mechanical proper­ties of structural materials.

The thermodynamic factor in eqn [12] is propor­tional to the second derivative of the Gibbs free energy G of the alloy, with respect to the molar fraction of one of the components. It can be calcu­lated on the basis of thermodynamic data. A database such as CALPHAD61 builds free-energy composition functions of the alloy phases from thermodynamic measurements (specific heats, activities, etc.). When available, the phase diagrams are used to refine and/or to assess the thermodynamic model. Although the CALPHAD free-energy functions are sophisti­cated functions of temperature and composition, it is interesting to study the simple case of a regular solu­tion model. In the case of a binary alloy A1—CBC with a clustering tendency, the Gibbs free energy is equal to

G = 2kBTcC(1 — C) + kBTC ln(C)

+ kBT(1 — C) ln(1 — C) [15]

where Tc is the critical temperature and C is the alloy composition. The regular solution approximation leads to a concentration-dependent thermodynamic factor equal to

F = 1 — 4C(1 — C)T [16]

where concentration C now corresponds to a local concentration of B atoms, which varies in space and time.

As FPs are produced in isolation during electron irradiation, the properties of single point defects and their interactions with impurities and sinks can be systematically investigated. An example is shown in Figure 11(a), where the results of low — temperature isochronal annealing of Cu are shown following 1.4 MeV electron irradiation at 6K.2 Recovery is observed to occur in ‘stages.’ These stud­ies have revealed that interstitial atoms become mobile at very low temperatures, always below 100 K, in so-called Stage I, while vacancies become mobile at higher temperatures, Stage III. The various substages IA-IE seen in Figure 11(a) arise from the interaction between interstitial-vacancy pairs, which are produced in close proximity. Stage IE refers to the free migration of interstitials in the lattice, away from its own vacancy, and annihilation at distant vacancies; these interstitials are freely migrating as discussed earlier. For comparison, Stage I annealing of Cu following neutron irradiation is shown in Figure 11(b). Notice that the close pair substages are suppressed during neutron irradiation, illustrat­ing the dramatic difference in the defect production process for these types of irradiation. Similarly, annealing studies on electron-irradiated Al doped with Mg or Ga impurities are shown in Figure 12.22 For these, it is observed that Stage I recovery is suppressed as interstitials trap at impurities and do not recombine. The recovery at higher temperature, in Stage II, reveals distinct subannealing stages.

These annealing stages are generally attributed to either the interstitial dissociating from the impurity, or the interstitial-impurity complex migrating to a vacancy or a defect sink. Migrating interstitial-solute complexes lead to segregation. A compilation of the properties of point defects for many metals, and their interactions with impurities can be found in Ehrhart.2 This information has played a crucial role in develop­ing an understanding of radiation damage in more complex engineering alloys and under more complex irradiation conditions.

In the simplest physical terms, MD may be charac­terized as a method of ‘particle tracking.’ Operation­ally, it is a method for generating the trajectories of a system of N particles by direct numerical integration

of Newton’s equations of motion, with appropriate specification of an interatomic potential and suitable initial and boundary conditions. MD is an atomistic modeling and simulation method when the particles in question are the atoms that constitute the material of interest. The underlying assumption is that one can treat the ions and electrons as a single, classical entity. When this is no longer a reasonable approxi­mation, one needs to consider both ion and electron motions. One can then distinguish two versions of MD, classical and ab initio, the former for treating atoms as classical entities (position and momentum) and the latter for treating separately the electronic and ionic degrees of freedom, where a wave func­tion description is used for the electrons. In this chapter, we are concerned only with classical MD. The use of ab initio methods in nuclear materials research is addressed elsewhere (Chapter 1.08, Ab Initio Electronic Structure Calculations for Nuclear Materials). Figure 2 illustrates the MD sim­ulation system as a collection of N particles contained in a volume O. At any instant of time t, the particle coordinates are labeled as a 3N-dimensional vector, r3N(t) = {ri(t), r2(t), …, rN(t)}, where r, repre­sents the three coordinates of atom i. The simula­tion proceeds with the system in a prescribed initial configuration, r3N(t0), and velocity, r3N(t0), at time t = tu As the simulation proceeds, the particles evolve through a sequence of time steps, r3N (t0) ! r3N (t1) ! r3N(t2) ! ••• ! r3N(tL), where tk = t0 + kAt,

k = 1,2,…, L, and At is the time step ofMD simulation. The simulation runs for L number of steps and covers a time interval of LAt. Typical values of L can range from 104 to 108 and At ~ 10~15s. Thus, nominal MD simulations follow the system evolution over time inter­vals not more than ~ 1-10 ns.

The simulation system has a certain energy E, the sum of the kinetic and potential energies of the particles, E = K + U where K is the sum of individual kinetic energies

1 ^

K = 2m v • v [1]

y=1

and U = U(r3N) is a prescribed interatomic interac­tion potential. Here, for simplicity, we assume that all particles have the same mass m. In principle, the potential U is a function of all the particle coordi­nates in the system if we allow each particle to interact with all the others without restriction. Thus, the dependence of U on the particle coordinates can be as complicated as the system under study demands. However, for the present discussion we introduce an approximation, the assumption of a two-body or pair-wise additive interaction, which is sufficient to illustrate the essence of MD simulation.

To find the atomic trajectories in the classical version of MD, one solves the equations governing the particle coordinates, Newton’s equations of motion in mechanics. For our N-particle system with potential energy U, the equations are

d2 r

m =-Vr;U (r3N ); j = 1; …; N [2]

where m is the particle mass. Equation [2] may look deceptively simple; actually, it is as complicated as the famous N-body problem that one generally cannot solve exactly when N is >2. As a system of coupled second-order, nonlinear ordinary differential equations, eqn [2] can be solved numerically, which is what is carried out in MD simulation.

Equation [2] describes how the system (particle coordinates) evolves over a time period from a given initial state. Suppose we divide the time period ofinter — est into many small segments, each being a time step of size At. Given the system conditions at some initial time t0, r3N(t0), and r3N(t0), integration means we advance the system successively by increments of At,

r3N (to) ! r3N (tl) ! r3N (t2) ! • • • ! r3N (tL) [3]

where L is the number of time steps making up the interval of integration.

How do we numerically integrate eqn [3] for a given U? A simple way is to write a Taylor series expansion,

r, (to + At) = r, (to) + v, (to)At

2 [4]

+ l/2ay (to)(At) + •••

and a similar expansion for ry (t0 — At). Adding the two expansions gives

r, (to + At)= — rj (to — At) + 2ry (to)

+ ay (to)(At) + •••

Notice that the left-hand side of eqn [5] is what we want, namely, the position of particle j at the next time step to + At. We already know the positions at to and the time step before, so to use eqn [5] we need the acceleration of particle y at time to. For this we substitute Fy(r3N(to))/m in place of acceleration ay(to), where Fy is just the right-hand side of eqn [2]. Thus, the integration of Newton’s equations of motion is accomplished in successive time incre­ments by applying eqn [5]. In this sense, MD can be regarded as a method of particle tracking where one follows the system evolution in discrete time steps. Although there are more elaborate, and therefore more accurate, integration procedures, it is impor­tant to note that MD results are as rigorous as classical mechanics based on the prescribed interatomic poten­tial. The particular procedure just described is called the Verlet (leapfrog)lci method. It is a symplectic integrator that respects the symplectic symmetry of the Hamiltonian dynamics; that is, in the absence of floating-point round-off errors, the discrete mapping rigorously preserves the phase space volume.1 , Symplectic integrators have the advantage of long­term stability and usually allow the use of larger time steps than nonsymplectic integrators. However, this advantage may disappear when the dynamics is not strictly Hamiltonian, such as when some thermostat — ing procedure is applied. A popular time integrator used in many early MD codes is the Gear predictor- corrector method13 (nonsymplectic) of order 5. Higher accuracy of integration allows one to take a larger value of At so as to cover a longer time interval for the same number of time steps. On the other hand, the trade-off is that one needs more computer mem­ory relative to the simpler method.

A typical flowchart for an MD code11 would look something like Figure 3. Among these steps, the part that is the most computationally demanding is the force calculation. The efficiency of an MD simulation therefore depends on performing the force calculation as simply as possible without compromising the phys­ical description (simulation fidelity). Since the force is calculated by taking the gradient of the potential U, the specification of U essentially determines the com­promise between physical fidelity and computational efficiency.

1.09.2 The Interatomic Potential

This is a large and open-ended topic with an exten­sive literature.14 It is clear from eqn [2] that the interaction potential is the most critical quantity in MD modeling and simulation; it essentially controls the numerical and algorithmic simplicity (or complex­ity) of MD simulation and, therefore, the physical fidelity of the simulation results. Since Chapter 1.10, Interatomic Potential Development is devoted to interatomic potential development, we limit our dis­cussion only to simple classical approximations to U(rb r2, . . . , rN).

Practically, all atomistic simulations are based on the Born-Oppenheimer adiabatic approximation, which separates the electronic and nuclear motions.15 Since electrons move much more quickly because of their smaller mass, during their motion one can treat the nuclei as fixed in instantaneous positions, or equivalently the electron wave functions follow the nuclear motion adiabatically. As a result, the electrons are treated as always in their ground state as the nuclei move.

For the nuclear motions, we consider an expansion of U in terms of one-body, two-body, … N-body interactions:

NN

U (r3N )= V1(rj )+ V2(ri; rj )

}=1 i<J

[6]

N

+ V3(ri; rj ; rkH———-

i<j <k

The first term, the sum of one-body interactions, is usually absent unless an external field is present to couple with each atom individually. The second sum is the contribution ofpure two-body interactions (pairwise additive). For some problems, this term alone is sufficient to be an approximation to U. The third sum represents pure three-body interactions, and so on.

Although the various glue-type potentials attribute different aspect of the physics to the N-body and pair­wise terms in the potential, if one has complete free­dom in choosing the functions for V(r), F(p), and f, then it is possible to move energy between the two terms. Johnson and Oh noted that the EAM potential

U = X Vj(rj)+f [ X fj(r’j).

is invariant under a transformation

ъ [X f (rj) ! F [X f (rj) + A X f (rj)

For an alloy,45

1 fb (r) f a (r)

Vab(r) = ‘ TVa(r) + Vb(r)

Thus, it is possible to choose a ‘gauge’ for the poten­tial, for example, by setting F'(p0) = 0 for some ref­erence density p. The advantage of the gauge transformation is that it simplifies fitting the potential. It eliminates terms in F’(p0) for pressure and elastic moduli at the equilibrium volume: these terms are nonlinear in the fitting parameters. Thus, the fitting process can be done by linear algebra.

The downside of the gauge transformation is that it destroys the physical intuition behind the form of the many-body term. Moreover, the gauge is deter­mined by a particular reference configuration, a sim­ple concept for elements, but one which does not transfer readily to alloys.

The FS potentials do not have this freedom, because the function F is predefined as a square root. However, they introduced the ‘effective pair potential’

Veff(rj) = V(r) — f (r)/v’p0

where p0 is a reference configuration (typically the equilibrium crystal structure). Many of the equilibrium properties which they used for fitting depend only on this quantity.

In addition to gauge transformation, MD depends only on the derivative of the total energy. Energy can be partitioned between atoms in any way one likes, without changing the physical results. However, on-atom prop­erties, such as the magnetic moment in magnetic
potentials, typically do depend on the partition of energy between atoms. Such quantities do not have the gauge-invariance property.

In this section we consider some other cases, such as those involved in reaction R5 ignored above, and refocus some conclusions already made on other
reactions. We note here that so far the need to inves­tigate DD in irradiated metals was led mainly by the need to create multiscale modeling tools for predict­ing changes in mechanical properties due to irradia­tion. This is desirable for practical estimations in engineering support of real nuclear devices. How­ever, there is another consequence of dislocation activity that is related to microstructure changes that occur. While this may not be important in post­mortem experiments on irradiated materials, it can be important in real devices operating under irradiation. It is obvious that internal and external stresses can accumulate during irradiation ofcomplicated devices due to high temperature, radiation growth, swelling, and transition periods of operation, for example, shut down and restart. Creep is usually taken into account but microstructure changes due to dislocation activ­ity during irradiation are not. This activity affects the whole process of microstructure evolution and should therefore be taken into account in predicting the effects of irradiation. The validity of this state­ment is demonstrated by recent in-reactor straining experiments on some bcc and fcc metals and alloys.9 It is not possible at the moment to formulate unam­biguous conclusions and more experimental work will be necessary for this. Nevertheless, it is clear that dislocation activity during irradiation directly affects the radiation damage process. In the following section, we describe some mechanisms that can con­tribute to this at the atomic-scale level.

First, consider reaction R5, omitted in Section 1.12.4.2, for dislocation obstacles. Drag of glissile interstitial loops by moving edge dislocations was first observed in the MD modeling of bcc and fcc metals.89 It is well known that SIA clusters in the form of small perfect DLs exhibit thermally activated glide in the direction of their Burgers vector.90,91 It is characterized by very low activation energy ~0.01— 0.10 eV. An edge dislocation, having a long-range

elastic field, can interact with such clusters and, if bL is parallel to the dislocation glide plane, can drag or push them as it moves under applied stress/strain. The dynamics of this process have been investigated in detail and correlations between cluster and disloca­tion mobility analysed.92-94 Additional friction due to cluster drag reduces dislocation velocity, an effect that is stronger in fcc than bcc metals because offeatures of cluster structure.89 The maximum speed at which a dislocation can drag an SIA cluster is achieved by a compromise between dislocation-cluster interaction force and cluster friction and varies at T = 300 K from 180ms-1 in fcc Cu to >1000ms-1 in bcc Fe for loops containing a few tens of SIAs.89 An important consequence of this drag process is that a moving dislocation can sweep glissile clusters and transport them through the material.

Other reactions that may affect microstructure evolution involve both inclusion — and dislocation­like obstacles. The relevant reaction for the latter is denoted R3 in Table 1. In this case, an edge disloca­tion climbs and the formation of superjogs by defect absorption changes its structure and total line length. This changes its mobility and its cross-section for interaction with other defects, such as point defects, their clusters and impurities, and this in turn affects microstructure evolution. Reactions of type R2 can also be important for they change properties of obsta­cles. In thermal aging without stress, obstacles such as voids, SIA clusters and SFTs evolve towards their equilibrium low-energy state, that is voids/precipi — tates into faceted near-spherical shapes, SFTs into regular tetrahedron shape, and so on, whereas shear­ing creates interface steps on voids/precipitates and creates ledges on SFT faces. These surface features change the properties of the defects by putting them into a higher energy state.

Reaction R4 introduces another mechanism of mass transport, for the helical turn (representing the absorbed defects) can only extend or translate in the direction of its Burgers vector, that is, along the screw dislocation line. The case of dislocation-SFT interac­tion in a thin film considered in Section 1.12.4.2.1 has demonstrated that this may introduce completely new mechanisms. This effect can also play a role when a dislocation ends on an internal interface where it can cross-slip. Reaction R4 also orders the orientation of DLs left behind by a gliding screw dislocation for it changes bL of these loops to b of the dislocation, irrespective of their initial orientation.

Finally, we describe a case when several of the above mechanisms may have a significant effect on microstructure changes if operating at the same time on different defects. It is known from experimental studies6 that under neutron irradiation, Cu accumu­lates a high density of SFTs and this density saturates with dose at a high level (^1024m~3), close to con­ditions under which displacement cascades overlap. Taking into account that SIA clusters are necessarily accumulated in the system, this high saturation den­sity implies that annihilation reactions between the vacancy population in SFTs and SIA loops is sup­pressed. MD modeling in which an SIA cluster was placed between two SFTs about 10 nm apart and intersecting the loop glide cylinder has confirmed that annihilation reaction does not occur even after 50 ns at T< 900 K.95 The result is not surprising for each vacancy of an SFT is distributed over the four faces of the tetrahedron within the stacking faults.

However, simulations show that an annihilation reaction can be promoted by the involvement of a gliding dislocation. Two cases have been considered. In one, an edge dislocation under applied stress dragged a 1/2(110) SIA loop toward an SFT placed 7 nm below the dislocation slip plane and intersecting the glide cylinder of the loop. The overlapping frac­tions of SFT and dragged SIA loop annihilated by recombination. Different obstacle sizes (SFT from 45 to 61 vacancies and SIA loops from 37 to 91 SIAs) and geometries with different levels of overlap were simulated and recombination occurred in all cases. In the other situation, the same SFT and SIA loop 10-20 nm apart were placed on the slip plane of a screw dislocation. On approaching the obstacles, the dislocation absorbed a portion of each to form two helical turns (reaction R4). The turns of vacancy and interstitial character had opposite sign and the smal­ler was annihilated by recombination with part of the larger. On moving farther ahead, the dislocation released the unrecombined portion of the remaining helix to leave a small defect. This was usually part of the original SIA loop because a loop can be completely absorbed as a helical turn, whereas only a part of an SFT can be absorbed in this way. Thus, the overall result of interstitial loop drag under applied stress was annihilation of a significant part of both clusters by a reaction between helices of opposite signs.

As discussed in the previous section, the PBM changed the concept of RDT by recognizing that qualitatively different mechanisms operate in materials when the initial damage is in the form of only FPs and under neutron irradiation, when ther­mally stable glissile SIA clusters are continuously produced in cascades. The successful applications of the PBM have been limited to low irradiation doses (<1dpa) and pure metals (e. g., Cu). Furthermore, it predicts the saturation of void size with increas­ing irradiation dose. Thus, it fails to account for the most important observation under neutron or heavy-ion irradiation: continuous increase in void swelling.

The observed continuous void growth may be explained by the development of spatial correlations between voids and other lattice defects. Such as, precipi­tates and dislocations, that shadow voids from the SIA clusters (see Figures 8-10). It has been argued that this must be the case and the very absence of a void lattice (i. e., a particular case of spatial correla­tion, which is between voids) must be an indication that spatial correlations with other defects prevail.35

To account for this effect, a new parameter, ^c, has been introduced, called the 1 correlation-screening factor, which is equal to unity in the absence of shadowing effects and zero when voids are screened completely from the SIA clusters. The swelling rate is then given by

dS

df

where F is a proportionality coefficient, which is a function of all parameters involved.35 Experimental evidence on the association of large voids with vari­ous precipitates (G, ^, Laves, etc.)120,160-162 and the compression side of edge dislocations163,164 has been available for a long time. More recent evidence can be attributed to Portnykh et al.165 who studied the microstructure of 20% cold-worked 16Cr-15Ni- 2Mo-2Mn austenitic steel irradiated up to ^100 dpa in a BN-600 fast reactor in the temperature range from 410 to 600 °C. TEM studies revealed voids of three types: a-type associated with dislocations, b-type associated with G-phase precipitates and c-type distributed homogeneously. The c-type voids were the smallest and made practically no contribution to swelling, while the a-type voids were the largest. Such spatial correlations must be a common feature under cascade irradiation.

As discussed in Barashev and Golubov,35 the experimental data on void swelling can be fit by eqn

[143] with an appropriate choice of the dependence of on the irradiation dose (see Figure 10). At high doses, the voids must be completely shielded from the SIA clusters: = 0, and the steady-state swelling

rate of ~1% per dpa observed in austenitic steels33 can be interpreted as being equal to about half of the production bias, that is, the fraction of SIAs produced as 1D mobile clusters:

dS 1 g dfNRT ~ 2 Wi

where esurv = (1 — er) « 0.1 is the survival fraction of defects in displacement cascade. The weak depen­dence on steel composition observed is probably because the final defect structure is defined by early stages of cascades, when the energies involved are much higher than the binding energies of defects with solute atoms. The observed correlation of the incubation period prior to swelling with the forma­tion of a dislocation network may be connected with an increase of the volume for the nucleation and growth of voids in which voids are screened from the SIA clusters. Higher dislocation density also cor­responds to a smaller dislocation climb rate, which might be essential for preserving void-dislocation correlations.

Another distinguishing feature of neutron irradia­tion is transmutation of atoms, which transform even pure metals into alloys with increasing irradiation dose. The atmospheres of solute (or transmuted) elements near voids may repel SIA clusters and, hence, assist or even solely explain the unlimited void growth. It was shown (see, e. g., Golubov,166 Golubov etal.,167 and references therein) that RIS can provide an additional mecha­nism of preferential absorption of mobile defects even in the framework of FP3DM, causing a ‘segre­gation’ bias, which must be different for immobile defects (e. g., voids) and mobile defects, such as dis­locations. In the PBM, the interaction of the mobile SIA clusters with different defects may even be more important. Solute atoms may also decrease the mobil­ity of SIA clusters, thereby increasing the recombi­nation rate with migrating vacancies. In the case of very high binding energy of SIA clusters with impu­rity atoms, the ‘Singh—Foreman catastrophe’18 discussed in Section 1.13.6.2.1 may occur.

Thus, two additional features beyond those already in the PBM distinguish the microstructure evolution under neutron compared to electron irradiation at high enough doses: transmutation of atoms and devel­opment of spatial correlations. A fully predictive theory must account for these effects.

The development of a predictive theory requires revisiting all its essential elements: nucleation, growth, movement of voids, and other lattice defects in the presence of spatial correlations, etc. Carefully planned experiments spanning different tempera­tures, defect production rates, etc., must be a central part of these future studies. Development of the RIS theory for accounting for the SIA clusters is necessary for understanding the sensitivity of micro­structure to material composition. Generally, the challenge is to create a theory, where the mean-field approach in its conventional form is abandoned, a task not attempted before.

Acknowledgments

The authors would like to acknowledge the fruitful collaboration and discussions on the physics of radiation damage for many years with Professor Yu. V. Konobeev (IPPE, Russia), Drs. B. N. Singh (Riso National Laboratory, Denmark), H. Trinkaus (Forschungscentrum Julich, Germany), SJ. Zinkle and Yu. N. Osetsky (Oak Ridge National Laboratory, USA), and Professor D. J. Bacon (The University of Liverpool, UK). Various aspects of the author’s research discussed in this chapter were supported by the Division of Material Science and Engineering and the Office of Fusion Energy Sciences, Depart­ment ofEnergy and the Office ofNuclear Regulatory Research, US (SIG and RES) and the UK Engineer­ing and Physical Science Research Council (AVB).