**3. Multiscale materials models**

262 Nuclear Reactors

energy neutrons. Concurrently, high concentrations of fission products (in fuels) and transmutants (in cladding and structural materials) are generated, which can cause pronounced effects in the overall chemistry of the material, especially at high burnup. The primary knock-on atoms, as well as recoiling fission products and transmutant nuclei quickly lose kinetic energy through electronic excitations (that are not generally believed to produce atomic defects) and a chain of atomic collision displacements, generating cascade of vacancy and self-interstitial defects. High-energy displacement cascades evolve over very short times, 100 picoseconds or less, and small volumes, with characteristic length scales of 50 nm or less, and are directly amenable to molecular dynamics (MD) simulations if accurate potential functions are available and chemical reactions are not occurring. If change in electronic structure need to be included, then ab initio MD is needed and this is beyond

In order to simulate the appropriate reactor conditions for all models, it is important to connect the parameters of the atomistic models with reactor conditions and the type of irradiation encountered. The radiation damage event is composed of several distinct processes concluded by a displacement cascade (collection of point defect due to the PKA) and by the formation of an interstitial –which occurs when the PKA comes to rest-. In order to simulate the radiation effects, it is important to determine the type of energetic particle interaction we wish to model. In nuclear reactors, neutrons and charged fission product particles are the dominant energetic species produced (Beta and Gamma rays are also produced, but these create less damage than the neutrons and charged particles). The type of reactor determines the nature of the dominant energetic particle interaction. The proposed study will focus on neutrons and He ions. The additional aspect to consider concerns the energy of the PKA, which in the case of D-T fusion reaction can reach the order of ~1MeV. These energies are out of reach of atomistic simulations. Nonetheless cascade event simulations at lower energies –ranging from 5to 45 KeV- can yield significant insight on the evolution of defect type and number as a function of PKA energy. Was (2007) and Olander (1981) have extensively documented how it is possible to determine the primary damage state due to irradiation by energetic particles. The simplest model is one that approximates the event as colliding hard spheres with displacement occurring when the

transferred energy is high enough to knock the struck atom off its lattice site.

The physics of primary damage production in high-energy displacement cascades has been extensively studied with MD simulations(Was 2007). The key conclusions from those MD studies of cascade evolution have been that i) intra-cascade recombination of vacancies and self-interstitial atoms (SIAs) results in ~30% of the defect production expected from displacement theory, ii) many-body collision effects produce a spatial correlation (separation) of the vacancy and SIA defects, iii) substantial clustering of the SIAs and to a lesser extent, the vacancies occurs within the cascade volume, and iv) high-energy displacement cascades tend to break up into lobes, or sub-cascades which may also enhance recombination(Calder and Bacon 1993; Calder and Bacon 1994; Bacon, Calder et al. 1995;

It is the subsequent transport and evolution of the defects produced during displacement cascades, in addition to solutes and transmutant impurities, that ultimately dictate radiation effects in materials, and changes in material microstructure(Odette et al. 2001; Wirth et al. 2004). Spatial correlations associated with the displacement cascades continue to play an important role over much larger scales, as do processes including defect recombination,

current capabilities.

Phythian, Stoller et al. 1995).

Fig. 2. Illustration of the length and time scales (and inherent feedback) involved in the multiscale processes responsible for microstructural changes in irradiated materials

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 265

where, T is the total energy transferred, Ei is the energy of the incident neutron, θ is the angle of collision, and A is the mass of the lattice nucleus. With the assumption that scattering is isotropic in the center of mass system, the average energy transferred over all angles can be shown to be the average of the minimum and maximum possible transfer energies, i.e., .

For iron, the energy required to displace an atom is about 40 eV, depending on the direction from which it is struck. So, a neutron needs a minimum energy of about 581 eV to displace an iron atom. Neutrons produced from fission of uranium carry around 2 MeV of kinetic energy and so have potential to cause damage as they slow down. Additionally, deuterium-

The first attempt to create a model for defect production based on PKA energy comes from Kinchin and Pease (1955). In this model, above a certain threshold Ed, energy is lost only to electron excitation, while below it, energy is lost only in hard-sphere elastic scattering. Norgett, Robinson et al. (1975) proposed a revised model, taking into account a more realistic energy transfer cross section, based on binary collision model simulations, . Here, ND is the number of Frenkel pairs surviving relaxation and the damage energy ED is the amount of energy available for creating displacements through elastic collisions and is a function of T. Since some of the energy of the cascade is lost to electronic excitation, ED will be less than T; for the energy range considered in this paper ED can be estimated as equal to T. This model is frequently used as the standard for estimating DPA, but many molecular dynamics simulations have shown that it tends to strongly overestimate the actual damage efficiency. Bacon, Gao et al. (2000)proposed an empirical relationship between ND and T,where A and n are weakly temperature dependent constants fit to particular materials,

The study of irradiation damage cascades has been a popular topic over the last fifteen or so years. A through literature review of the many different of damage cascade simulations, such as binary collision approximation and kinetic Monte Carlo, that have been performed in a variety of materials is beyond the scope of this paper and readers are referred to (Was 2007). The following brief review will concentrate solely on molecular dynamics simulations in α-Fe. A thorough review of results from many papers was written by Malerba (2006).

Malerba (2006) states that the first published MD study in alpha-iron was performed by Calder and co-workers (Calder and Bacon 1993; Calder and Bacon 1994). Eighty cascades with PKAs of up to 5 keV were analyzed for properties such as percent of defects surviving relaxation, channeling properties, temperature dependence, and clustering. The interatomic potential used was developed by (Finnis and Sinclair 1984) and stiffened by Calder and Bacon to treat small interatomic distance properly. This article established a large base of

Following this initial study, many papers came out which utilized both the modified FS potential mentioned above and competing multi-body potentials including those from Johnson and Oh (1989), Harrison, Voter et al. (1989), and Simonelli, Pasianot et al. (1994). These papers had three main motivations: to generate data from a new potential, to compare data between two or more potentials, or to compare damage in α-Fe with that in another material such as copper. The main difficulties in comparing results from different authors are defining what makes up a cluster of defects and non-reporting of exactly how cascades were generated. Many authors contributed to generate databases; some papers of note are described here. Stoller, Odette et al. (1997), using the FS potential modified by Calder and Bacon, ran a number

tritium fusion reactions produce neutrons with energy of 14.1 MeV.

respectively equal to 5.57 and 0.83 for Fe at 100 K, and T is in keV.

data for future papers to compare with.

A hierarchy of models is employed in the theory and simulation of complex systems in materials science and condensed matter physics: macroscale continuum mechanics, macroscale models of defect evolution, molecular scale models based on classical mechanics, and various techniques for representing quantum-mechanical effects. These models are classified according to the spatial and temporal scales that they resolve (Figure 2). In this figure, individual modeling techniques are identified within a series of linked process circles showing the overlap of relevant length and timescales. The modeling methodology includes *ab initio* electronic structure calculations, molecular dynamics (MD), accelerated molecular dynamics, kinetic Monte Carlo (KMC), phase field equations or rate theory simulations with thermodynamics and kinetics by passing information about the controlling physical mechanisms between modeling techniques over the relevant length and time scales. The key objective of such an approach is to track the fate of solutes, impurities and defects during irradiation and thereby, to predict microstructural evolution. Detailed microstructural information serves as a basis for modeling the physical behavior through meso (represented by KMC, dislocation dynamics, and phase field methods) and continuum scale models, which must be incorporated into constitutive models at the continuum finite element modeling scale to predict performance limits on both the test coupons and components.

#### **4. Modeling and simulation examples**

To span length and temporal scales, these methods can be linked into a multi-scale simulation. The objective of the lower length scale modeling is to provide constitutive properties to higher length scale, continuum level simulations, whereas these higher length scale simulations can provide boundary conditions to the lower length scale models, as well as input regarding the verification/validity of the predicted constitutive properties. Here some examples are provided of models and simulations of materials under irradiation. The emphasis is on the development of the model, the assumptions and the underlying physics that goes into model development.

#### **4.1 Atomistic simulations of radiation damage cascades**

Damage to materials caused by neutron irradiation is an inherently multiscale phenomenon. Macroscopic properties, such as plasticity, hardness, brittleness, and creep behavior, of structural reactor materials may change due to microstructural effects of radiation. Atomistic models can be a useful tool to generate data about the structure and development of defects, on length and time scales that experiments cannot probe. Data from simulations can be fed into larger scale models that predict the long term behavior of materials subject to irradiation.

The amount of energy that an incident particle can transfer to a lattice atom is a function of their masses and the angle at which the collision occurs. Energy can be lost through inelastic collisions, (n, 2n) or (n, γ) reactions, and, most importantly, elastic collisions. Elastic collisions between neutrons and nuclei can be treated within the hard sphere model with the following equations:

$$T = \frac{\mathcal{Y}}{2} E\_i (1 - \cos \theta) \tag{1}$$

$$\gamma = \frac{4A}{A+1} \tag{2}$$

A hierarchy of models is employed in the theory and simulation of complex systems in materials science and condensed matter physics: macroscale continuum mechanics, macroscale models of defect evolution, molecular scale models based on classical mechanics, and various techniques for representing quantum-mechanical effects. These models are classified according to the spatial and temporal scales that they resolve (Figure 2). In this figure, individual modeling techniques are identified within a series of linked process circles showing the overlap of relevant length and timescales. The modeling methodology includes *ab initio* electronic structure calculations, molecular dynamics (MD), accelerated molecular dynamics, kinetic Monte Carlo (KMC), phase field equations or rate theory simulations with thermodynamics and kinetics by passing information about the controlling physical mechanisms between modeling techniques over the relevant length and time scales. The key objective of such an approach is to track the fate of solutes, impurities and defects during irradiation and thereby, to predict microstructural evolution. Detailed microstructural information serves as a basis for modeling the physical behavior through meso (represented by KMC, dislocation dynamics, and phase field methods) and continuum scale models, which must be incorporated into constitutive models at the continuum finite element modeling scale to predict performance limits on both the test coupons and components.

To span length and temporal scales, these methods can be linked into a multi-scale simulation. The objective of the lower length scale modeling is to provide constitutive properties to higher length scale, continuum level simulations, whereas these higher length scale simulations can provide boundary conditions to the lower length scale models, as well as input regarding the verification/validity of the predicted constitutive properties. Here some examples are provided of models and simulations of materials under irradiation. The emphasis is on the development of the model, the assumptions and the underlying physics

Damage to materials caused by neutron irradiation is an inherently multiscale phenomenon. Macroscopic properties, such as plasticity, hardness, brittleness, and creep behavior, of structural reactor materials may change due to microstructural effects of radiation. Atomistic models can be a useful tool to generate data about the structure and development of defects, on length and time scales that experiments cannot probe. Data from simulations can be fed into larger scale models that predict the long term behavior of materials subject to irradiation. The amount of energy that an incident particle can transfer to a lattice atom is a function of their masses and the angle at which the collision occurs. Energy can be lost through inelastic collisions, (n, 2n) or (n, γ) reactions, and, most importantly, elastic collisions. Elastic collisions between neutrons and nuclei can be treated within the hard sphere model with the

> (1 ) <sup>2</sup> *T E cos <sup>i</sup>* γ

> > 4 1 *A A* γ

θ

(1)

<sup>=</sup> + (2)

= −

**4. Modeling and simulation examples** 

**4.1 Atomistic simulations of radiation damage cascades** 

that goes into model development.

following equations:

where, T is the total energy transferred, Ei is the energy of the incident neutron, θ is the angle of collision, and A is the mass of the lattice nucleus. With the assumption that scattering is isotropic in the center of mass system, the average energy transferred over all angles can be shown to be the average of the minimum and maximum possible transfer energies, i.e., .

For iron, the energy required to displace an atom is about 40 eV, depending on the direction from which it is struck. So, a neutron needs a minimum energy of about 581 eV to displace an iron atom. Neutrons produced from fission of uranium carry around 2 MeV of kinetic energy and so have potential to cause damage as they slow down. Additionally, deuteriumtritium fusion reactions produce neutrons with energy of 14.1 MeV.

The first attempt to create a model for defect production based on PKA energy comes from Kinchin and Pease (1955). In this model, above a certain threshold Ed, energy is lost only to electron excitation, while below it, energy is lost only in hard-sphere elastic scattering. Norgett, Robinson et al. (1975) proposed a revised model, taking into account a more realistic energy transfer cross section, based on binary collision model simulations, . Here, ND is the number of Frenkel pairs surviving relaxation and the damage energy ED is the amount of energy available for creating displacements through elastic collisions and is a function of T. Since some of the energy of the cascade is lost to electronic excitation, ED will be less than T; for the energy range considered in this paper ED can be estimated as equal to T. This model is frequently used as the standard for estimating DPA, but many molecular dynamics simulations have shown that it tends to strongly overestimate the actual damage efficiency. Bacon, Gao et al. (2000)proposed an empirical relationship between ND and T,where A and n are weakly temperature dependent constants fit to particular materials, respectively equal to 5.57 and 0.83 for Fe at 100 K, and T is in keV.

The study of irradiation damage cascades has been a popular topic over the last fifteen or so years. A through literature review of the many different of damage cascade simulations, such as binary collision approximation and kinetic Monte Carlo, that have been performed in a variety of materials is beyond the scope of this paper and readers are referred to (Was 2007). The following brief review will concentrate solely on molecular dynamics simulations in α-Fe. A thorough review of results from many papers was written by Malerba (2006).

Malerba (2006) states that the first published MD study in alpha-iron was performed by Calder and co-workers (Calder and Bacon 1993; Calder and Bacon 1994). Eighty cascades with PKAs of up to 5 keV were analyzed for properties such as percent of defects surviving relaxation, channeling properties, temperature dependence, and clustering. The interatomic potential used was developed by (Finnis and Sinclair 1984) and stiffened by Calder and Bacon to treat small interatomic distance properly. This article established a large base of data for future papers to compare with.

Following this initial study, many papers came out which utilized both the modified FS potential mentioned above and competing multi-body potentials including those from Johnson and Oh (1989), Harrison, Voter et al. (1989), and Simonelli, Pasianot et al. (1994). These papers had three main motivations: to generate data from a new potential, to compare data between two or more potentials, or to compare damage in α-Fe with that in another material such as copper. The main difficulties in comparing results from different authors are defining what makes up a cluster of defects and non-reporting of exactly how cascades were generated.

Many authors contributed to generate databases; some papers of note are described here. Stoller, Odette et al. (1997), using the FS potential modified by Calder and Bacon, ran a number

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 267

The atomistic simulations provide information on the number of surviving defects after initial damage and are not able to simulate large time scale or length scales. They provide a good atomistic picture of the unit processes affecting the formation of defects and the evolution of the primary damage state. Experiments cannot yet access this small time and length scale of radiation damage processes, therefore, experimental corroboration is hard to find for such simulations. Results of atomistic simulations of radiation cascades can be used to develop higher length and time scale theories and simulation of radiation effects in materials. Parameters that can be passed to other simulations/theories include the number and spatial distribution of defects created at the conclusion of the radiation cascade phase.

While atomistic methods can probe the primary damage state with great detail, they can also be used to probe the interactions of the defects formed with the underlying microstructure. An example is the case of creep due to irradiation in materials. Creep of metals and alloys under irradiation has been the subject of many experimental and theoretical studies for more than 30 years. Although a vast amount of knowledge of irradiation creep has accumulated, the database on irradiation creep comes from many relatively small experiments, and there were often differences in experimental conditions from one study to the next. Theoretical models are based on linear elasticity. Among the many theories that exist to describe the driving force for irradiation creep, the most

Stress Induced Preferential Nucleation of loops (SIPN) is based on the idea that the application of external stress will result in an increased number of dislocation loops nucleating on planes of preferred orientations. Interstitial loops will tend to be oriented perpendicular to the applied tensile stress, while vacancy loops will prefer to be oriented parallel to the stress. The net result is elongation of the solid in the direction of applied stress. While there is some experimental support of this theory, it is thought that it cannot

An alternative theory is Stress Induced Preferential Absorption/Attraction (SIPA). The essential idea behind SIPA is that interstitials are preferentially absorbed by dislocations of particular orientations, resulting in climb; this is described by an elastic interaction between the stress fields of the defect and dislocation. A variant on SIPA that accounts for anistropic diffusion is SIPA-AD. This theory uses the full diffusion equations, derived by Dederichs and Schroeder (1978), to take into account anisotropic stress fields. Savino and Tome developed this theory and found that it generally gives a larger contribution to dislocation climb than the original SIPA (Tome, Cecatto et al.). A thorough review of many dislocation

These models go a long way towards explaining irradiation creep due to dislocations. However, all models based on linear elasticity break down near a dislocation core due to the 1/r terms in the stress and strain field expressions. Atomistic calculations do not suffer from this problem, so they can be used to verify the range of validity of theoretical expressions

Molecular statics calculations can be performed in order to understand the interactions between vacancies and interstitials and line dislocations in bcc iron. These can be compared

**4.2 Molecular dynamics calculations of dislocation defect interactions** 

important are the SIPN, SIPA, and SIPA-AD effects.

creep models was prepared by Matthews and Finnis (1988).

and successfully predict true behavior at the core.

account fully for creep seen in materials.

of cascades at energies up to 40 keV. They found evidence for vacancy clustering, a feature not seen in previous works. Bacon et al. (2000) performed a study comparing the cascade characteristics of bcc, hcp, and fcc metals. They found that there were no major differences in interstitial and vacancy production, so concluded that any differences observed experimentally must be due to evolution of the microstructure following the primary damage event. Caturla, Soneda et al. (2000) compared bcc Fe with fcc Cu, finding that clustering in Fe was at least an order of magnitude less than in Cu. Terentyev, Lagerstedt et al. (2006) produced a study looking solely at differences between four available potentials by applying the same defect counting criteria to each. They found that the stiffness of a potential, a somewhat arbitrary feature, was the most important factor in determining cascade properties.

In all primary damage cascade simulations, first, the incident radiation has an interaction with an atom in the crystal lattice, transferring enough energy to remove the atom from its site. This atom, the primary knock-on atom (PKA), goes on to interact with other atoms in the crystal, removing them from their sites and generating a displacement cascade in the thermal spike phase. Atoms that are removed from their perfect lattice sites and come to rest between other atoms are known as interstitials; the empty lattice sites they leave are called vacancies. At some time shortly after the PKA is created, some peak number of Frenkel pairs, Np, will exist in the crystal, where a Frenkel pair is defined as one vacancy plus one interstitial. After this point, the defects will begin to recombine as the energy is dissipated. After a few picoseconds, only a few defects, Nd will remain. This generally results in a core of vacancies surrounded by a shell of interstitials. A profile of the number of defects over time in a typical cascade can be seen in Figure 3.

Fig. 3. Initial stages of radiation damage cascade, the number of vacancies and interstitials as a functions of time (from Hayward and Deo 2010)

of cascades at energies up to 40 keV. They found evidence for vacancy clustering, a feature not seen in previous works. Bacon et al. (2000) performed a study comparing the cascade characteristics of bcc, hcp, and fcc metals. They found that there were no major differences in interstitial and vacancy production, so concluded that any differences observed experimentally must be due to evolution of the microstructure following the primary damage event. Caturla, Soneda et al. (2000) compared bcc Fe with fcc Cu, finding that clustering in Fe was at least an order of magnitude less than in Cu. Terentyev, Lagerstedt et al. (2006) produced a study looking solely at differences between four available potentials by applying the same defect counting criteria to each. They found that the stiffness of a potential, a somewhat arbitrary feature, was the most important factor in determining cascade properties. In all primary damage cascade simulations, first, the incident radiation has an interaction with an atom in the crystal lattice, transferring enough energy to remove the atom from its site. This atom, the primary knock-on atom (PKA), goes on to interact with other atoms in the crystal, removing them from their sites and generating a displacement cascade in the thermal spike phase. Atoms that are removed from their perfect lattice sites and come to rest between other atoms are known as interstitials; the empty lattice sites they leave are called vacancies. At some time shortly after the PKA is created, some peak number of Frenkel pairs, Np, will exist in the crystal, where a Frenkel pair is defined as one vacancy plus one interstitial. After this point, the defects will begin to recombine as the energy is dissipated. After a few picoseconds, only a few defects, Nd will remain. This generally results in a core of vacancies surrounded by a shell of interstitials. A profile of the number of defects over

Fig. 3. Initial stages of radiation damage cascade, the number of vacancies and interstitials as

time in a typical cascade can be seen in Figure 3.

a functions of time (from Hayward and Deo 2010)

The atomistic simulations provide information on the number of surviving defects after initial damage and are not able to simulate large time scale or length scales. They provide a good atomistic picture of the unit processes affecting the formation of defects and the evolution of the primary damage state. Experiments cannot yet access this small time and length scale of radiation damage processes, therefore, experimental corroboration is hard to find for such simulations. Results of atomistic simulations of radiation cascades can be used to develop higher length and time scale theories and simulation of radiation effects in materials. Parameters that can be passed to other simulations/theories include the number and spatial distribution of defects created at the conclusion of the radiation cascade phase.

### **4.2 Molecular dynamics calculations of dislocation defect interactions**

While atomistic methods can probe the primary damage state with great detail, they can also be used to probe the interactions of the defects formed with the underlying microstructure. An example is the case of creep due to irradiation in materials. Creep of metals and alloys under irradiation has been the subject of many experimental and theoretical studies for more than 30 years. Although a vast amount of knowledge of irradiation creep has accumulated, the database on irradiation creep comes from many relatively small experiments, and there were often differences in experimental conditions from one study to the next. Theoretical models are based on linear elasticity. Among the many theories that exist to describe the driving force for irradiation creep, the most important are the SIPN, SIPA, and SIPA-AD effects.

Stress Induced Preferential Nucleation of loops (SIPN) is based on the idea that the application of external stress will result in an increased number of dislocation loops nucleating on planes of preferred orientations. Interstitial loops will tend to be oriented perpendicular to the applied tensile stress, while vacancy loops will prefer to be oriented parallel to the stress. The net result is elongation of the solid in the direction of applied stress. While there is some experimental support of this theory, it is thought that it cannot account fully for creep seen in materials.

An alternative theory is Stress Induced Preferential Absorption/Attraction (SIPA). The essential idea behind SIPA is that interstitials are preferentially absorbed by dislocations of particular orientations, resulting in climb; this is described by an elastic interaction between the stress fields of the defect and dislocation. A variant on SIPA that accounts for anistropic diffusion is SIPA-AD. This theory uses the full diffusion equations, derived by Dederichs and Schroeder (1978), to take into account anisotropic stress fields. Savino and Tome developed this theory and found that it generally gives a larger contribution to dislocation climb than the original SIPA (Tome, Cecatto et al.). A thorough review of many dislocation creep models was prepared by Matthews and Finnis (1988).

These models go a long way towards explaining irradiation creep due to dislocations. However, all models based on linear elasticity break down near a dislocation core due to the 1/r terms in the stress and strain field expressions. Atomistic calculations do not suffer from this problem, so they can be used to verify the range of validity of theoretical expressions and successfully predict true behavior at the core.

Molecular statics calculations can be performed in order to understand the interactions between vacancies and interstitials and line dislocations in bcc iron. These can be compared

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 269

Fig. 4c. Intermediate position of interstitial near dislocation creates a jog in the dislocation

Fig. 4d. Final position of the dislocation with the interstitial absorbed in the core

Figures 4a-d shows the evolution of a defect in the vicinity of the dislocation. In this example an interstitial is positioned near a dislocation core and the energy of the system is

to similar results given by dipole tensor calculations based in linear elasticity theory. Results from two methods are used to calculate the interaction energy between a dislocation and a point defects in bcc iron are compared. For vacancies and a variety of self-interstitial dumbbell configurations near both edge and screw dislocation cores, there are significant differences between direct calculations and atomistics. For vacancies some interaction is seen with both edge and screw dislocations where none is predicted. Results for interstitials tended to have a strong dependence on orientation and position about the core. Particularly for the screw, continuum theory misses the tri-fold splitting of the dislocation core which has a large influence on atomistic results.


Fig. 4a. Initial position of the interstitial near the dislocation core

Fig. 4b. Intermediate position of the interstitial as it folds into the dislocation

to similar results given by dipole tensor calculations based in linear elasticity theory. Results from two methods are used to calculate the interaction energy between a dislocation and a point defects in bcc iron are compared. For vacancies and a variety of self-interstitial dumbbell configurations near both edge and screw dislocation cores, there are significant differences between direct calculations and atomistics. For vacancies some interaction is seen with both edge and screw dislocations where none is predicted. Results for interstitials tended to have a strong dependence on orientation and position about the core. Particularly for the screw, continuum theory misses the tri-fold splitting of the dislocation core which

has a large influence on atomistic results.

Fig. 4a. Initial position of the interstitial near the dislocation core

Fig. 4b. Intermediate position of the interstitial as it folds into the dislocation

Fig. 4c. Intermediate position of interstitial near dislocation creates a jog in the dislocation

Fig. 4d. Final position of the dislocation with the interstitial absorbed in the core

Figures 4a-d shows the evolution of a defect in the vicinity of the dislocation. In this example an interstitial is positioned near a dislocation core and the energy of the system is

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 271

Fig. 5. Shows the basic mechanisms of helium and vacancy activity in single crystal bcc iron. Large filled circles represent iron, large open circles represent vacancies, small filled circles represent helium atoms and small open circles represent the octahedral bcc sites. (a) Helium migration on the octahedral sublattice (b) Vacancy and self interstitial migration in bcc iron (c) Dissociation of helium from an embryonic bubble (d) dissociation of vacancy from an

lowest-energy migration path of the SIA corresponds to a nearest-neighbor translationrotation jump in the <111> direction. The rates of migration of the point defect entities are

exp

*i i migration*

where the superscript *i* refers to the helium, self interstitial atoms and the vacancy point

constant and the temperature respectively. Two point defect entities are considered to be in a cluster when the distance between them is less than a0, which is the lattice constant of bcc iron. Interstitial atoms are not allowed to dissociate from clusters. Dissociation of the helium and the vacancy from the cluster is described in Figures 2(c) and 2(d) respectively. The rate

*migration migration*

*r = ν*

defect entities. The rate of migration of the point defect entity is *<sup>i</sup>*

*i*

*E*

*migration v* , the migration barrier is *<sup>i</sup> Emigration* , while *kB* and *T* are the Boltzmann

*B*

*k T* <sup>−</sup> 

, (3)

*migration r* , the attempt

embryonic bubble.

calculated as

frequency is *<sup>i</sup>*

minimized. The interstitial moves into the dislocation core and forms an extended jog in the dislocation core structure. Such a process is not captured by linear elasticity calculations which fail to capture the core structure and the core-defect interactions.

The relationship between crystal plasticity and dislocation behavior in materials has motivated a wide range of experimental and computational studies of dislocation behavior (Vitek 1976; Osetsky, Bacon et al. 1999). Computational studies of dislocation activity can be performed at several different length and time scales. In some cases, a multi-scale modeling approach is adopted(Ghoniem et al. 2003). Core properties and atomic mechanisms are simulated using first principles calculations and molecular dynamics simulations. These results can be used to form the rules that govern large-scale Dislocation Dynamics (DD) simulations (Wen et al. 2005) that account for the activity of a large number of dislocation segments. Polycrystalline plasticity models are then developed that utilize the information at the atomistic scale to parameterize partial differential equations of rate dependent viscoplasticity(Deo, Tom et al. 2008). While understanding dislocation creep processes, dislocation climb rates and hence, the interaction of dislocations with point defects is an important quantity to be calculated. Here, we show how the dislocation core affects the interaction energy between the dislocation and the point defect using both linear elasticity as well as atomistic calculations.

#### **4.3 Kinetic Monte Carlo simulations of defect evolution**

Kinetic Monte Carlo models and simulations have been employed to study cascade ageing and defect accumulation at low doses using the input from atomistic simulations of cascades and defect migration properties (Barashev, Bacon et al. 1999; Gao, Bacon et al. 1999; Barashev, Bacon et al. 2000; Heinisch, Singh et al. 2000; Heinisch and Singh 2003; Domain, Becquart et al. 2004; Becquart, Domain et al. 2005; Caturla, Soneda et al. 2006). These have mostly focused on intrinsic defect (interstitial and vacancy) migration, clustering and annihilation under irradiation conditions. In this example, we focus on extrinsic gas atoms such as helium or hydrogen migration in irradiated materials.

For the Fe-He system, modeling of helium clusters has been performed by Morishita, Wirth and co-workers (Morishita, Sugano et al. 2003a; Morishita, Sugano et al. 2003b) and Bringa, Wirth et al. (2003)using semi-empirical potentials. In one paper, Morishita, Sugano et al. (2003a) performed molecular dynamics (MD) calculations to evaluate the thermal stability of helium–vacancy clusters in Fe In another paper, Morishita, Sugano et al. (2003b)have looked at dissolution of helium-vacancy clusters as a function of temperature increase using the empirical potentials for the Fe-He system. Wirth and Bringa (Wirth and Bringa 2004) have simulated the motion of one single 2He-3Vac cluster at 1000 K using the same potential system. First principles calculations of helium atoms in interstitial and substitutional sites has been performed by Fu and Willaime 2005 and can be used to provide an input parameter set for kinetic Monte Carlo calculations.

The kinetic Monte Carlo model(Deo, Srivilliputhur et al. 2006; Deo, Okuniewski et al. 2007; Deo, Srinivasan et al. 2007) consists of helium interstitials on the octahedral sublattice and vacancies on the bcc iron lattice. The migration of the free (not clustered) helium is shown in figure 5(a) and that of the vacancies and self interstitial atoms is shown in Figure 5(b). The

minimized. The interstitial moves into the dislocation core and forms an extended jog in the dislocation core structure. Such a process is not captured by linear elasticity calculations

The relationship between crystal plasticity and dislocation behavior in materials has motivated a wide range of experimental and computational studies of dislocation behavior (Vitek 1976; Osetsky, Bacon et al. 1999). Computational studies of dislocation activity can be performed at several different length and time scales. In some cases, a multi-scale modeling approach is adopted(Ghoniem et al. 2003). Core properties and atomic mechanisms are simulated using first principles calculations and molecular dynamics simulations. These results can be used to form the rules that govern large-scale Dislocation Dynamics (DD) simulations (Wen et al. 2005) that account for the activity of a large number of dislocation segments. Polycrystalline plasticity models are then developed that utilize the information at the atomistic scale to parameterize partial differential equations of rate dependent viscoplasticity(Deo, Tom et al. 2008). While understanding dislocation creep processes, dislocation climb rates and hence, the interaction of dislocations with point defects is an important quantity to be calculated. Here, we show how the dislocation core affects the interaction energy between the dislocation and the point defect using both linear elasticity

Kinetic Monte Carlo models and simulations have been employed to study cascade ageing and defect accumulation at low doses using the input from atomistic simulations of cascades and defect migration properties (Barashev, Bacon et al. 1999; Gao, Bacon et al. 1999; Barashev, Bacon et al. 2000; Heinisch, Singh et al. 2000; Heinisch and Singh 2003; Domain, Becquart et al. 2004; Becquart, Domain et al. 2005; Caturla, Soneda et al. 2006). These have mostly focused on intrinsic defect (interstitial and vacancy) migration, clustering and annihilation under irradiation conditions. In this example, we focus on extrinsic gas atoms

For the Fe-He system, modeling of helium clusters has been performed by Morishita, Wirth and co-workers (Morishita, Sugano et al. 2003a; Morishita, Sugano et al. 2003b) and Bringa, Wirth et al. (2003)using semi-empirical potentials. In one paper, Morishita, Sugano et al. (2003a) performed molecular dynamics (MD) calculations to evaluate the thermal stability of helium–vacancy clusters in Fe In another paper, Morishita, Sugano et al. (2003b)have looked at dissolution of helium-vacancy clusters as a function of temperature increase using the empirical potentials for the Fe-He system. Wirth and Bringa (Wirth and Bringa 2004) have simulated the motion of one single 2He-3Vac cluster at 1000 K using the same potential system. First principles calculations of helium atoms in interstitial and substitutional sites has been performed by Fu and Willaime 2005 and can be used to provide an input

The kinetic Monte Carlo model(Deo, Srivilliputhur et al. 2006; Deo, Okuniewski et al. 2007; Deo, Srinivasan et al. 2007) consists of helium interstitials on the octahedral sublattice and vacancies on the bcc iron lattice. The migration of the free (not clustered) helium is shown in figure 5(a) and that of the vacancies and self interstitial atoms is shown in Figure 5(b). The

which fail to capture the core structure and the core-defect interactions.

as well as atomistic calculations.

**4.3 Kinetic Monte Carlo simulations of defect evolution** 

such as helium or hydrogen migration in irradiated materials.

parameter set for kinetic Monte Carlo calculations.

Fig. 5. Shows the basic mechanisms of helium and vacancy activity in single crystal bcc iron. Large filled circles represent iron, large open circles represent vacancies, small filled circles represent helium atoms and small open circles represent the octahedral bcc sites. (a) Helium migration on the octahedral sublattice (b) Vacancy and self interstitial migration in bcc iron (c) Dissociation of helium from an embryonic bubble (d) dissociation of vacancy from an embryonic bubble.

lowest-energy migration path of the SIA corresponds to a nearest-neighbor translationrotation jump in the <111> direction. The rates of migration of the point defect entities are calculated as

$$\boldsymbol{\nu}\_{migration}^{i} = \boldsymbol{\nu}\_{migration}^{i} \exp\left(-\frac{\boldsymbol{E}\_{migration}^{i}}{k\_B T}\right) \tag{3}$$

where the superscript *i* refers to the helium, self interstitial atoms and the vacancy point defect entities. The rate of migration of the point defect entity is *<sup>i</sup> migration r* , the attempt frequency is *<sup>i</sup> migration v* , the migration barrier is *<sup>i</sup> Emigration* , while *kB* and *T* are the Boltzmann constant and the temperature respectively. Two point defect entities are considered to be in a cluster when the distance between them is less than a0, which is the lattice constant of bcc iron. Interstitial atoms are not allowed to dissociate from clusters. Dissociation of the helium and the vacancy from the cluster is described in Figures 2(c) and 2(d) respectively. The rate

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 273

and clustered vacancies are mobile. Vacancies can dissociate from a vacancy cluster as well as a helium-vacancy cluster. Helium atoms migrate on the octahedral sublattice as well as part of helium-vacancy bubbles. A substitutional helium is considered as a 1-1 He1V1 and is mobile. If a bubble has a helium-vacancy ratio greater than 5, it emits a self interstitial atom. Small bubbles migrate according to Equation 3 and Table 1 while large bubbles migrate according to Eq. 3. Self interstitial atoms and vacancies annihilate when they meet either as point defects or in a cluster. The boundary of the simulation cell acts as a sink for the point

**Entity Event E (eV)** *v***0 (s-1) Remarks** 

Vacancy Migration 0.65 6e13 Vacancy migration on the

He-V clusters 3D migration 1.1 1e14 Clusters containing up to 3

He-V clusters 3D migration --- --- Diffusivity of clusters containing

Table 1. A table of the events included in the kinetic Monte Carlo model. Migration energies

At each kMC step, the system is monitored to identify a clustering event. When two point defect entities (helium-helium, vacancy-vacancy, helium-vacancy) are in a cluster the simulation creates a mapping between the entities and the cluster such that for each cluster there are at least two entities associated with the cluster. The event catalog is updated with the new rates of event occurrence and the transition probabilities for the next kMC event are

Helium Migration 0.078 1e14 Helium migrates on the interstitial

SIA Migration 0.3 6e13 SIA migration on the substitutional

sublattice

sublattice

2.0 6e13 Vacancy dissociation from the

0.30 1e14 Helium dissociation from the

0.20 6e13 Vacancy dissociation from the vacancy cluster

to this rate

1D migration 0.1 6e13 Interstitial clusters up to size 4 are

cluster

2.0 1e14 Helium dissociation from the He-V

substitutional sublattice

helium-vacancy cluster

helium-helium cluster

considered mobile

mechanisms (Eq. 3)

vacancies atoms migrate according

more than 3 vacancies calculated by considering surface diffusion

defect entities.

Helium Dissociation

Vacancy Dissociation

Helium Dissociation

Vacancy Dissociation

Interstitial clusters

from HenVm

from HenVm

from Hen

from Vm

and attempt frequencies are provided where applicable.

calculated using Equations (4-6) using the parameters from Table 1.

of dissociation of a point defect entity (*i* = helium or vacancy) from a cluster into the bulk lattice is considered to be thermally activated and is calculated as:

$$\sigma\_{\text{dissociation}}^i = \upsilon\_{\text{dissociation}}^i \exp\left(-\frac{E\_{\text{dissociation}}^i}{k\_B T}\right) \tag{4}$$

where *<sup>i</sup> dissociation r* is the rate or dissociation, *<sup>i</sup> dissocation v* is the attempt frequencies, *<sup>i</sup> Edissociation* is the energy of dissociation. The dissociation energy *<sup>i</sup> Edissociation* of a point defect from a cluster is taken to be the sum of the energy to bind a point defect entity to the cluster and *<sup>i</sup> Emigration* . Small bubbles migrate with a Arrhenius migration rate parameterized by Table 1. Larger bubbles migrate by surface diffusion at the bubble-matrix interface (Cottrell 2003; Evans 2004),

$$D\_B = D\_S \left(\frac{3\Omega^{4/3}}{2\pi r^4}\right) \tag{5}$$

where DB and Ds are the bubble and surface diffusivities respectively, Ω is the atomic volume and r the radius of the bubble.

Morishita et al. (Morishita et al. 2003) and Fu et al (Fu and Willaime 2005) have calculated the migration energies of helium and vacancies as well as the binding energies of some heliumvacancy clusters. We employ these migration barriers to calculate the rate of migration of the point defect entities. Parameters used to calculate these quantities are described in Table 1. These parameters are used to calculate the rates migration and dissociation events (Equations 1-3) in the system and build the event catalog for the kMC simulation.

The event catalog is generated by calculating the rates of migration or dissociation of the point defect entities using Equations (3-6). The kMC event catalog consists of the migration, clustering and dissociation of the point defect entities, helium, self interstitial atoms and vacancies. The transition probability of each event is proportional to the rate of event occurrence, calculated by the Equations (3-6). We follow the well established kMC simulation algorithm (Bortz, Kalos et al. 1975; Fichthorn and Weinberg 1991) which is a stochastic, atomicscale method to simulate the time-evolution of defects and nano/microstructural evolution that focuses on individual defects and not on atomic vibrations.

Reaction pathways available in the system are tabulated in table 1. Rates of migration of events are calculated at each kMC step. Parameters are obtained both from literature using first principles calculations (Fu, Willaime et al. 2004; Fu and Willaime 2005) and also from molecular statics calculations using semi-empirical potentials as employed by Morishita et al (Morishita et al. 2003) and Wirth and Bringa (Wirth and Bringa 2004) (using the Ackland Finnis–Sinclair potential, the Wilson–Johnson potential and the Ziegler–Biersack–Littmark– Beck potential for describing the interactions of Fe–Fe, Fe–He and He–He, respectively). Helium atoms are introduced at random positions at the beginning of the simulation.

Self interstitial atoms (SIA) are produced in the simulation along with vacancies as Frenkel pairs as cascade debris. Self interstitial atoms are mobile and cluster. Self interstitial clusters up to size 5 migrate one dimensionally. Self interstitials of higher size are stationary. Mono-

of dissociation of a point defect entity (*i* = helium or vacancy) from a cluster into the bulk

*i i dissociation*

the energy of dissociation. The dissociation energy *<sup>i</sup> Edissociation* of a point defect from a cluster is taken to be the sum of the energy to bind a point defect entity to the cluster and *<sup>i</sup> Emigration* . Small bubbles migrate with a Arrhenius migration rate parameterized by Table 1. Larger bubbles migrate by surface diffusion at the bubble-matrix interface (Cottrell 2003; Evans

exp

4/3 4

3 2

where DB and Ds are the bubble and surface diffusivities respectively, Ω is the atomic

Morishita et al. (Morishita et al. 2003) and Fu et al (Fu and Willaime 2005) have calculated the migration energies of helium and vacancies as well as the binding energies of some heliumvacancy clusters. We employ these migration barriers to calculate the rate of migration of the point defect entities. Parameters used to calculate these quantities are described in Table 1. These parameters are used to calculate the rates migration and dissociation events (Equations

The event catalog is generated by calculating the rates of migration or dissociation of the point defect entities using Equations (3-6). The kMC event catalog consists of the migration, clustering and dissociation of the point defect entities, helium, self interstitial atoms and vacancies. The transition probability of each event is proportional to the rate of event occurrence, calculated by the Equations (3-6). We follow the well established kMC simulation algorithm (Bortz, Kalos et al. 1975; Fichthorn and Weinberg 1991) which is a stochastic, atomicscale method to simulate the time-evolution of defects and nano/microstructural evolution

Reaction pathways available in the system are tabulated in table 1. Rates of migration of events are calculated at each kMC step. Parameters are obtained both from literature using first principles calculations (Fu, Willaime et al. 2004; Fu and Willaime 2005) and also from molecular statics calculations using semi-empirical potentials as employed by Morishita et al (Morishita et al. 2003) and Wirth and Bringa (Wirth and Bringa 2004) (using the Ackland Finnis–Sinclair potential, the Wilson–Johnson potential and the Ziegler–Biersack–Littmark– Beck potential for describing the interactions of Fe–Fe, Fe–He and He–He, respectively).

Helium atoms are introduced at random positions at the beginning of the simulation.

Self interstitial atoms (SIA) are produced in the simulation along with vacancies as Frenkel pairs as cascade debris. Self interstitial atoms are mobile and cluster. Self interstitial clusters up to size 5 migrate one dimensionally. Self interstitials of higher size are stationary. Mono-

π*r* <sup>Ω</sup> <sup>=</sup> 

*D D B S*

*i*

*E*

*B*

*dissocation v* is the attempt frequencies, *<sup>i</sup> Edissociation* is

, (5)

*k T* − 

(4)

lattice is considered to be thermally activated and is calculated as:

*r =v*

1-3) in the system and build the event catalog for the kMC simulation.

that focuses on individual defects and not on atomic vibrations.

*dissociation r* is the rate or dissociation, *<sup>i</sup>*

volume and r the radius of the bubble.

where *<sup>i</sup>*

2004),

*dissociation dissociation*

and clustered vacancies are mobile. Vacancies can dissociate from a vacancy cluster as well as a helium-vacancy cluster. Helium atoms migrate on the octahedral sublattice as well as part of helium-vacancy bubbles. A substitutional helium is considered as a 1-1 He1V1 and is mobile. If a bubble has a helium-vacancy ratio greater than 5, it emits a self interstitial atom. Small bubbles migrate according to Equation 3 and Table 1 while large bubbles migrate according to Eq. 3. Self interstitial atoms and vacancies annihilate when they meet either as point defects or in a cluster. The boundary of the simulation cell acts as a sink for the point defect entities.


Table 1. A table of the events included in the kinetic Monte Carlo model. Migration energies and attempt frequencies are provided where applicable.

At each kMC step, the system is monitored to identify a clustering event. When two point defect entities (helium-helium, vacancy-vacancy, helium-vacancy) are in a cluster the simulation creates a mapping between the entities and the cluster such that for each cluster there are at least two entities associated with the cluster. The event catalog is updated with the new rates of event occurrence and the transition probabilities for the next kMC event are calculated using Equations (4-6) using the parameters from Table 1.

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 275

This example demonstrates the use of the kinetic Monte Carlo method to simulate defect diffusion in irradiated materials. The method depends on the ability to calculate rates of migration events of individual defects. The advantage of the method is that all atoms do not necessarily need to be simulated and large time scales are accessible to the simulation

Fig. 6. is a plot of the concentration of bubbles (cB) as a function of helium concentration (expressed in terms of appm/dpa) for two damage levels, 0.1 dpa (filled triangles) and 1 dpa (open squares). The line is a fit to the data assuming a square root dependence of the bubble density on helium concentration. The simulation temperature is 543K (0.3 Tm, the melting

The previous examples have looked at the atomistic and mesoscale of radiation damage and defect formation. This information can be used by plasticity models and microstructural mechanics models of the effect of radiation on materials properties. Here an example is presented where the atomistic calculations are used to parameterize a viscoplasticity

Hardening and embrittlement are controlled by interactions between dislocations and irradiation induced defect clusters. Radiation hardening and embrittlement that occurs in metals irradiated at low temperatures (below ~0.3 *Tm*, where *Tm* is the melting temperature) is a an important technical challenge for advanced nuclear energy systems(Zinkle and Matsukawa 2004). In this example, the Visco Plastic Self Consistent (VPSC) polycrystalline code (Lebensohn and Tome 1993) is employed in order to model the yield stress dependence

**4.4 Microstructural mechanics of irradiation hardening** 

treatment of hardening in materials due to irradiation.

provided the underlying physics remains invariant.

temperature for iron).

Simulations were performed for a damage energy of 100keV. A production rate of randomly distributed cascades is assumed such that the damage is introduced at a rate of 10-6 dpa/s. The simulation cell size is 400a0x400a0x400a0, where a0 is the lattice parameter of iron. The number of Frenkel pairs introduced in the simulation cell are calculated by the Norgett-Robinson-Torrens relationship,

$$\text{displacements per cascade} = \frac{0.8E\_D}{2E\_d} \tag{6}$$

where ED is the energy of incident particles while Ed is the threshold energy (40eV for iron (1994)). Initial concentration of helium is a parameter in the simulation and is varied from 1 to 25 appm/dpa. Damage measured in dpa is another parameter in the system and determines the length of the simulation. The incident energy is also a simulation control parameter that determines the number of defects introduced in a single cascade. The simulation is performed until sufficient damage is accumulated in the simulation cell.

Simulations are performed for values of the He concentrations varying from 1 to 25 appm/dpa (appm = atomic parts per million). Initial concentrations of both point defect entities are calculated from the NRT formula (Equation 4) using an incident energy of 100keV. Overall damage is introduced at the rate of 10-6 dpa/s and the simulation is carried out until the required amount of damage has accumulated in the system.

Figure 3 shows the concentration of bubbles as a function of He concentration (appm/dpa) after damage equal to 0.1 and 1 dpa is introduced. The simulation temperature is 0.3Tm. The bubble density increases with increasing He/dpa ratio. The bubble density can be described by a power law expression,

$$\mathcal{L}c\_B = K \left(\mathcal{c}\_{He}\right)^m \quad \text{\$\mathcal{L}\_B = K\$}\left(\mathcal{c}\_{He}\right)^m\$ \tag{7}$$

where is the bubble density, is the helium concentration expressed as appm He/dpa and *K* and *m* are constants determined by the kinetic Monte Carlo simulations. We find that the exponent *m* is approximately 0.5. In Figure 6, the solid lines are fit to the data assuming a square root dependence of the bubble density on helium concentration. Thus the bubble density increases as the square root of the He/dpa ratio. While experimental evidence of this variation is difficult to find, such dependence has been suggested by rate theory calculations (SINGH and TRINKAUS 1992) for the case of cold helium implantation annealing. The square root dependence is also found for equilibrium bubbles containing an ideal gas(Marksworth 1973).

The bubble density increases with damage (expressed as dpa) as damage is accumulated till 1 dpa. At higher dpa ratios more vacancies are produced that may serve as nucleation sites for embryonic bubbles by trapping helium. Thus the bubble density scales directly with increasing damage. That bubble density increases with accumulating dpa is observed in experiments (Sencer, Bond et al. 2001; Sencer, Garner et al. 2002; Zinkle, Hashimoto et al. 2003). These experiments suggest much smaller values of bubble density than those calculated in the present simulations. Embryonic bubbles are submicroscopic and difficult to estimate from experimental observations; while the kMC simulations have a large fraction of nanometer size bubbles; making comparison with experiment difficult.

Simulations were performed for a damage energy of 100keV. A production rate of randomly distributed cascades is assumed such that the damage is introduced at a rate of 10-6 dpa/s. The simulation cell size is 400a0x400a0x400a0, where a0 is the lattice parameter of iron. The number of Frenkel pairs introduced in the simulation cell are calculated by the Norgett-

where ED is the energy of incident particles while Ed is the threshold energy (40eV for iron (1994)). Initial concentration of helium is a parameter in the simulation and is varied from 1 to 25 appm/dpa. Damage measured in dpa is another parameter in the system and determines the length of the simulation. The incident energy is also a simulation control parameter that determines the number of defects introduced in a single cascade. The simulation is performed until sufficient damage is accumulated in the simulation cell.

Simulations are performed for values of the He concentrations varying from 1 to 25 appm/dpa (appm = atomic parts per million). Initial concentrations of both point defect entities are calculated from the NRT formula (Equation 4) using an incident energy of 100keV. Overall damage is introduced at the rate of 10-6 dpa/s and the simulation is carried

Figure 3 shows the concentration of bubbles as a function of He concentration (appm/dpa) after damage equal to 0.1 and 1 dpa is introduced. The simulation temperature is 0.3Tm. The bubble density increases with increasing He/dpa ratio. The bubble density can be described

where is the bubble density, is the helium concentration expressed as appm He/dpa and *K* and *m* are constants determined by the kinetic Monte Carlo simulations. We find that the exponent *m* is approximately 0.5. In Figure 6, the solid lines are fit to the data assuming a square root dependence of the bubble density on helium concentration. Thus the bubble density increases as the square root of the He/dpa ratio. While experimental evidence of this variation is difficult to find, such dependence has been suggested by rate theory calculations (SINGH and TRINKAUS 1992) for the case of cold helium implantation annealing. The square root dependence is also found for equilibrium bubbles containing an

The bubble density increases with damage (expressed as dpa) as damage is accumulated till 1 dpa. At higher dpa ratios more vacancies are produced that may serve as nucleation sites for embryonic bubbles by trapping helium. Thus the bubble density scales directly with increasing damage. That bubble density increases with accumulating dpa is observed in experiments (Sencer, Bond et al. 2001; Sencer, Garner et al. 2002; Zinkle, Hashimoto et al. 2003). These experiments suggest much smaller values of bubble density than those calculated in the present simulations. Embryonic bubbles are submicroscopic and difficult to estimate from experimental observations; while the kMC simulations have a large fraction of

nanometer size bubbles; making comparison with experiment difficult.

( )*<sup>m</sup> B He c Kc* <sup>=</sup> , ( )*<sup>m</sup>*

out until the required amount of damage has accumulated in the system.

0.8 2 *D d <sup>E</sup> displacements per cascade <sup>E</sup>* <sup>=</sup> (6)

*B He c Kc* = (7)

Robinson-Torrens relationship,

by a power law expression,

ideal gas(Marksworth 1973).

This example demonstrates the use of the kinetic Monte Carlo method to simulate defect diffusion in irradiated materials. The method depends on the ability to calculate rates of migration events of individual defects. The advantage of the method is that all atoms do not necessarily need to be simulated and large time scales are accessible to the simulation provided the underlying physics remains invariant.

Fig. 6. is a plot of the concentration of bubbles (cB) as a function of helium concentration (expressed in terms of appm/dpa) for two damage levels, 0.1 dpa (filled triangles) and 1 dpa (open squares). The line is a fit to the data assuming a square root dependence of the bubble density on helium concentration. The simulation temperature is 543K (0.3 Tm, the melting temperature for iron).

#### **4.4 Microstructural mechanics of irradiation hardening**

The previous examples have looked at the atomistic and mesoscale of radiation damage and defect formation. This information can be used by plasticity models and microstructural mechanics models of the effect of radiation on materials properties. Here an example is presented where the atomistic calculations are used to parameterize a viscoplasticity treatment of hardening in materials due to irradiation.

Hardening and embrittlement are controlled by interactions between dislocations and irradiation induced defect clusters. Radiation hardening and embrittlement that occurs in metals irradiated at low temperatures (below ~0.3 *Tm*, where *Tm* is the melting temperature) is a an important technical challenge for advanced nuclear energy systems(Zinkle and Matsukawa 2004). In this example, the Visco Plastic Self Consistent (VPSC) polycrystalline code (Lebensohn and Tome 1993) is employed in order to model the yield stress dependence

1998) are

clusters.

where τ α

is the initial CRSS,

τ0

(yield) stress is given by the average over all orientations.

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 277

metals might cause one to overestimate the reported barrier strength for the visible defect

It is possible to introduce a hardening law that is a function of the strain to describe the threshold resolved shear stress required to activate dislocations. In the present application, however, evolution is not simulated and only the initial threshold is required. We assume that the initial critical resolved shear stress (CRSS) in each grain is affected by irradiation according to the dispersed barrier hardening law and follows the Orowan expression,

the same meaning as in Equations 8 and 9. Observe that the Taylor factor is not included in Eq. 10, since the geometric crystal orientation effects are accounted for by the polycrystal model. The critical stress τ is assigned to the 12 (110)[111] and the 12 (110)[112] slip systems of the BCC structure. The initial texture of the rolled ferritic steel is represented using 1000 crystallographic orientations. Each orientation is treated as an ellipsoidal inclusion embedded in and interacting with the effective medium that represents the aggregate. An incremental strain is enforced along the rolling direction, while leaving the lateral strains unconstrained. The stress and the strain is different from grain to grain, and the macroscopic

Through Eq. 10 the model includes a dependence of the yield stress on the damage created due to radiation. Radiation damage is usually expressed as a statistical quantity describing the average number of displacements for each atom (dpa). The dpa influences the yield stress by determining the number density and the size of the defect clusters (obstacles) that impede the

It has commonly been assumed that the defect cluster density in irradiated metals increases linearly with increasing dose, up to the onset of cascade overlap which causes a saturation in the cluster density (Makin, Whapman et al. 1962; Koppenaal and Arsenault 1971; Trinkaus, Singh et al. 1996). However, in several pure FCC metals the defect accumulation as measured by electrical resistivity (Makin et al. 1962; Zinkle 1987) or transmission electron microscopy (Zinkle 1987; Muroga, Heinisch et al. 1992) often appears to exhibit an intermediate dose regime where the defect cluster density is proportional to the square root of dose. The defect accumulation behavior was found to be linear at very low doses (<0.0001 dpa, where the probability of uncorrelated point defect recombination is negligible), and proportional to the square root of dose at higher doses. According to simple kinetic models such as the unsaturable trap model (Thompson, Youngblood et al. 1973; Theis and Wollenberger 1980), the critical dose for transition from linear to square root behavior depends on specimen purity. In this model, the transition to square-root accumulation behavior can be delayed up to high doses if impurity trapping of migrating interstitial- type defects is dominant compared to interstitial– interstitial or interstitial–vacancy reactions.

The dependence on irradiation dose (expressed as dpa) of the defect cluster density (N) and the defect diameter (d) are taken from atomic level kinetic Monte Carlo (kMC) simulations and experimental observations (Deo et al. 2006; Deo, Baskes et al. 2007) The kMC model

path of the dislocations and increase the critical stress required to move the dislocation.

0 τ τ αμ

= 0.4 or higher. It is possible that hardening from atomic scale voids in the BCC

= + *b Nd* , (10)

is the unirradiated initial CRSS and the other parameters have

in ferritic steels on the irradiation dose. The dispersed barrier hardening model is implemented in the VPSC code by introducing a hardening law, function of the strain, to describe the threshold resolved shear stress required to activate dislocations. The size and number density of the defect clusters varies with the irradiation dose in the model. Such modeling efforts can both reproduce experimental data and also guide future experiments of irradiation hardening.

In order to describe the nature of the yield stress dependence on the irradiation dose, we implemented a new microstructural model at the grain level in the VPSC code. The model assumes that hardening is affected by the presence of the defects and defect clusters produced by irradiation. These defects interact with the pre-existing dislocations in the microstructure leading to an increase in the critical stress necessary to move the dislocations. This leads to an increase in the overall yield stress of the material.

Defects are treated as barriers to the motion of dislocations. Two approximate dislocation barrier models have historically been used to describe radiation hardening in metals (Zinkle and Matsukawa 2004) and are reviewed in (Koppenaal and Arsenault 1971; Kocks 1977). The dispersed barrier model (Seeger, Diehl et al. 1957) is based on straightforward geometrical considerations for obstacles intersecting the dislocation glide plane and it is most appropriate for strong obstacles. An alternative hardening relationship was developed by Friedel–Kroupa–Hirsch (FKH) for weak obstacles (Friedel 1955; Kroupa and Hirsch 1964), where the effective interparticle spacing is increased compared to the planar geometric spacing due to less extensive dislocation bowing prior to obstacle breakaway. Using the simple approximation for dislocation line tension, the functional dependence of polycrystalline yield strength increase on defect cluster size and density for these two limiting cases is given by the following equations:

$$
\Delta \sigma = M a \sharp b \circ \text{N} d \,, \tag{8}
$$

$$
\Delta \sigma = \frac{1}{8} M \mu b N^{\frac{2}{3}} d\,, \tag{9}
$$

Equation 8 corresponds to the dispersed barrier hardening model and Equation 9 to the FKH model. In the two equations, Δσ is the change in the yield stress, M is the Taylor factor (3.06 for non-textured BCC and FCC metals), α is the defect cluster barrier strength, μ is the shear modulus, *b* is the Burgers vector of the primary glide dislocations, and *N* and *d* are the defect cluster density and diameter.

Most radiation hardening studies have used the dispersed barrier model (Equation 8) for data interpretation, and in this work we find that it provides a better representation of our experimental results. However, the FKH model (Equation 9) may be more appropriate for many radiation-induced small defect clusters which are weak obstacles to dislocation motion. According to some early analyses (Kocks 1977), the FKH model is adequate for barrier strengths up to 1/4 of the Orowan (impenetrable obstacle) limit, i.e., α < 0:25, and the dispersed barrier model is more appropriate for barrier strengths of 0.25<α <1. Typical experimental values of the defect cluster barrier strength for copper and austenitic stainless steel neutron-irradiated and tested near room temperature are α = 0.15–0.2 (Zinkle 1987). The reported barrier strengths for the visible defect clusters in BCC metals (Rice and Zinkle

in ferritic steels on the irradiation dose. The dispersed barrier hardening model is implemented in the VPSC code by introducing a hardening law, function of the strain, to describe the threshold resolved shear stress required to activate dislocations. The size and number density of the defect clusters varies with the irradiation dose in the model. Such modeling efforts can both reproduce experimental data and also guide future experiments

In order to describe the nature of the yield stress dependence on the irradiation dose, we implemented a new microstructural model at the grain level in the VPSC code. The model assumes that hardening is affected by the presence of the defects and defect clusters produced by irradiation. These defects interact with the pre-existing dislocations in the microstructure leading to an increase in the critical stress necessary to move the dislocations.

Defects are treated as barriers to the motion of dislocations. Two approximate dislocation barrier models have historically been used to describe radiation hardening in metals (Zinkle and Matsukawa 2004) and are reviewed in (Koppenaal and Arsenault 1971; Kocks 1977). The dispersed barrier model (Seeger, Diehl et al. 1957) is based on straightforward geometrical considerations for obstacles intersecting the dislocation glide plane and it is most appropriate for strong obstacles. An alternative hardening relationship was developed by Friedel–Kroupa–Hirsch (FKH) for weak obstacles (Friedel 1955; Kroupa and Hirsch 1964), where the effective interparticle spacing is increased compared to the planar geometric spacing due to less extensive dislocation bowing prior to obstacle breakaway. Using the simple approximation for dislocation line tension, the functional dependence of polycrystalline yield strength increase on defect cluster size and density for these two

> Δ = σ

Δσ

barrier strengths up to 1/4 of the Orowan (impenetrable obstacle) limit, i.e.,

the dispersed barrier model is more appropriate for barrier strengths of 0.25<

steel neutron-irradiated and tested near room temperature are

 αμ

> μ

Equation 8 corresponds to the dispersed barrier hardening model and Equation 9 to the

α

shear modulus, *b* is the Burgers vector of the primary glide dislocations, and *N* and *d* are the

Most radiation hardening studies have used the dispersed barrier model (Equation 8) for data interpretation, and in this work we find that it provides a better representation of our experimental results. However, the FKH model (Equation 9) may be more appropriate for many radiation-induced small defect clusters which are weak obstacles to dislocation motion. According to some early analyses (Kocks 1977), the FKH model is adequate for

experimental values of the defect cluster barrier strength for copper and austenitic stainless

The reported barrier strengths for the visible defect clusters in BCC metals (Rice and Zinkle

8 Δ = σ

2 <sup>3</sup> 1

*M b Nd* , (8)

*M bN d* , (9)

μis the

< 0:25, and

<1. Typical

α

α

= 0.15–0.2 (Zinkle 1987).

is the change in the yield stress, M is the Taylor factor

is the defect cluster barrier strength,

α

This leads to an increase in the overall yield stress of the material.

limiting cases is given by the following equations:

FKH model. In the two equations,

defect cluster density and diameter.

(3.06 for non-textured BCC and FCC metals),

of irradiation hardening.

1998) are α = 0.4 or higher. It is possible that hardening from atomic scale voids in the BCC metals might cause one to overestimate the reported barrier strength for the visible defect clusters.

It is possible to introduce a hardening law that is a function of the strain to describe the threshold resolved shear stress required to activate dislocations. In the present application, however, evolution is not simulated and only the initial threshold is required. We assume that the initial critical resolved shear stress (CRSS) in each grain is affected by irradiation according to the dispersed barrier hardening law and follows the Orowan expression,

$$
\pi = \pi\_0 + \alpha \mu b \sqrt{\text{Nd}} \tag{10}
$$

where τ is the initial CRSS, τ0 is the unirradiated initial CRSS and the other parameters have the same meaning as in Equations 8 and 9. Observe that the Taylor factor is not included in Eq. 10, since the geometric crystal orientation effects are accounted for by the polycrystal model. The critical stress τ is assigned to the 12 (110)[111] and the 12 (110)[112] slip systems of the BCC structure. The initial texture of the rolled ferritic steel is represented using 1000 crystallographic orientations. Each orientation is treated as an ellipsoidal inclusion embedded in and interacting with the effective medium that represents the aggregate. An incremental strain is enforced along the rolling direction, while leaving the lateral strains unconstrained. The stress and the strain is different from grain to grain, and the macroscopic (yield) stress is given by the average over all orientations.

Through Eq. 10 the model includes a dependence of the yield stress on the damage created due to radiation. Radiation damage is usually expressed as a statistical quantity describing the average number of displacements for each atom (dpa). The dpa influences the yield stress by determining the number density and the size of the defect clusters (obstacles) that impede the path of the dislocations and increase the critical stress required to move the dislocation.

It has commonly been assumed that the defect cluster density in irradiated metals increases linearly with increasing dose, up to the onset of cascade overlap which causes a saturation in the cluster density (Makin, Whapman et al. 1962; Koppenaal and Arsenault 1971; Trinkaus, Singh et al. 1996). However, in several pure FCC metals the defect accumulation as measured by electrical resistivity (Makin et al. 1962; Zinkle 1987) or transmission electron microscopy (Zinkle 1987; Muroga, Heinisch et al. 1992) often appears to exhibit an intermediate dose regime where the defect cluster density is proportional to the square root of dose. The defect accumulation behavior was found to be linear at very low doses (<0.0001 dpa, where the probability of uncorrelated point defect recombination is negligible), and proportional to the square root of dose at higher doses. According to simple kinetic models such as the unsaturable trap model (Thompson, Youngblood et al. 1973; Theis and Wollenberger 1980), the critical dose for transition from linear to square root behavior depends on specimen purity. In this model, the transition to square-root accumulation behavior can be delayed up to high doses if impurity trapping of migrating interstitial- type defects is dominant compared to interstitial– interstitial or interstitial–vacancy reactions.

The dependence on irradiation dose (expressed as dpa) of the defect cluster density (N) and the defect diameter (d) are taken from atomic level kinetic Monte Carlo (kMC) simulations and experimental observations (Deo et al. 2006; Deo, Baskes et al. 2007) The kMC model

Multiscale Materials Modeling of Structural Materials for Next Generation Nuclear Reactors 279

studies that can inform future atomic-scale studies of dislocation – obstacle interactions. The VPSC model could also incorporate experimentally observed defect cluster distributions, number densities to assess the effect of multiple defect types and distributions. A detailed multiscale study, wherein the dislocation-obstacle strength and the number density and size of defects are correlated to the increasing strain and each other, would then further explain the

The VPSC calculations provide a means to link atomistic first principles calculations to macroscopic observables. The formulation of the irradiation hardening law allows for the introduction of parameters such as the defect size and number density that can be calculated from evolution models and simulations such as the kinetic Monte Carlo method. The interaction of the dislocation with defect clusters can be investigated by using atomistic molecular dynamics calculations. In this document we have provided a framework for performing physically based modeling and simulations of hardening behavior observed during irradiation. Such modeling efforts can both reproduce experimental data and also guide future experiments of irradiation hardening. Performing modeling and simulation studies before initiating an expensive neutron or proton beam experiment would prove

In this chapter, an overview of multiscale materials modeling tools used to simulate structural materials in irradiation conditions is presented. Next generation nuclear reactors will require a new generation of materials that can survive and function in extreme environments. Advanced modeling and simulation tools can study these materials at various length and time scale. Such varied methods are needed as radiation damage affects materials in excess of 10 orders of magnitude in length scale from the sub-atomic nuclear to structural component level, and span 22 orders of magnitude in time from the subpicosecond of nuclear collisions to the decade-long component service lifetimes. The inherently wide range of time scales and the "rare-event" nature of the controlling mechanisms make modeling radiation effects in materials extremely challenging and experimental characterization is often unattainable. Thus, modeling and simulation of such materials holds great promise if coupled with suitably designed experiments in order to

The author would like to thank Erin Hayward, Carlos Tome, Richardo Lebensohn, Stuart Maloy, Maria Okuniewski, Michael Baskes, Srinivasan Srivilliputhur and Michael James for useful discussions and acknowledge funding sources from the Department of Energy (DOE NEUP award DE-AC07-05ID14517 09-269) and a Nuclear Regulatory Commission Faculty

Arsenlis, A., B. D. Wirth, et al. (2004). "Dislocation density-based constitutive model for the

mechanical behaviour of irradiated Cu." Philosophical Magazine 84(34): 3617-3635.

effect of irradiation on mechanical properties of ferritic steels.

develop and sustain materials for advanced nuclear energy.

invaluable and cost-effective.

**6. Acknowledgements** 

**7. References** 

Development Grant NRC-38-08-938

**5. Conclusions** 

takes atomic level information of the migration energies and jump attempt frequencies of irradiation induced defects (interstitials, vacancies) and transmutation products (e.g., helium under high energy proton irradiation), and evolves the microstructure according to the rates of migration of these defects. The defects are allowed to cluster, and new irradiation damage is introduced during the simulation according to the irradiation dose rate. Our kMC simulations predict that the number density varies as the square root of the displacements per atom for the case of bcc iron irradiated up to 1 dpa by high energy proton irradiation.

The size dependence on irradiation dose is more complicated as the kMC simulations provide an entire distribution of defect cluster sizes. A single value of d as a function of dose is still a simplification of the kMC results. *.* The defect size usually increases with increasing dose (dpa) and can be fit by a power law; however the exponent of the power law expression can vary from 0 to 0.5 depending on initial simulation conditions (dose rate, temperature) and the defect cluster size considered. At low dpa, the exponent of the power law dependence is small for all defect sizes and increases at higher dpa.

The density of defects N is assumed to vary as the square root of the dpa while two cases of size dependence are considered, one in which the size is invariant with the dose (dpa) while the other in which the defect size varies as the square root of the dose. Additional systematic work is needed to confirm the presence and to understand the physical mechanisms responsible for this square root fluence-dependent defect cluster accumulation regime.

The link between the atomic level simulations and the VPSC calculations was established using the dispersed barrier hardening model. In this model, the vacancy /interstitial clusters produced in radiation cascades are assumed to act as barriers to the gliding dislocation in the slip plane and are therefore taken to be the main source of radiation hardening. A different model of radiation hardening postulates the formation of defect clouds along the length of the grown-in dislocation( see [4,5] for review). These clouds prevent the dislocation from acting as Frank Read dislocation sources and emitting more dislocations. Singh, Golubov et al. (1997) proposed the cascade induced source hardening model which accounted for interstitial cluster formation during radiation cascade formation. Such cluster formation has been observed in molecular dynamics simulations. In the CISH model, glissile loops produced directly in cascades are assumed to decorate grown-in dislocations so they cannot act as dislocation sources. The yield stress is related to the breakaway stress which is necessary to pull the dislocation away from the clusters/loops decorating it. Various aspects of the model (main assumptions and predictions) have been investigated by these researchers using analytical calculations, 3-D dislocation dynamics and molecular dynamics simulations It is possible to investigate such recent radiation hardening mechanisms by including them to develop the links between the atomic level understanding of defect sizes and concentrations and the VPSC model of polycrystalline hardening. Such mechanisms may also be investigated by atomic level simulations of single dislocation motion in the presence of defect impurities.

In a manner similar to the approach of Arsenelis and co-workers (Arsenlis, Wirth et al. 2004), the VPSC model can be used to combine microstructural input from both experimental observations and model predictions to evaluate the contributions from multiple defect cluster types. Although not all of the relevant parameters are currently known, such parameterstudies that can inform future atomic-scale studies of dislocation – obstacle interactions. The VPSC model could also incorporate experimentally observed defect cluster distributions, number densities to assess the effect of multiple defect types and distributions. A detailed multiscale study, wherein the dislocation-obstacle strength and the number density and size of defects are correlated to the increasing strain and each other, would then further explain the effect of irradiation on mechanical properties of ferritic steels.

The VPSC calculations provide a means to link atomistic first principles calculations to macroscopic observables. The formulation of the irradiation hardening law allows for the introduction of parameters such as the defect size and number density that can be calculated from evolution models and simulations such as the kinetic Monte Carlo method. The interaction of the dislocation with defect clusters can be investigated by using atomistic molecular dynamics calculations. In this document we have provided a framework for performing physically based modeling and simulations of hardening behavior observed during irradiation. Such modeling efforts can both reproduce experimental data and also guide future experiments of irradiation hardening. Performing modeling and simulation studies before initiating an expensive neutron or proton beam experiment would prove invaluable and cost-effective.
