Category Archives: Comprehensive nuclear materials

Radiation Damage Theory

1.13.1 Introduction

The study of radiation effects on the structure and properties of materials started more than a century ago,1 but gained momentum from the development of fission reactors in the 1940s. In 1946, Wigner2 pointed out the possibility of a deleterious effect on material properties at high neutron fluxes, which was then confirmed experimentally.3 A decade later, Konobeevsky et al4 discovered irradiation creep in fissile metallic uranium, which was then observed in stainless steel.5 The discovery of void swelling in neutron-irradiated stainless steels in 1966 by Cawthorne and Fulton6 demonstrated that radiation effects severely restrict the lifetime ofreactor materials and that they had to be systematically studied.

The 1950s and early 1960s were very productive in studying crystalline defects. It was recognized that atoms in solids migrate via vacancies under thermal — equilibrium conditions and via vacancies and self­interstitial atoms (SIAs) under irradiation; also that the bombardment with energetic particles generates high concentrations of defects compared to equilib­rium values, giving rise to radiation-enhanced diffu­sion. Numerous studies revealed the properties of point defects (PDs) in various crystals. In particular, extensive studies of annealing of irradiated samples resulted in categorizing the so-called ‘recovery stages’ (e. g., Seeger7), which comprised a solid basis for understanding microstructure evolution under irradiation.

Already by this time, which was well before the discovery of void swelling in 1966, the process of interaction of various energetic particles with solid

targets had been understood rather well (e. g., Kinchin and Pease8 for a review). However, the primary dam­age produced was wrongly believed to consist of Frenkel pairs (FPs) only. In addition, it was com­monly believed that this damage would not have serious long-term consequences in irradiated materi­als. The reasoning was correct to a certain extent; as they are mobile at temperatures of practical interest, the irradiation-produced vacancies and SIAs should move and recombine, thus restoring the original crys­tal structure. Experiments largely confirmed this sce­nario, most defects did recombine, while only about 1% or an even smaller fraction survived and formed vacancy and SIA-type loops and other defects. How­ever small, this fraction had a dramatic impact on the microstructure of materials, as demonstrated by Cawthorne and Fulton.6 This discovery initiated extensive experimental and theoretical studies of radiation effects in reactor materials which are still in progress today.

After the discovery of swelling in stainless steels, it was found to be a general phenomenon in both pure metals and alloys. It was also found that the damage accumulation takes place under irradiation with any particle, provided that the recoil energy is higher than some displacement threshold value, £d, (^30-40 eV in metallic crystals). In addition, the microstructure of different materials after irradiation was found to be quite similar, consisting of voids and dislocation loops. Most surprisingly, it was found that the microstructure developed under irradiation with ~1 MeV electrons, which produces FPs only, is similar to that formed under irradiation with fast neutrons or heavy-ions, which produce more compli­cated primary damage (see Singh et a/.1). All this created an illusion that three-dimensional migrating (3D) PDs are the main mobile defects under any type of irradiation, an assumption that is the foundation of the initial kinetic models based on reaction rate the­ory (RT). Such models are based on a mean-field approximation (MFA) of reaction kinetics with the production of only 3D migrating FPs. For conve­nience, we will refer to these models as FP produc­tion 3D diffusion model (FP3DM) and henceforth this abbreviation will be used. This model was devel­oped in an attempt to explain the variety of phenom­ena observed: radiation-induced hardening, creep, swelling, radiation-induced segregation (RIS), and sec­ond phase precipitation. A good introduction to this theory can be found, for example, in the paper by Sizmann,9 while a comprehensive overview was pro­duced by Mansur,10 when its development was already completed. The theory is rather simple, but its general methodology can be useful in the further development of radiation damage theory (RDT). It is valid for ^1MeV electron irradiation and is also a good introduction to the modern RDT, see Section

1.13.5.

Soon after the discovery of void swelling, a number of important observations were made, for example, the void super-lattice formation11-14 and the microm­eter-scale regions of the enhanced swelling near grain boundaries (GBs).15 These demonstrated that under neutron or heavy-ion irradiation, the material micro­structure evolves differently from that predicted by the FP3DM. First, the spatial arrangement of irradia­tion defects voids, dislocations, second phase particles, etc. is not random. Second, the existence of the micrometer-scale heterogeneities in the microstruc­ture does not correlate with the length scales accounted for in the FP3DM, which are an order of magnitude smaller. Already, Cawthorne and Fulton6 in their first publication on the void swelling had reported a nonrandomness of spatial arrangement of voids that were associated with second phase precipi­tate particles. All this indicated that the mechanisms operating under cascade damage conditions (fast neu­tron and heavy-ion irradiations) are different from those assumed in the FP3DM. This evidence was ignored until the beginning of the 1990s, when the production bias model (PBM) was put forward by Woo and Singh.16,17 The initial model has been changed and developed significantly since then18-2 and explained successfully such phenomena as high swelling rates at low dislocation density (Section

1.13.6.2.2) , grain boundary and grain-size effects in void swelling, and void lattice formation (Section

1.13.6.2.3) . An essential advantage of the PBM over the FP3DM is the two features of the cascade dam­age: (1) the production of PD clusters, in addition to single PDs, directly in displacement cascades, and (2) the 1D diffusion of the SIA clusters, in addition to the 3D diffusion of PDs (Section 1.13.3). The PBM is, thus, a generalization of the FP3DM (and the idea of intracascade defect clustering intro­duced in the model by Bullough et a/. (BEK29)). A short overview of the PBM was published about 10 years ago.1 Here, it will be described somewhat differently, as a result of better understanding of what is crucial and what is not, see Section 1.13.6.

From a critical point of view, it should be noted that successful applications of the PBM have been limited to low irradiation doses (< 1 dpa) and pure metals (e. g., copper). There are two problems that prevent it from being used at higher doses. First, the PBM in its present form1 predicts a saturation of void size (see, e. g., Trinkaus et at19 and Barashev and Golubov30 and Section 1.13.6.3.1). This originates from the mixture of 1D and 3D diffusion-reaction kinetics under cascade damage conditions, hence from the assumption lying at the heart of the model. In contrast, experiments demonstrate unlimited void growth at high doses in the majority of materials and conditions (see, e. g., Singh et at.,31 Garner,32 Garner et at.,33 and Matsui et at.34). An attempt to resolve this contradiction was undertaken23,25,27 by including thermally activated rotations of the SIA-cluster Burgers vector; but it has been shown25 that this does not solve the problem. Thus, the PBM in its present form fails to account for the important and common observation: the indefinite void growth under cascade irradiation. The second problem of the PBM is that it fails to explain the swelling satura­tion observed in void lattices (see, e. g., Kulchinski et at.13). In contrast, it predicts even higher swelling rates in void lattices than in random void arrange — ments.25 This is because of free channels between voids along close-packed directions, which are formed during void ordering and provide escape routes for 1D migrating SIA clusters to dislocations and GBs, thus allowing 3D migrating vacancies to be stored in voids.

Resolving these two problems would make PBM self-consistent and complete its development. A solu­tion to the first problem has recently been proposed by Barashev and Golubov35,36 (see Section 1.13.7). It has been suggested that one of the basic assumptions of all current models, including the PBM, that a random arrangement of immobile defects exists in the material, is correct at low and incorrect at high doses. The analysis includes discussion of the role of RIS and provides a solution to the problem, making the PBM capable of describing swelling in both pure metals and alloys at high irradiation doses. The solu­tion for the second problem of the PBM mentioned above is the main focus of a forthcoming publication by Golubov et at.37

Because of limitations of space, we only give a short guide to the main concepts of both old and more recent models and the framework within which radiation effects, such as void swelling, and hardening and creep, can be rationalized. For the same reason, the impact of radiation on reactor fuel materials is not considered here, despite a large body of relevant experimental data and theoretical results collected in this area.

Modeling Challenges to Predict Irradiation Effects on Materials

The effect of irradiation on materials is a classic example of an inherently multiscale phenomenon, as schematically illustrated in Figure 1. Pertinent processes span over more than 10 orders of magni­tude in length scale from the subatomic nuclear to

Nano/microstructure and local chemistry changes; nucleation and growth of extended defects and precipitates I

image889

ф

 

Подпись: He and H generation

Подпись: Atomic-nm nm-pm pm-mm mm-m
Подпись: Lengthscale
Подпись: Figure 1 Illustration of the length and time scales (and inherent feedback) involved in the multiscale processes responsible for microstructural changes in irradiated materials. Reproduced from Wirth, B. D.; Odette, G. R.; Marian, J.; Ventelon, L.; Young, J. A.; Zepeda-Ruiz, L. A. J. Nucl. Mater. 2004, 329-333, 103.

Defect

recombination, clustering, and migration Primary defect production and short-term annealing

Underlying
microstructure
(preexisting and evolving)
impacts defect and solute
fate

structural component level, and span 22 orders of magnitude in time from the subpicosecond of nuclear collisions to the decade-long component service life­times.1,2 Many different variables control the mix of nano/microstructural features formed and the corresponding degradation of physical and mechanical properties of nuclear fuels, cladding, and structural materials. The most important variables include the initial material composition and microstructure, the thermomechanical loads, and the irradiation history. While the initial material state and thermomechanical loading are of concern in all materials performance — limited engineering applications, the added complexity introduced by the effects of radiation is clearly the distinguishing and overarching concern for materials in advanced nuclear energy systems.

At the smallest scales, radiation damage is continu­ally initiated by the formation of energetic primary knock-on atoms (PKAs) primarily through elastic col­lisions with high-energy neutrons. Concurrently, high concentrations of fission products (in fuels) and trans­mutants (in cladding and structural materials) are gen­erated and can cause pronounced effects in the overall chemistry of the material, especially at high burnup. The PKAs, 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 col­lision displacements, generating a cascade of vacancy and self-interstitial defects. High-energy displacement cascades evolve over very short times, 100 ps or less, and small volumes, with characteristic length scales of 50 nm or less, and are directly amenable to MD simula­tions if accurate potentials are available.

The physics of primary damage production in high-energy displacement cascades has been exten­sively studied with MD simulations.3-8 (see Chapter

1.11, Primary Radiation Damage Formation) The key conclusions from the MD studies of cascade evo­lution have been that (1) intracascade recombination of vacancies and self-interstitial atoms (SIAs) results in ^30% of the defect production expected from displacement theory, (2) many-body collision effects produce a spatial correlation (separation) of the va­cancy and SIA defects, (3) substantial clustering of the SIAs and to a lesser extent the vacancies occur within the cascade volume, and (4) high-energy displace­ment cascades tend to break up into lobes or subcas­cades, which may also enhance recombination.4-7

Nevertheless, it is the subsequent diffusional trans­port and evolution of the defects produced during displacement cascades, in addition to solutes and transmutant impurities, that ultimately dictate radia­tion effects in materials and changes in material micro — structure.1, Spatial correlations associated with the displacement cascades continue to play an important role in much larger scales as do processes including defect recombination, clustering, migration, and gas and solute diffusion and trapping. Evolution of the underlying materials structure is thus governed by the time and temperature kinetics of diffusive and reactive processes, albeit strongly influenced by spatial correlations associated with the microstructure and the continuous production of new radiation damage.

The inherently wide range of time scales and the ‘rare-event’ nature of many of the controlling mech­anisms make modeling radiation effects in materials extremely challenging and experimental characteri­zation often unattainable. Indeed, accurate models of microstructure (point defects, dislocations, and grain boundaries) evolution during service are still lacking. To understand the irradiation effects and microstructure evolution to the extent required for a high fidelity nuclear materials performance model will require a combination of experimental, theoreti­cal, and computational tools.

Furthermore, the kinetic processes controlling defect cluster and microstructure evolution, as well as the materials degradation and failure modes may not entirely be known. Thus, a substantial challenge is to discover the controlling processes so that they can be included within the models to avoid the detri­mental consequences of in-service surprises. High performance computing can enable such discovery of class simulations, but care must also be taken to assess the accuracy of the models in capturing critical physical phenomena. The remainder of this chapter will thus focus on a description of KMC modeling, along with a few select examples of the application of KMC models to predict irradiation effects on materials and to identify opportunities for additional research to achieve the goal of accelerating the devel­opment of advanced computational approaches to simulate nucleation, growth, and coarsening of mi­crostructure in complex engineering materials.

Treatment of Solutions

Whether it is the nonstoichiometry of fluorite — structure UO2±x or variable composition ortho­rhombic or tetragonal U-Zr alloy fuel, the accurate thermochemical description of these phases has

 

image422

image1007

Подпись: Figure 1Diagram illustrating the computer coupling of phase diagrams and thermochemistry approach after Zinkevich.2

been through the use of solution models. Solid and liquid solution modeling from simple highly dilute systems to more complex interstitial and substitu­tional solutions with multiple lattices has been a rich field for some time, yet it is far from fully developed. To accurately describe the energetics of solutions will eventually require bridging atomistic models with the mesoscale treatments currently being used. Recent approaches have begun to deal with defect structures in phases, although only in a very constrained manner which limits clustering and other phenomena. Yet, when they are coupled with accurate data that allow fitting of the model para­meters, the resulting representations have been highly predictive of phase relations and chemistry.

A number of texts provide useful descriptions of solution models from the simple to the complex.12-15 While there are several relatively accurate but rather intricate approaches, such as the cluster variation method, the discussion in this chapter is confined to
simpler models that are easily used in thermochemi­cal equilibrium computational software and thus applicable to large, multicomponent systems of inter­est in nuclear fuels.

The simplest model is the ideal solution where the constituents are assumed to mix randomly with no structural constraints and no interactions (bonding or short — or long-range order). The standard Gibbs free energy and ideal mixing entropy contributions are

G° = X »<G?’

Gid = RT^2"j ln n, [8]

where G is the weighted sum of the standard Gibbs free energy for the constituents in the j phase solution, Gid is the contribution from the entropy increase due to randomly mixing the constituents, which is the configurational entropy, and R is the ideal gas law constant. For an ideal solution, the sum of the two represents the Gibbs free energy of the system.

In cases where there are significant interactions (bonding or repulsive interaction energies) among mixing constituents, an energetic term or terms need to be added to the solution free energy. The inclusion of a simple compositionally weighted excess energy term, Gex, accounts for the additional solution energy for what is historically termed a regular solution

Gex = XX XXEj [9

i=1 j>1

where X is the mole fraction and Ej the interaction energy between the components i and j. The system free energy is thus

G = G° + Gid + Gex [10]

Challenges of the RIS Continuous Models

Ab initio calculations have become a very powerful tool for RIS simulations. They have been shown to be able to provide not only the variation of the atomic jump frequencies with local concentration,125 but also new diffusion mechanisms.7 From a precise knowledge of the atomic jump frequencies and the recent develop­ment of diffusion models,98 a quantitative description of the flux coupling is expected to be feasible even in concentrated alloys. A unified description of flux coupling in dilute and concentrated alloys would allow the simultaneous prediction of two different mechanisms leading to RIS: solute drag by vacancies, and an IK effect involving the major elements.

An RIS segregation profile spreads over nan­ometers. Cell sizes of RIS continuous models are then too small for the theory of TIP to be valid. A mesoscale approach that includes interface energy between cells, such as the Cahn-Hilliard method, would be more appropriate. A derivation of quantita­tive phase field equations with fluctuations has recently been published.62 The resulting kinetic equations are dependent on the local concentration and also cell-size dependent. However, the diffusion mechanism involved direct exchanges between atoms. The same method needs to be developed for a system with point defect diffusion mechanisms.

Although it has been suggested since the first publications30 that the vacancy flux opposing the set up of RIS could slow down the void swelling, the change of microstructure and its coupling with RIS have almost never been modeled. Only very recently, phase field methods have tackled the kinetics of a concentration field and its interaction with a cavity population (see Chapter 1.15, Phase

Field Methods). The development of a simulation tool able to predict the mutual interaction between the point defect microstructure and the flux coupling is quite a challenge.

Motivation for Using Ion Beams to Study Radiation Damage

In the 1960s and 1970s, heavy ion irradiation was developed for the study of radiation damage pro­cesses in materials. As ion irradiation can be conducted at a well-defined energy, dose rate, and temperature, it results in very well-controlled experi­ments that are difficult to match in reactors. As such, interest grew in the use of ion irradiation for the purpose of simulating neutron damage in support of the breeder reactor program.1-3 Ion irradiation and simultaneous He injection were also used to simulate the effects of 14MeV neutron damage in conjunction with the fusion reactor engineering program. The application of ion irradiation (defined here as irradi­ation by any charged particle, including electrons) to the study of neutron irradiation damage caught the interest of the light water reactor community to address issues such as swelling, creep, and irradiation assisted stress corrosion cracking of core structural materials.4-6 Ion irradiation was also being used to understand the irradiated microstructure of reactor pressure vessel steels, Zircaloy fuel cladding, and materials for advanced reactor concepts.

There is significant incentive to use ion irradiation to study neutron damage as this technique has the potential for yielding answers on basic processes in addition to the potential for enormous savings in time and money. Neutron irradiation experiments are not amenable to studies involving a wide range of conditions, which is precisely what is required for investigations of the basic damage processes. Simula­tion by ions allows easy variation of the irradiation parameters such as dose, dose rate, and temperature over a wide range of values.

One of the prime attractions of ion irradiation is the rapid accumulation of end of life doses in short periods of time. Typical neutron irradiation experi­ments in thermal test reactors may accumulate dam­age at a rate of 3-5 dpayear-1. In fast reactors, the rates can be higher, on the order of 20 dpa year. For low dose components such as structural components in boiling water reactor (BWR) cores that typically have an end-of-life damage of 10 dpa, these rates are acceptable. However, even the higher dose rate of a fast reactor would require 4-5 years to reach the peak dose of ^80 dpa in the core baffle in a pressurized water reactor (PWR). For advanced, fast reactor con­cepts in which core components are expected to receive 200 dpa, the time for irradiation in a test reactor becomes impractical.

In addition to the time spent ‘in-core,’ there is an investment in capsule design and preparation as well as disassembly and allowing for radioactive decay, add­ing additional years to an irradiation program. Analysis of microchemical and microstructural changes by atom probe tomography (APT), Auger electron spec­troscopy (AES) or microstructural changes by energy dispersive spectroscopy via scanning transmission electron microscopy (STEM-EDS) and mechanical property or stress corrosion cracking (SCC) evalua­tion can take several additional years because of the precautions, special facilities, and instrumentation required for handling radioactive samples. The result is that a single cycle from irradiation through micro­analysis and mechanical property/SCC testing may require over a decade. Such a long cycle length does not permit for iteration of irradiation or material conditions that is critical in any experimental research program. The long cycle time required for design and irradiation also reduces flexibility in altering irradia­tion programs as new data become available. The requirement of special facilities, special sample handling, and long irradiation time make the cost for neutron irradiation experiments very high.

In contrast to neutron irradiation, ion (heavy, light, or electrons) irradiation enjoys considerable advantages in both cycle length and cost. Ion irradia­tions of any type rarely require more than several tens of hours to reach damage levels in the 1-100 dpa range. Ion irradiation produces little or no residual radioactivity, allowing handling of samples without

Подпись: VNRT(T )Подпись: [1]Подпись:Подпись: dpathe need for special precautions. These features translate into significantly reduced cycle length and cost. The challenge then is to verify the equivalency between neutron and ion irradiation in terms of the changes to the microstructure and properties of the material. The key question that needs to be answered is how do results from neutron and charged particle irradiation experiments compare? How, for example, is one to compare the results of a component irradiated in-core at 288 °C to a fluence of 1 x 1021n cm~ (E > 1 MeV) over a period of one year, with an ion irradiation experiment using 3 MeV protons at 400 °C to 1 dpa (displacements per atom) at a dose rate of 10~5dpas_1 (~1day), or 5MeV Ni2+ at 500°C to 10dpa at a dose rate of 5 x 10~3 dpas-1 (~1 h)?

The first question to resolve is the measure of radia­tion effect. In the Irradiation assisted stress corrosion cracking (IASCC) problem in LWRs, concern has cen­tered on two effects of irradiation: radiation-induced segregation of major alloying elements or impurities to grain boundaries, which may cause embrittlement or enhance the intergranular stress corrosion cracking (IGSCC) process, and hardening of the matrix that results in localized deformation and embrittlement. The appropriate measure of the radiation effect in the former case would then be the alloy concentration at the grain boundary or the amount of impurity segregated to the grain boundary. This quantity is measurable by analytical techniques such as AES, APT, or STEM-EDS. For the latter case, the measure of the radiation effect would be the nature, size, density, and distribution of dislocation loops, black dots, and the total dislocation network, and how they impact the deformation of the alloy. Hence, specific and mea­surable effects of irradiation can be determined for both neutron and ion irradiation experiments.

The next concern is determining how ion irradia­tion translates into the environment describing neu­tron irradiation. That is, what are the irradiation conditions required for ion irradiation to yield the same measure of radiation effect as that for neutron irradiation? This is the key question, for in a postirra­diation test program, it is only the final state of the material that determines equivalence, not the path taken. Therefore, if ion irradiation experiments could be devised that yielded the same measures of irradiation effects as observed in neutron irradiation experiments, the data obtained in postirradiation experiments will be equivalent. In such a case, ion irradiation experiments can provide a direct substitute for neutron irradiation. While neutron irradiation will always be required to qualify materials for reactor
application, ion irradiation provides a low-cost and rapid means of elucidating mechanisms and screening materials for the most important variables.

A final challenge is the volume of material that can be irradiated with each type of radiation. Neutrons have mean free paths on the order of centimeters in structural materials. One MeV electrons penetrate about 500 pm, 1 MeV protons penetrate about 10 pm, and 1 MeV Ni ions have a range of less than 1 pm. Thus, the volume of material that can be irradiated with ions from standard laboratory-sized sources (TEMs, accelerators), is limited.

Ab Initio for Irradiation

Irradiation damage, especially cascade modeling, is usually preferentially dealt by larger scale methods such as molecular dynamics with empirical potentials rather than ab initio calculations. However, recently ab initio studies that directly tackle irradiation pro­cesses have appeared.

1.08.3.3.1 Threshold displacement energies

First, the increase in computer power has allowed the calculations of threshold displacement energies by ab initio molecular dynamics. We are aware of studies in GaN27 and silicon carbides.28,29 The procedure is the same as that with empirical potentials: one initi­ates a series of cascades of low but increasing energy and follows the displacement of the accelerated atom. The threshold energy is reached as soon as the atom does not return to its initial position at the end of the cascade. Such calculations are very promising as empirical potentials are usually im­precise for the orders of energies and interatomic distances at stake in threshold energies. However, they should be done with care as most pseudopo­tentials and basis sets are designed to work for moderate interatomic distances, and bringing two atoms too close to each other may lead to spurious results unless the pseudopotentials are specifically designed.

1.08.3.3.2 Electronic stopping power Second, recent studies have been published in the ab initio calculations of the electronic stopping power for high-velocity atoms or ions. The framework best suited to address this issue is time-dependent DFT (TD-DFT). Two kinds of TD-DFT have been applied to stopping power studies so far.

The first approach relies on the linear response of the system to the charged particle. The key quan­tity here is the density-density response function that measures how the electronic density of the solid reacts to a change in the external charge density. This observable is usually represented in reciprocal space and frequency, so it can be confronted directly with energy loss measurements. The density — density response function describes the possible excitations of the solid that channel an energy trans­fer from the irradiating particle to the solid. Most noticeably the (imaginary part of the) function vanishes for an energy lower than the band gap and shows a peak around the plasma frequency. Integrating this function over momentum and energy transfers, one obtains the electronic stopping power. Campillo, Pitarke, Eguiluz, and Garcia have implemented this approach and applied to some simple solids, such as aluminum or silicon.3 — They showed that there is little difference between the usual approximations of TD-DFT: the random phase approximation, which means basically no exchange correlation included, or adiabatic LDA, which means that the exchange correlation is local in space and instantaneous in time. The influence of the band structure of the solid accounts for noticeable deviations from the homogeneous elec­tron gas model.

The second approach is more straightforward conceptually but more cumbersome technically. It proposes to simply monitor the slowing down of the charged irradiated particle in a large box in real space and real time. The response of the solid is hence not limited to the linear response: all orders are automat­ically included. However, the drawback is the size of the simulation box, which should be large enough to prevent interaction between the periodic images. Fol­lowing this approach, Pruneda and coworkers33 cal­culated the stopping power in a large band gap insulator, lithium fluoride, for small velocities of the impinging particle. In the small velocity regime, the nonlinear terms in the response are shown to be important.

Unfortunately, whatever the implementation of TD-DFT in use, the calculations always rely on very crude approximations for the exchange-correlation effects. The true exchange-correlation kernel (the second derivative of the exchange-correlation energy with respect to the density) is in principle nonlocal (it is indeed long ranged) and has memory. The use of novel approximations of the kernel was recently introduced by Barriga-Carrasco but for homogeneous electron gas only.34,35

LJ Phase Diagram

In practice, pair potentials are cut off at a certain range, which can have a surprising effect on stability as shown in Figure 1

While the LJ fluid is very well studied, the finite temperature crystal structure has only recently been resolved. The problem is that the fcc and hcp structures are extremely close in energy (see Figure 1(b)), so the entropy must be calculated extremely accurately.

Подпись: Figure 1 Energy difference between hcp and fcc for the Lennard-Jones potential at 0 K, as a function of cutoff (rc) with either simple truncation or with the potential shifted to remove the energy discontinuity at the cutoff. Without truncation the difference is 0.0008695 e, with hcp more stable.
This has been done by Jackson using ‘lattice — switch Monte Carlo3 (see Figure 2). The equivalent phase diagram for the Morse potential remains unsolved. The LJ potential has been used extensively for fcc materials, and it still comes as a surprise to many researchers that fcc is not the ground state.

Influence of free surfaces

The rationale for investigating the impact of free surfaces on cascade evolution is the existence of an influential body of experimental data provided by experiments in which thin foils are irradiated by high-energy electrons and/or heavy ions.98-106 In most cases, the experimental observations are carried out in situ by TEM and the results of MD simulations are in general agreement with the data from these experiments. For example, some material-to-material differences observed in the MD simulations, such as differences in in-cascade clustering between bcc iron and fcc copper, also appear in the experimental data.59,107,108 However, the yield of large point defect clusters in the simulations is lower than would be expected from the thin foil irradiations, particularly for vacancy clusters. It is desirable to investigate the source of this difference because of the influence this data has on our understanding of cascade damage formation.

Both simulations81,97,109,110 and experimental work105,106 indicate that the presence of a nearby free surface can influence primary damage formation. For example, interesting effects of foil thickness
have been observed in some experiments.105 Unlike cascades in bulk material, which produce vacan­cies and interstitials in equal numbers, the number of surviving vacancies in surface-influenced cascades can exceed the number of interstitials because of interstitial transport to the surface. This could lead to the formation of larger vacancy clusters and account for the differences in visible defect yield observed between the results of MD cascade simula­tions conducted in bulk material and the thin-film, in situ experiments. Initial modeling work reported by Nordlund and coworkers81 and Ghaly and Averback109 demonstrated the nature of effects that could occur, and Stoller and coworkers97,100 subse­quently conducted a study involving a larger number of simulations at 10 and 20keV to determine the magnitude of the effects.

To carry out the simulations,97,100 a free surface was created on one face of a cubic simulation cell containing 250 000 atom sites. Atoms with sufficient kinetic energy to be ejected from the free surface (sputtered) were frozen in place just above the surface. Periodic boundary conditions are otherwise imposed. Two sets of nine 100 K simulations at 10keV were carried out to evaluate the effect of the free surface on cascade evolution. In one case, all the PKAs selected were surface atoms and, in the other, PKA were chosen from the atom layer 10a0 below the free surface. The PKA in eight 20keV, 100 K simulations were all surface atoms. Several PKA directions were used, with each of these directions slightly more than 10° off the [001] surface normal.

Подпись:

Подпись: x
Подпись: 1.588E-11
Подпись: Figure 25 Defect evolution in typical 10 keV cascade initiated by a surface atom: (a) peak damage state at ~1.1 ps, and (b) final damage state at ~15 ps.

Figure 25 provides a representative example of a cascade initiated at the free surface. The peak damage state at 1.1 ps is shown in (a), with the final damage state at ~15 ps shown in (b). The large number of apparent vacancies and interstitials in

Figure 25(a) is due to the pressure wave from the cascade reaching the free surface. With the constrain­ing force of the missing atoms removed, this pressure wave is able to displace the near-surface atoms by more than 0.3a0, which is the criterion used to choose atom locations to be displayed. As mentioned above, a similar pressure wave occurs in bulk cascades, making the maximum number of displaced atoms much greater than the final number of displacements. Most of these displacements are short-lived, as shown in Figure 26, in which the time dependence of the defect popula­tion is shown for three typical bulk cascades, one surface-initiated cascade, and one cascade initiated 10a0 below the surface. The effect of the pressure wave persists longer in surface-influenced cascades, and may contribute to stable defect formation.

The number of surviving point defects (normalized to NRT displacements) is shown in Figure 2 7 for both bulk and surface cascades, with error bars indicating the standard error of the mean. The results are simi­lar at 10 and 20keV. Stable interstitial production in surface cascades is not significantly different than in bulk cascades; the mean value is slightly lower for the 10 keV surface cascades and slightly higher for the 20 keV case. However, there is a substantial increase in the number of stable vacancies produced, and the change is clearly significant. It is particularly worth noting that the number of surviving interstitials and vacancies is no longer equal for cascades initiated at the surface because interstitials can be lost by sputter­ing or the diffusion of interstitials and small glissile
interstitial clusters to the surface. Reducing the num­ber of interstitials leads to a greater number of sur­viving vacancies, as less recombination can occur.

In-cascade clustering of interstitials is also rela­tively unchanged in the surface cascades (e. g., see Figures 4 and 5 in Stoller11 ). The effect on in­cascade vacancy clustering was more substantial. The vacancy clustering fraction per NRT (based on the fourth NN criterion discussed above) increased from ~-0.15 to 0.18 at 10 keV and from ^0.15 to 0.25 at 20 keV. Moreover, the vacancy cluster size distribu­tions changed dramatically, with larger clusters pro­duced in the surface cascades. The free surface effect on the vacancy cluster size distributions obtained at 20 keV bulk is shown in Figure 28. The largest vacancy cluster observed in the bulk cascades con­tained only six vacancies, while the surface cascades had clusters as large as 21 vacancies. This latter size is near the limit of visibility in TEM, with a diameter of almost 1.5 nm. Overall, these results imply that cascade defect production in bulk material is different from that observed in situ using TEM. More research such as that by Calder and coworkers111 is required to fully assess these phenomena, particularly for higher cascade energies, in order to improve the ability to make quantitative comparisons between simulations and experiments.

Sink strength of voids

Consider 3D diffusion of mobile defects near a spherical cavity of radius R, which is embedded in a lossy-medium of the sink strength k2:

G — k2 D(C — Ceq) — VJ = 0 [46]

where Ceq is the thermal-equilibrium concentration of mobile defects and the defect flux is

J = — d(vC + kCrVtf) [47]

Here, D is the diffusion coefficient, U is the interac­tion energy of the defect with the void, kB is the Boltzmann constant, and T the absolute temperature. The boundary conditions for the defect concentra­tion, C, at the void surface and at infinity are

C(R) = Ceq [48]

G

C1 = Ceq + _ [49]

k2D

Equation [49] follows from eqn [46] and the require­ment that the gradients vanish at large distances. Here, all other sinks in the system, voids, dislocations, etc. are considered in the MFA and contribute to the total sink strength k2. This procedure is self­consistent.

The interaction energy of a defect with the void in eqn [47] is small and usually neglected. The solution of eqn [46] for a void located at the origin of the coordinate system, r = 0, is then

C(r) = Ceq + (c1 — Ceq) 1 — Rexp[-k(r — r)] [50]

 

image787

In order to predict the evolution of mobile PDs and their impact on immobile defects, one needs to know the sink strength of different defects for vacancies and SIAs and the rate of their mutual recombination. The reaction kinetics of 3D migrating defects is con­sidered to be of the second order because the rate equations contain terms with defect concentrations to the second power.40 An important property of such kinetics is that the leading term in the sink strength of any individual defect depends on the characteristics of this defect only. Thus,

N

ka = k2j [45]

j =1

 

The total defect flux, I, through the void surface S = 4nR2 is given by

I = — SJ(R) = kC(R)D(C1 — Ceq) [51]

where the void sink strength is

kC (R) = 4pR(1 + kR) [52]

 

The sink strength of all voids in the system is obtained by integrating over the SDF, f (R):

 

where a = v, i and N is the total number of sinks per unit volume. For example, the total sink strength of an ensemble of voids of the same radius, R, is equal to ka = Nka(R). The individual sink strength such as a void or a dislocation loop may be obtained from a solution to the PD diffusion equation. In the

 

image788

k2

kC =

 

[53]

 

where NC = dRf(R) is the void number density, (R) is the void mean radius and (R2) is the mean radius squared. Typically, k2 « 1014m-2, that is,

 

image789

k—1 « 100 nm, while the void radii are much smaller, so that one can omit the term proportional to the radius squared:

kC = 4k(R)Nc [54]

Equation [52] is derived by neglecting the interaction between the void and mobile defect. There is a differ­ence between the interaction of SIAs and vacancies with voids due to differences in the corresponding dilatation volumes. As a result, the void capture radius for an SIA is slightly larger than that for a vacancy (see, e. g., Golubov and Minashin94). However, this difference is usually negligible compared to that for an edge dislocation, which is described below.

 

where p, is the dislocation density and Zd the cap­ture efficiency. The capture efficiencies for vacancies and SIAs, Za and Za, are different because of the difference in their dilatation volumes (see eqn [56])

 

2p

ln(1/kR°[)

 

Zd

 

[60]

 

where a = v, i and

 

mb l + n eg

3p 1 — n 4kB T

 

image790

[61]

 

DOa

 

The dilatation volume of SIAs is larger than that of vacancies, hence RD > RD and the absorption rate of dislocations is higher for SIAs: Zd > Z,. This is the reason for void swelling, which is shown below in Section 1.13.5.2.1. A more detailed analysis of the sink strengths of dislocations and voids for 3D diffus­ing PDs can be found in a recent paper by Wolfer.96

 

image791

General Principles and Applications of PF Modeling

The first step in PF modeling lies in choosing and defining the fields of interest. These continuous variables are functions of space and time, and they are in most cases scalar fields, such as temperature, or the concentration of some chemical species of interest. In systems with solid-liquid interfaces, a phenomenological field variable is introduced in such a way that it varies continuously from 0 to 1 as one goes from a fully solid to a fully liquid phase. Multidimensional fields can be used as well, for instance, to describe the local composition of a multicomponent alloy, the local degree of chemical order, or the local crystallographic orientation of grains. These multidimensional fields may transform like vectors under symmetry operations, thus leading to a vectorial representation of the system and ten — sorial expressions for mobilities (as will be discussed later), but there are cases for which the multidimen­sional fields cannot be reduced to vectors.27 In all cases, an averaging procedure is necessary to define continuous field variables for systems that are intrin­sically discrete at the atomic scale. Various averages can be used, including (1) a spatial average over representative volume elements, which will corre­spond to the cells used for evolving the PF variables; (2) a spatial and temporal average; or (3) a spatial and ensemble average. The spatial averaging method is used most often, although in many cases the exact conditions of the averaging procedure are not defined. Section 1.15.3 will cover a model where this coarse-graining is performed explicitly and rigor­ously. The last two averaging procedures are rarely explicitly invoked, although one oftheir advantages is that a smaller volume can be used for the spatial average, thanks to the additional averaging per­formed either in time or in the configuration space of a system ensemble.

Turning now to the kinetic equations used to describe the evolution of these field variables, an important distinction is whether the field variable is conserved or nonconserved. For the sake of simplic­ity, the following discussion focuses on alloy systems. Let us consider two simple examples, one where the field variable is the local composition, C(r, t), in a binary A-B alloy system, and a second example, also for a binary alloy system, but this time with a fixed composition and where chemical ordering takes place. The degree of chemical order is described by the field S(r, t). For the sake of convenience, one

Подпись: d C (r,t) dt Подпись:Подпись: MVПодпись: [2]Подпись: [3]Подпись: F{C(r, t)}Подпись: [4]Подпись: dVПодпись:Подпись: amay normalize that field such that S(r, t) = 0 corre­sponds to a fully disordered state and S(r, t) = ±1 to a fully ordered state. The first field variable C(r, t) is globally conserved — assuming here that the system of interest is not exchanging matter with its environ­ment. This imposes the constraint that the time evolu­tion of the field variable at r is balanced by the divergence of the flux of species exchanged between the representative volume centered on r and the remainder of the system:

Чг =-VJ (V) [1]

One then makes use of linear response theory in the context of thermodynamics of irreversible processes28 to linearly relate the flux J(r, t) to the driving force responsible for this flux. Here, this driving force is the gradient of the chemical potential m(r, t) = 8F/8C(r, t), where F is the free energy of the system for a compositional field given by C(r, t). The resulting evolution equation is thus

‘ 8 F 8C(r, t)

where M is a mobility coefficient. In contrast, for the nonconserved order parameter S(r, t), its evolution is directly related to the free energy change as S(r, t) varies, so that by making use of linear response theory again

d S (r, t) r 8 F

dt = ~L 8S(r, t)

where L is the mobility coefficient for the nonconserved field S(r, t). Two important consequences of eqns [2] and [3] are worth noticing. First, although all extrema of the system free energy (i. e., minima, maxima, saddle points) are stationary states, often in practice only the minima can be obtained at steady state due to numerical errors. Second, the stationary state reached from some initial state may not correspond to the absolute minimum of the free energy. In order to overcome this problem, noise can be added to transform these deterministic equations into stochastic (Langevin) equations, as will be discussed in Section 1.15.3.

Following the work of Cahn and coworkers2-4 and Landau and Ginzburg,6 the free energy F is decomposed into a homogeneous contribution and an heterogeneous contribution. Treating the inhomo­geneity contribution as a perturbation of a homoge­neous state, one finds that, in the limit of small amplitude and long wavelength for this perturbation, the lowest order correction to the homogeneous free
energy is proportional to the square of the gradient of the field variable. For instance, returning to the sim­ple example of an alloy described by the concentra­tion field C(r, t), the total free energy can be written as dV [f (C) + k(VC)2] where f (C) is the free-energy density of a homoge­neous alloy for the composition C, and к the gradient energy coefficient, which is positive for an alloy system with a positive heat of mixing. A similar expression can be used in the case of a nonconserved order parameter, for example, S(r, t), or more generally, in the case of an alloy described by nC conserved order parameters and nS nonconserved order parameters f (Cb C2 ••• CnC ; Sb S2 • •• SnS )

nS

+ Kp (VCp) + Vij ViSqVjSq

p=1 q=1

An implicit summation over the indices i and j is assumed in the last term of eqn [5]. The number of nonzero and independent gradient energy coefficients V j for the nonconserved order parameters is dictated by the symmetry of the ordered phase. Specific exam­ples, for instance for the L12 ordered structure, can be found in Braun eta/.27 and Wang eta/.29 The free energy can also be augmented to include other contributions, in particular those coming from elastic fields using the elasticity theory of multiphase coherent solids pio­neered by Khachaturyan,30 in the homogeneous mod­ulus case approximation. This makes it possible to take into account the effect of coherent strains imposed by phase transformations or by a second phase, for exam­ple, a substrate onto which a thin film is deposited.31

Two important interfacial quantities, the excess interface free energy and the interface width, can be derived from eqn [5] for a system at equilibrium. We follow here the derivation given by Cahn and Hilliard.2 Considering the case of a binary alloy where two phases may form, referred to as a and p, and with respective B atom concentrations Ca and Cp, the existence of an interface between these two phases results in an excess free energy a dv[f(c) + k(vc)2 — CmB — (1 — C)mA] [6] where mA and mB are the chemical potential of A and B species when the two phases a and p coexist at

Подпись: dV K(V С)2 V Подпись:Подпись: sПодпись: Cp 2 dC KVC ca Подпись:Подпись: [12]Подпись: 2С - 1 Cap Подпись:Подпись:equilibrium. At equilibrium, this excess free energy is minimum. A homogeneous free energy Dfreferenced to the equilibrium mixture of a — and p-phases is introduced as

Df(c)=f(c)-[CmB + (1 — c)mA]

= сMC) — mB] + (1 — c)[ma(c) — mA] [7]

(Note that the ‘D’ symbol in Df in eqn [7] does not refer to a Laplacian.) The variational derivative of this excess energy with respect to the concentration field is given by

8s dDf дк 2 r і

8c = 1>C — 2kDC — дс(гс) [8]

At equilibrium, the excess free energy s is minimum, and the concentration field must be such that 8s/8C = 0. Thus,

= 2kDC + (vc)2 = (k(VC)2) [9]

дС 3l ; дС

Equation [9] must hold locally for any value of the concentration field along the equilibrium profile join­ing the a — and p-phases, and this can only be satisfied if

k(VC)2 = Df [10]

for all values of C(r). It is interesting to note that eqn [10] means that the equilibrium concentration profile is such that, at any point on this profile, the homoge­neous and inhomogeneous contributions to the total free energy are equal. The interfacial excess free energy is thus given by 2 dV Df = 2

V

This last integral over the spatial coordinates can be rewritten as an integral over the concentration field. Assuming a one-dimensional (1D) system for simplicity,

Cp

2 dCy/KAf

Ca

In order to proceed further, it is necessary to assume a functional shape for the concentration profile or for the homogeneous free energy Df. Expanding the free energy near the critical point Tc yields a symmetric double-well potential for the homogeneous free energy,2 which we write here as 2С — 1

Cap
with Cap = Cp — Ca = 1 — 2Ca = 2Cp — 1. Using eqn [10], the equilibrium concentration profile can now be obtained:

ОД="Т +1 [14|

Integration along this equilibrium profile from eqn [11] yields the interfacial energy

s 3Cap у8 KA/max [15]

Furthermore, the width of the equilibrium profile we, which is defined as the length scale entering the argu­ment of the hyperbolic tangent function in eqn [14], is given by

w=crvлЕ [16]

In this conventional approach to PFMs, Dfmax and к are phenomenological coefficients. Equations [15] and [16] play an important role in assigning values to these coefficients for a specific alloy system. The excess interfacial energy may be known experimen­tally or it may be calculated separately, for instance by ab initio calculations.32 If the interfacial width we is also known, one can obtain Dfmax and к from an inverse solution of eqns [15] and [16] (note that Cap is given by the equilibrium phase diagram). Even if we is not known, values for Dfmax and к can be chosen to yield a prescribed value for s. In all cases, it is important to recognize that any microstructural fea­ture that develops during the simulations is expressed in units of we. At elevated temperatures, as T! Tc, s vanishes2 while we goes to infinity, and therefore, at high enough temperatures, interfaces are diffuse, thus meeting this essential requirement underlying the PF method.

The PF eqns [2] and [3] are usually solved numer­ically on a uniform mesh with an explicit time inte­gration, using periodic boundary conditions, when surface effects are not of interest. When the free energy contains an elastic energy contribution, it is quite advantageous to use semi-implicit Fourier — spectral algorithms (see Chen9 and Feng et a/.33 for details). Variable meshing can also be employed, in particular to better resolve interfaces when they tend to be sharp, for instance at low temperatures.

A few examples selected from the literature serve to illustrate the capacity of PFMs to successfully reproduce a wide range of phenomena. In particular, Khachaturyan30 and his collaborators34-38 proposed a microelasticity theory of multiphase coherent solids,
which has been widely used to include a strain energy in the overall free energy. A method for systems with strong elastic heterogeneity has been proposed by Hu and Chen,39 which includes higher order terms that are usually neglected in Khachaturyan’s approach. Figure 1 illustrates the anisotropic morphology of Al2Cu precipitates growing in an Al-rich matrix.3 Bulk-free energies were calculated using a mixed- space cluster expansion technique, with input from first-principle calculations for about 40 different ordered structures with full atomic relaxations. Interfacial energies were calculated at T = 0 K from first-principle calculations as well, using configura­tions where the Al-rich solid solution and the tetrag­onal 0′-Al2Cu coexist. For the elastic strain energy calculations, the elastic constants of 0′-Al2Cu were calculated ab initio. An important feature of this sys­tem is that both elastic and interfacial energies are strongly anisotropic, and the PF approach makes it possible to include these anisotropies. Furthermore, when the high-aspect-ratio 0′-phase forms, its growth kinetics will be anisotropic as well, which can be included in a phenomenological way by introducing a dependence of the mobility on the orientation of the precipitate-matrix interface. Figure 1 illustrates that these three anisotropies, interfacial, elastic, and kinetic, are required to reproduce the morphology of 0′ precipitates.

Figure 2 illustrates another effect of coherency stress on microstructural evolution, this time for an A1-L10 order-disorder transition in a Co-Pt alloy.40 The tetragonal distortion accompanying the ordering reaction leads to the formation of self-organized tweed patterns of coexisting (cubic) A1- and (tetrago­nal) L10-phases. As seen from Figure 2, the agreement between experimental and simulated microstructures is remarkable. Ni and Khachaturyan proposed recently that, in order to minimize elastic energy during trans­formations involving symmetry changes and lattice strain, a pseudospinodal decomposition is likely to take place, leading to 3D chessboard patterns.41

PF modeling has also been used extensively to study martensitic transformations,34-38,42-44 phase transformations in ferroelectrics45-57 (see also the recent review by Chen58 on that topic), transforma­tions in thin films,47,59-65 grain growth and recrystal­lization,66-81 and microstructural evolution in the presence of cracks or voids.82-84 A recent extension of PFMs has been the inclusion of dislocations in the models,85-87 by taking advantage of the equivalence between dislocation loops and coherent misfitting platelet inclusions.88 This approach has been applied, for instance, to study the interaction between moving dislocations and solute atoms,89 or to study the influ­ence of dislocation arrays on spinodal decomposition in thin films.61 Rodney et a/.87 have pointed out,

image206

Подпись: IsotropicПодпись: Interface onlyimage929Strain only

(c)

image930
Int +strain + kinetics Experiment

(b) (c)

Подпись:Time

(d) (e) (f)

Figure 2 Comparison between transmission electron microscopy experimental observations (a-c) and phase field modeling (d-f) of formation of chessboard pattern in Co395Pt605 cooled from 1023 K to (a) 963 K, (b) 923 K, and (c) 873 K. The scale bar corresponds to 30 nm. Reproduced from Le Bouar, Y.; Loiseau, A.; Khachaturyan, A. G. Acta Mater. 1998,46(8), 2777-2788.

image932however, that the artificially wide dislocation cores required by the above approach lead to weak short — range interactions. These authors have introduced a different PFM for dislocations, which allows for nar­row dislocation cores. As an illustration of that model, Figure 3 shows the development of dislocation loops and their interaction with hard precipitates in a 3D y/y’ single crystal. It is interesting to note that dislo­cation loop initially expands by gliding in the soft y channels, until the local stresses are large enough for the dislocation to shear the hard y’-phase.

The above presentation of the PF equations leaves certain questions open. First, the maximum homoge­neous free energy difference А/max involves the free energy of the unstable state separating the two minima at Ca and Cp. It is thus been questioned90 whether this quantity can be rigorously defined from thermody­namic principles. If one employs mean field techni­ques such as the cluster variation method (CVM)91-95 to derive the homogeneous free energy of an alloy, А/max is in fact very sensitive to the approximation used, and generally decreases as the size of the largest cluster used in the CVM increases.96 Kikuchi97-99 has argued that, in order to resolve this paradox, А/ should not be considered as the free energy of any homogeneous state, but that it should be understood as the local contribution to the free energy of the system along the equilibrium composition profile.

The second set of questions relates to the gradient energy coefficients к and q in eqn [5]. In many applications of PF modeling, these coefficients are taken as phenomenological constants that can be adjusted at will, as long as the microstructures are scaled in units of к1/2 or q1/2. Such an approach, however, is problematic for many reasons. First, when one scalar field variable is employed, for instance C(r/), a regular solution model,2, 0 or equiv­alently a Bragg-Williams approximation,101 estab­lishes that the gradient energy coefficient is not arbitrary but that it is directly proportional to the interaction energy between atoms, that is, to the heat of mixing of the alloy. Furthermore, in the most general case, к should in fact be composition and temperature dependent. Starting from an atomistic model, rigorous calculations of к are possible by monitoring the intensity of composition fluctuations as a function of their wave vector, and using the fluctuation-dissipation theorem.1 0 In the case of a simple Ising-like binary alloy, it is observed that к varies as C(1 — C), where C is the local composition of the alloy. 0 Furthermore, when more than one field variable is employed, care should be taken to consider all possible contributions of field heterogeneities to the free energy of the system, as the different fields may be coupled. Symmetry considerations are important to identify the nonvanishing terms, but it

may remain challenging to assign values to these nonvanishing terms that are consistent with the ther­modynamics of the alloy considered.

Another important point is that interfacial ener­gies are in general anisotropic. In order to obtain realistic morphological evolution, it is often impor­tant, and sometimes even absolutely necessary, to include this anisotropy, for example, in the modeling of dendritic solidification. The symmetry of the mesh chosen for numerically solving the PF equations introduces interfacial anisotropy but in an unphysical and uncontrolled way. One possible approach to introduce interfacial anisotropy is to let к vary with the local orientation of the interface with respect to crystallographic directions.11, Another approach is to rely on symmetry constraints27,30 to determine the number of independent coefficients in a general expression of the inhomogeneity term, see eqn [5]. In both approaches, the different coefficients entering
the interfacial anisotropy can be fitted to experiments or to atomistic simulations.

Let us now return to the mobility coefficients M and L introduced in eqns [2] and [3]. For the sake of simplicity, many PF calculations are performed while assigning an arbitrary constant value to these coeffi­cients. An improvement can be made by relating the mobility to a diffusion coefficient. In the case of M, for instance, in order to make eqn [1] consistent with Fick’s second law for an ideal binary alloy sys­tem, one should choose

C)

-Б where C is the average solute concentration and D the interdiffusion coefficient. In both cases, the simulated times are expressed in arbitrary units of M-1 or L — , thus precluding a direct connection with experimental kinetics. This problem is also directly related to the lack of absolute physical length scales in these simulations. Moreover, using a 1D Bragg-Williams model composed of atomic planes, Martin101 showed that M is not a constant but is in fact a function of the local composition along the equilibrium profile. A complete connection between atomistic dynamics and M will be made in Section 1.15.3. Similar to the discussion on coupling between various fields for the gradient energy terms, kinetic coupling is also expected in general. The kinetic couplings between composi­tion (a conserved order parameter) and chemical ordering (a nonconserved order parameter) are revealed by including sublattices into Martin’s 1D model and deriving the macroscopic evolution of the fields from the microscopic dynamics. In that case, atoms jump between adjacent planes.102,103 As a result, instead of the mere superposition of eqns [2] and [3], the kinetic evolution of coupled con­centration and chemical order in a binary alloy is given by

dC(r, t) dt

+ V( M-’V( Щ5))]

+ V(L-’V

+ V L3 V where L1 is a mobility coefficient, and L2, L3, M1, and M2 are second-rank mobility tensors, since they

image422

relate diffusional fluxes (vectors) to chemical potential gradients (vectors). In the case of cubic crystalline phases, second-rank tensors reduce to scalars, but in many ordering reactions, noncubic phases form, thus leading to anisotropic mobility. Vaks and coworkers104 have also derived PFMs for simultaneous ordering and decomposition starting from microscopic models. These works, however, illustrate the fact that it would be quite difficult, especially for multidimensional field variables, to assign correct values to the kinetic coefficients for a given alloy system by relying solely on a phenomeno­logical approach.