Etude de la diffusion du carbone dans le zirconium et la zircone en volume des gaines de combustible usées par simulations multi-échelles

Dans le cadre de la gestion des déchets nucléaires, les coques et embouts des gaines de combustibles sont coupés, compactés et empilés dans les colis CSD-C (Colis Standards de Déchets Compactés). Actuellement entreposés à la Hague, ces déchets seront stockés en milieu géologique profond dans des alvéoles en béton. Ces déchets MA-VL contiennent des RN d'intérêt dont le carbone-14 qui provient de l'activation neutronique de l'azote-14 et de l'oxygène-17 présents dans le Zircaloy oxydé. L'objectif de cette thèse, qui est menée en collaboration avec EDF et AREVA, est de répondre à la question du devenir et du comportement du carbone-14 dans ZrO 2 et Zr-métal. Cette thèse consistera à utiliser la simulation numérique, au moyen de la Théorie de la Fonctionnelle de la Densité (DFT), afin de modéliser le comportement du carbone-14 dans les CSD-C.

Les simulations sont effectuées avec une approche multi-échelles : 1) A l'échelle atomique, nous avons optimiséles modèles atomiques permettant dereprésenter la phase hexagonale de Zr et les différentes phases de ZrO 2 . Par l'identification des différents sites d'insertion possibles pour l'atome de carbone, deux types de sites interstitiels possibles ont été caractérisés: le site Basal Tétragonal déformé noté, site BT_d, et le site Octahedrique, noté que site O. Le site BT_d est moins stable que le site O par 1,45 eV. Pour C dans la phase monoclinique, huit types de sites interstitiels possibles pour l'atome de carbone ont été caractérisés. Les énergies de formation pour l'atome de carbone en substitution à un atome d'oxygène ou de zirconium ont été calculées pour comparer les stabilités relatives de l'atome carbone. La modélisation des différents chemins d'énergie minimale, associés à la migration d'un atome de carbone d'un site à un autre, nous permettra d'obtenir la trajectoire du carbone par mécanismes interstitiels et lacunaires.

2) A l'échelle macroscopique, les simulations par l'approche Kinetic Monte-Carlo (KMC) dans ZrO 2 monoclinique pure et Zr pure ont été réalisées en utilisant les probabilités de saut. A partir des simulations KMC à des températures différentes, la trajectoire du carbone ainsi que le coefficient de diffusion (en fonction de la température) ont été calculés. Il est constaté que l'énergie d'activation effective est 2.02 eV pour la diffusion du carbone dans zirconium, ce qui correspond au chemin de migration avec l'énergie la plus basse: site O siteBT_d siteBT_d (n fois) site O siteBT_d ... (par rapport au chemin site O site O). Le coefficient de diffusion du carbone dans les cristaux de Zr pur, à la température de stockage (environ 50 °C) est calculé pour être très faible (~ 10 -34 cm 2 /s). Par conséquent, ces premières investigations montrent qu'il sera très difficile pour un atome de carbone de diffuser dans un cristal parfait de Zr dans les conditions de stockage. D'après les simulations macroscopiques Monte-Carlo cinétique pour l'atome carbone dans la zircone monoclinique, on a constaté qu'il est très difficile pour cet atome de diffuser par un mécanisme interstitiel (avec une énergie d'activation effective de 1,19 eV). Le coefficient de diffusion du carbone dans la phase Zr monoclinique, à la température de stockage (environ 50 °C) est calculée pour être très faible (~ 10 -22 cm 2 /s).A titre d'exemple, un atome de carbone présent au centre d'une couche d'oxyde de 50µm aura besoin de 10 11 ans pour la quitter. defects . ............................................................................................................................ O1, .................................................................................................................................................

II.1.A Vacancy of

II.2 Electronic

General Introduction Introduction

In a pressurized water reactor (PWR), the core contains fuel assemblies that are cooled by water circulated using powered pumps. PWRs contain between 150-200 fuel assemblies. In each fuel assembly, there are about 264 rods placed in a square array. For the 4 meter tall rod, the fuel consists of uranium dioxide and MOX (uranium and plutonium oxide) which are loaded in 272 pellets. The fuel rods' claddings are made of zirconium alloys (such as zircaloy-4) which is principally zirconium. The composition of zircaloy-4 is controlled with the weight percent listed in Table 1. The zirconium has a hexagonal close packed structure (α -Zr) which is stable until 800 °C. The composition elements' characterization and their crystallographic structures in the zirconium matrix have been described in the case of Zy-4 [X. Y. Meng 1985].

Element

Content (wt %) Zr 97.37-98.17 The formation of oxide layer and Since there is water outside the claddings, the oxidation reaction with the zirconium forms an oxide layer which is zirconia (tetragonal and

The oxide process of the zirconium claddings has been studied in several projects [P. Bouvier 2000[START_REF] Gosmain | [END_REF], C. Roy 1970[START_REF] Garzarolli | [END_REF], J. S. Moyal Garcia 1999]. The thickness of the oxide layer is function of time [L. [START_REF] Gosmain | [END_REF]. For the 570 µm thickness cladding, the oxide layer is about 50 µm after 6 nuclear fuel cycles [ et al., 2005] (Figure 2). At the storage temperature ( structure is monoclinic zirconia which will be chosen in this work for later calculations [E. C. Subarao 1974, H. Boysen 1991]. The formation of oxide layer and 14 C

Since there is water outside the claddings, the oxidation reaction with the zirconium forms an oxide layer which is zirconia (tetragonal and monoclinic phases).

→ 2

The oxide process of the zirconium claddings has been studied in several projects [P. Bouvier 2000[START_REF] Gosmain | [END_REF], C. Roy 1970[START_REF] Garzarolli | [END_REF], J. S. Moyal Garcia 1999]. The thickness of the oxide layer is function of time [L. [START_REF] Gosmain | [END_REF]. For the 570 µm thickness cladding, the oxide layer is about 50 µm after 6 nuclear fuel cycles [ (Figure 2). At the storage temperature (max around 50° C), the stable zirconia structure is monoclinic zirconia which will be chosen in this work for later calculations [E. C. Subarao 1974, H. Boysen 1991]. Since there is water outside the claddings, the oxidation reaction with the zirconium forms

(1)

The oxide process of the zirconium claddings has been studied in several projects [P. Bouvier 2000[START_REF] Gosmain | [END_REF], C. Roy 1970[START_REF] Garzarolli | [END_REF], J. S. Moyal 2000, E. A. Garcia 1999]. The thickness of the oxide layer is function of time [L. [START_REF] Gosmain | [END_REF]. For the 570 µm thickness cladding, the oxide layer is about 50 µm after 6 nuclear fuel cycles [Bossis around 50° C), the stable zirconia structure is monoclinic zirconia which will be chosen in this work for later calculations [E. C.

contain numerous activation products such as 14 C. It comes from curies/g) and 14 N (8.7×10 -3

(2)

(3)

ZrO 2
This neutron irradiation during the reaction at high temperature causes the phase transition from tetragonal to monoclinic zirconia [M. C. Wittels 1956, M. C. Wittels 1959]. Moreover, it induces some defects in the zirconia which has been proved by simulation studies [D. Simeone 2000, D. Simeone 2002, C. G. Mougel 2000, and X. Iltis 1995]. In the oxide layer, the O vacancy concentration is controlled to be as low as possible by decrease the tin level [A.

M. Garde 1994].

The process of nuclear waste treatment

In France, the spent nuclear assemblies are sent to AREVA La Hague to perform the nuclear waste treatment. After a cooling period, the claddings are sheared, and dipped in a dissolving bath which is mainly concentrated hot nitric acid to dissolve and extract the remaining fuel and fission products. Then the claddings are cut (Figure 3 The CSD-C are intended to be stored in the cells of Intermediary Level Waste Long Lived (ILW-LL) deep storage. For this, four CSD-C are arranged in a concrete disposal packages with enhanced holding capacity (Figure 4).

The importance to study the carbon diffusion in the claddings

In order to demonstrate the safety and the CSD-C behavior in deep geological disposal, the behavior of all the radionuclides must be quantified.

The behavior of radionuclides (RN) contained in zirconia and zirconium is currently poorly characterized. Of penalizing manner, ANDRA takes no retention performance of these materials in its safety analysis. In the release model of ANDRA, 100% of the 14 C content in waste is considered to be released upon arrival in the storage cell (term major source during the reversibility period for releases -100 years). Furthermore, it could end 14 C in organic form (CO 2 or CH 4 ) with the arrival of the water in the storage thereby promoting migration. In the absence of data, ANDRA considers the entirety of 14 C content in waste is released in organic form.

Some studies are required to improve the waste package storage behavior model from the main family ILW-LL generated by the activities of EDF and AREVA NC (CSD-C) and thus limit its impact on the long-term safe storage.

Regarding 14 C, experimental analysis of used cladding is complicated because of its activity.

However, it is difficult to represent the used cladding conditions in laboratory. The introduction/implantation of carbon without generating non representative defects is difficult.

Moreover, the measure of carbon diffusion as well as zirconia dissolution are not easy because of the slow kinetic of these processes. Consequently, simulations could provide some insights and allow to better quantify the carbon behavior in zirconium and zirconia.

The aim of this thesis is to characterize, by the mean of a multi-scale approach based on ab initio calculations and kinetic Monte Carlo method, the carbon diffusion behavior in bulk zirconium and zirconia to help the estimation of the waste safety of the fuel rods' claddings in the storage condition.

Manuscript organization

To study the carbon diffusion behavior in the used fuel claddings, supercells of Zr and ZrO 2 were first optimized to then perform the carbon diffusion in the two structures.

The first chapter presents the simulation methods used in this work: ab initio calculations in the framework of the Density Functional Theory (DFT) to determine formation energies of different defects and migration energies during atomic jumps, and the Kinetic Monte Carlo (KMC) method at the macroscopic scale to simulate diffusion coefficient as a function of temperatures in order to predict diffusion coefficients at storage temperature.

In the second chapter, the zirconium metal is considered, the properties of carbon in interstitial positions are investigated, and the carbon diffusion by interstitial mechanism is calculated with the multi-scale approach. The vacancy effect on the diffusion is predicted by using the binding energies.

The monoclinic zirconia of the oxide layer is studied in the third chapter. All the possible interstitial sites of carbon atoms are identified with their formation energies. The atomic jumps by interstitial mechanism are performed and the carbon diffusion coefficient is calculated with the KMC approach. Then, the vacancy effect on the diffusion is also considered. The carbon diffusion by vacancy mechanism and its competition with selfdiffusion of oxygen atoms are studied. The density of states (DOS) calculations are performed to help a better understanding of the different defects and their effect on the carbon diffusion.

Finally, a general conclusion will give an overview of the carbon diffusion behavior in bulk of the two structures with and without vacancies by different mechanisms.

[D. Sherrill, Physical Review Letter, 3, 1976(1959) Chapter 1

Methodologies

Introduction of quantum mechanics

Ab initio means "from the beginning" or "from first principles". Ab initio quantum chemistry distinguishes itself from other computational methods in that it is based solely on established laws of nature: quantum mechanics.

The theoretical foundation for computational chemistry is the time-independent Schrodinger wave equation.

I The Schrödinger equation and Born-Oppenheimer approximation

In quantum mechanics, Schrödinger equation formulated by Erwin Schrodinger in 1926 [E. Nelson 1966, W. Kohn 1954, A. D. Mclachlan 1964] is used to calculate the system state. The ultimate goal of most approaches in solid state physics and quantum chemistry is the solution of the time-independent, non-relativistic Schrödinger equation although some progress has been made in dealing with fully relativistic Dirac equation. The energy can be expressed by this formula:

= (1)
where:

• H is its Hamiltonian which is associated with the observable energy.

• Ψ is the wave function. It is a function of the positions of nuclides and electrons.

• E is the total energy of the system. It is a scalar.

The wave equation is a postulate of quantum mechanics. H contains all the terms that contribute to the energy of a system:

= (2)
For a solid made of nucleons with positive charges and electrons with negative charges, the energy of the system is given by:

= ∑ - ℏ ∇ ! " ! # $% -∑ ℏ ∇ & " & ' (% ∑ ) & ) * + , &* ' (% (-. ∑ + / !0 # $% $-1 -∑ ∑ ) & + , & ! ' (% # $%
(3) where:

• ∑ - ℏ ∇ ! " 2 # $%
is kinetic energy of electrons noted as T e;

• -∑ ℏ ∇ & " & ' (% is kinetic energy of nucleons noted as T N; • ∑ ) & ) * + , &* ' (% (-.
is potential energy nucleon / nucleon noted as V NN ;

• ∑ + / !0 # $% $-1
is potential energy nucleon / nucleon noted as V ee ;

• ∑ ∑ ) & + , & ! ' (% # $%
is potential energy nucleon / electron noted as V Ne .

The ab initio calculation can solve the equation when it is time independent without knowing the values of specific parameters.

It is usual to invoke the so-called adiabatic or Born-Oppenheimer approximation [N. C. Handy 1986, O. L. Polyansky 1997, E. F. Valeev 2003]. Since the nucleon mass is much larger than the electron mass (mA >>mi, the ratio of the electron mass to the nuclei mass mi/mA is less than 10 -5 ), the electrons move much faster than the nuclei. Hence, the electrons can be thought as particles responding instantaneously to the motion of the nuclei without requiring a finite relaxation time. Alternatively, during the time of a cycle of electrons motion, the change in nuclear configuration is negligible. So, it can be assumed that the electrons are in the adiabatic equilibrium of their ground state with respect to the positions of nuclei at all the times. The approximation makes it possible to separate the wave function by nucleon wave function and electron wave function. Electrons adapt themselves thus almost immediately to the position of cores. After Born-Oppenheimer approximation, the total wave function is described as (4) and (5): 7 $ , 9 ( : = 7 $ : * 79 ( :

(4)

+ 7 $ : = + 7 $ :

(5)

So it can be solved when the nuclear cores were thought to be fixed and the electrons moved around them. However, Eq 5 is still complicated to solve, even for simple atoms or molecules.

So further simplifications are introduced to make the total energy calculations performed accurately and efficiently. These include the Hartree-Fock (HF) approximation or the density functional theory (DFT) method to perform the electron-electron interactions and the pseudopotential theory to deal with the electron-ion interactions.

II The density functional theory (DFT)

Density functional theory (DFT) [BOOK E. K. U. Gross 1995, R. G. Parr 1983, M. Orio 2009[START_REF] Car | [END_REF]] is a computational quantum mechanical modeling method used in physics, chemistry and materials science to investigate the electronic structure (the ground state) of many-body systems, in particular atoms, molecules, and the condensed phases. With this theory, the properties of a many-electron system can be determined by using functional spatially dependent electron density. The DFT simulation is presently one of the most successful approaches to compute the electronic structure of matter.

Its applicability ranges from atoms, molecules and solids. In its original formulation, the density functional theory provides the ground state properties of a system, and the electron density plays a key role.

III Kohn-Sham method

To simplify the Hamiltonian, there are different methods. One of them is proposed by Kohn andSham in 1964 [W. Kohn 1965]. The electronic density ρ(r) can describe all the characters of the total energy. When the density is in its fundamental state ρ(r) = ρ 0 (r), the energy of the system is the smallest and corresponds to the fundamental state. The system receives the external potential energy from the nucleons, and the total energy of the electron clouds is calculated by Eq.6:

<=> = <=> ++ <=> ?@ <=> A +?B 7 :=7 :Cᶟ (6)
where:

• <=> is the kinetic energy of the electron gas without interaction;

• ++ <=> is the classical interaction coulomb electron-electron term;

• ?@ <=> is the exchange and correlation term; • A +?B 7 :=7 :Cᶟ is the interaction of electrons with the external potential created by the nucleons.

IV Generalized Gradient Approximation (GGA)

As discussed above, the Kohn-Sham equation allows an exact treatment of most contributions to the electronic energy including the major fraction of the kinetic energy. So far, the DFT theory is exact. All the calculations were performed within the framework of the generalized gradient approximation (GGA) [J. P. Perdew 1992, J. P. Perdew 1996, Z. Wu 2006]. All remaining uncertainty about the universal functional are folded into the exchangecorrelation functional ?@ <=>. By comparison with the Hartree-Fock effective Hamiltonians, this term must contain the exchange energy, correlation energy and the difference between the kinetic energies of non-interacting and interacting systems. To generate the exchangecorrelation functional, approximations are necessary to be proposed. There are several methods to calculate the exchange and correlation term, such as LDA (Local Density Approximation) and GGA (Generalized Gradient Approximation). In our calculation, the GGA approximation was used.The first logical step to go beyond LDA is the use of not only the information about the density at a particular point ρ(r), but to supplement the density with information about the gradient of the charge density, ∇ρ(r), in order to account for the nonhomogeneity of the true electron density. So the exchange correlation term becomes a function of ρ and ∆ρ.

?@ <=> = A E ?@ 7 :<7=7 : + F=7 :>C G (7)

It provides a larger variation than LDA. In the later calculation, we use the functional PW-91 (exchange and correlation of Perdew and Wang) [J. P. Perdew 1992, J. P. Perdew 1996, K. Burke 1998]. GGA's approximations have reduced the LDA errors of atomization energies of standard set of small molecules by a factor 3-5. This improved accuracy has made DFT a significant component of quantum chemistry.

V Plane wave notions

According to Kohn-Sham method, the electronic density ρ(r) can describe all the characters of the total energy. To describe the electronic density, the Plane wave basis sets can be used.

In the periodic simulation, a unit cell of atoms is used to model the system by replicating the cell in all dimensions, using periodic approach. When the system has the force convergence, its total energy is locally minimized.

Because the system is periodic, according to Bloch's theorem, the mono electronic wavefunction in a regular repeating lattice (periodic boundary conditions) can be given by:

H 7 I: = $H JI .L JJJI * M H JI 7 I: (8)
where:

• H JI 7 I: is the electronic wave function at a point R by a direct vector; • N JI is a wave vector set of three quantum number (n x, n y , n z );

• $H JI ., JJJI is the phase factor.

• I is the position.

For any lattice vector 9 JI of the simulation cell, there is

H O I + 9 JI P = $H JI ., JJJI * H JI 7 I: (9)
where:

• M H JI 7 I + 9 JI :=M H JI 7 I:

If9 JI is a replication vector, the electronic wave function must have the same magnitude in the two places, but possibly offset by a plane wave phase factor $H JI ., JJJI . The quantum numbers are essentially continues in an infinite solid. But in practice, we use a cut-off energy to restrict the set to a sphere in reciprocal space. When E cut-off is too small, the wave function couldn't correctly describe the electronic density. But if E cut-off is too large, there would increase a lot of calculation time.

It should make an infinitive plane wave to produce the wave function exactly. In practice, it is not possible. The E cut-off limits the number of plane wave:

ℏ²| H JI ²ST I ²⃓ " < @WBXYZZ (10) where:

• m: mass of an electron;

• N JI is a wave vector set of three quantum number (n x, n y , n z );

• [ I are the reciprocal lattice vectors defined by the condition $T I ., JJJI =1.

VI Pseudopotentials

A pseudopotential or effective potential description of real systems.

The projected augmented wave (PAW) method [P. E. Blochl 1994, G. Kresse 1999], as implanted in the Vienna Ab initio

In this approximation, core electrons are replaced by a pseudopotential V expressed by the wave function cut-off radius r is used to separate the valence electrons from the core electrons. The VASP calculation will be only applied on the valence electrons. First derived from an atomic reference state, requiring that the pseudo and all eigenstates have the same energies and amplitude (and thus density) outside a chosen core cut-off radius r c . When r is increased to r the total energy of the atom.

The exact all electron wave function is replaced by a suitable "soft" pseudo wave function -Polynomials (Troullier and Martins (Ref N. Troullier 1993)) effective potential is used as an approximation for the simplified

The projected augmented wave (PAW) method [P. E. Blochl 1994, G. Kresse 1999], as

Ab initio Simulation Package (VASP) is used here.

In this approximation, core electrons are replaced by a pseudopotential V y the wave function Ψ pseudo. The cut-off energy which corresponds to the different off radius r is used to separate the valence electrons from the core electrons. The VASP calculation will be only applied on the valence electrons. First-principles pseu derived from an atomic reference state, requiring that the pseudo and all eigenstates have the same energies and amplitude (and thus density) outside a chosen core . When r is increased to r c , the plan wave is simplified by

The exact all electron wave function is Ψ(r) . Inside the radius r c , the exact wavefunction is replaced by a suitable "soft" pseudo wave function Ψ pseudo (r) (Figure 1).

potential of atoms

Many different basis sets have been proposed in the literature historically, and the two most Polynomials (Troullier and Martins (Ref N. Troullier 1993))

is used as an approximation for the simplified

The projected augmented wave (PAW) method [P. E. Blochl 1994, G. Kresse 1999], as Simulation Package (VASP) is used here.

In this approximation, core electrons are replaced by a pseudopotential V pseudo which is

M. Rappe1990)

The last one is the standard scheme for VASP pseudopotentials, which is also called "optimized". It is standard norm-conserving pseudopotential method.

All the PAW pseudopotentials used in the calculations are chosen in the library developed with the code (Table 1).An extension _h implies that the potential is harder than the standard potential and hence requires a greater energy cutoff. The extension _s means that the potential is softer than the standard version. For each chosen pseudopotential, all the possible combination sets will be calculated and compared. The results and conclusions will be presented in the next chapter 2 and 3. The principle characterization of the pseudopotentials and all the possible sets were listed in the

VII Periodic DFT code (VASP)

The Vienna Ab initio SimulationPackage, VASP, is a code for performing static and dynamics DFT simulations by using a periodic approach. By this approach, the system is defined by a supercell replicated infinitely in the three directions. It is well adapted to the solids and surfaces since it can simulate infinite systems and resolve the board problem. The process to do the calculations is as follow: • INCAR: It includes various keywords for the VASP calculations, the cut-of energy, the convergence parameters, spin-polarization parameters, the smearing parameters, density of states parameters… • POSCAR: This file contains the starting lattice geometry and ionic positions.

• POTCAR: List of pseudopotentials for each atom. The different parts for each element are in the same order as in the POSCAR file.

• KPOINTS: Number of k-points in the Brillouin zone.

Main analyzed output files:

• CONTCAR: It contains the lattice geometry and the ionic positions after optimization.

• OUTCAR: It records the information during the optimization… • OSZICAR: It includes the total energy of the supercell, and the total magnetic moments (if ISPIN=2).

• DOSCAR: List of density of states for the atoms of the supercell.

• PROCAR: It contains the spd and site projected wave function character of each band.

VIII Brillouin zone integration

In the VASP calculation, it is important to optimize the number of k-point in the Brillouin zone. So here in this part, the integration in the Brillouin zone will be presented. In solid state physics, the first Brillouin zone is a uniquely defined primitive cell in reciprocal space (Figure 4). The boundaries of this cell are given by planes related to points on the reciprocal lattice. It is found by the same method as for the Wigner-Seitz cell in the Bravais lattice. The importance of the Brillouin zone stems from the Bloch wave description of waves in a periodic medium, in which it is found that the solutions can be completely characterized by their behavior in a single Brillouin zone. The first Brillouin zone is the local points in reciprocal space that are closer to the origin of the reciprocal lattice than they are to any other reciprocal lattice points (see the derivation of the Wigner-Seitz cell). Another definition is as the set of points in k-space that can be reached from the origin without crossing any Bragg plane.

There are also second, third, etc., Brillouin zones, corresponding to a sequence of disjoint regions (all with the same volume) at increasing distances from the origin, but these are used less frequently. As a result, the first Brillouin zone is often called simply the Brillouin zone.

(In general, the n-th Brillouin zone consists of the set of points that can be reached from the origin by crossing exactly n -1 distinct Bragg planes.) In order to define the Brillouin zone we need to define first the reciprocal lattice. The set of all wave vectors k that yield plane waves with the periodicity of a given Bravais lattice is known as its reciprocal lattice. Analytically, lattice of points R, provided that the relation Eq. 11 holds for any lattice.

Factoring out, we can characterize the reciprocal lattice as the set of wave vectors satisfying Eq. 12for all R in the Bravais lattice.

There are also second, third, etc., Brillouin zones, corresponding to a sequence of disjoint regions (all with the same volume) at increasing distances from the origin, but these are used frequently. As a result, the first Brillouin zone is often called simply the Brillouin zone.

th Brillouin zone consists of the set of points that can be reached from the origin by crossing exactly n -1 distinct Bragg planes.) In order to define the Brillouin zone we need to define first the reciprocal lattice. The set of that yield plane waves with the periodicity of a given Bravais lattice is known as its reciprocal lattice. Analytically, k belongs to the reciprocal lattice of a Bravais , provided that the relation Eq. 11 holds for any r, and all

$\7]S^: = $\]
Factoring out, we can characterize the reciprocal lattice as the set of wave vectors in the Bravais lattice.

$\^= 1 30 There are also second, third, etc., Brillouin zones, corresponding to a sequence of disjoint regions (all with the same volume) at increasing distances from the origin, but these are used frequently. As a result, the first Brillouin zone is often called simply the Brillouin zone.

th Brillouin zone consists of the set of points that can be reached from the The reciprocal lattice is itself a Bravais lattice and its primitive vectors can be generated from the vectors of the direct lattice. Let `a, `b, `c be a set of primitive vectors, then the reciprocal lattice can be generated by the three primitive vectors:

d a = be `b × `c `a • 7`b × `c: d b = be `a × `c `a • 7`b × `c: d c = be `a×`b `a•7`b×`c: (13)
Using the relations between direct and reciprocal lattice, it can be shown that the reciprocal lattice of sc is sc (at k space), the reciprocal of bcc is fcc, and reciprocal of fcc is bcc.

The first Brillouin zone is defined to be the Wigner-Seitz primitive cell of the reciprocal lattice, or it could be defined as the set of points in k space that can be reached from the origin without crossing any Bragg plane. The second Brillouin zone is the set of points that can be reached from the first zone by crossing only one Bragg plane. The (n + 1) th Brillouin zone is the set of points not in the (n -1) th zone that can be reached from the n th zone by crossing (n -1) Bragg planes. Alternatively, the n th Brillouin zone can be defined as the set of points that can be reached from the origin by crossing (n -1) Bragg planes, but no fewer.

IX Density of States (DOS)

The total concentration of available states in certain energy range is provided by an interval of energy, which is

C = = +#+/hi 7 :C (14)
where:

• C is the interval of electron density presents in the energy range dE;

• = +#+/hi 7 :is the density of states.

So the total concentration of available states in the system between energy E1and E2 can be represented as

j Xj = A = +#+/hi 7 :C j j (15)
Obtaining = +#+/hi 7 : is accomplished through a density of states calculation. The resulting DOS is important since it can provide a better handle on the state distribution of the system.

Using N tot as the total number of electrons, the Fermi level is defined by:

BYB = A = +#+/hi 7 :C jZ X∞ (16)
where:

• E f is the Fermi energy.

For the VASP calculation, the DOS is integrated numerically in discrete intervals.

Different from the direct DOS plot, the derived integral DOS plot can present the localized states.

X Formation and binding energies of defects.

In the used fuel claddings, there are different point defects in the zirconium layer and oxide layer. So it is necessary to consider the effect of the point defects on the carbon diffusion.

Here in this part, two point defects are presented with their formation energies and the interaction by binding energies.

A point defect is formed when an atom is missing from a position, creating a vacancy, or when an atom occupies an interstitial site where no atom would ordinarily appear, causing an interstitial site. The point defects are showed in Figure 5. To understand the stabilities of the simple atom point defects, the formation energy of acarbon interstitial atom or a vacancy is calculated by Eq. 18 and Eq. 19:

m = nX@+oo -7 n + @+oo :

where:

• @+oo is the total energy of a supercell; • n is an atomic energy in a reference crystal (diamond structure); • nX@+oo is the total energy of a supercell with a carbon atom in an interstitial position.

If the formation energy of an interstitial C atom, m is positive, it means nX@+oo > n @+oo , the C atom is less stable in the interstitial position than in the reference crystal (diamond structure). lp@p#@i = @+oo7q$Br lp@p#@i ?:

? -@+oo (19)

where:

• @+oo7q$Br lp@p#@i ?: is the total energy of a supercell with a vacancy x in it; •

? is an atomic energy in a reference crystal; • @+oo is the total energy of a supercell.

If the formation energy lp@p#@i of a vacancy x is positive, it means @+oo7q$Br lp@p#@i ?:

? > @+oo , the atom x is more stable in the vacancy position than in the reference crystal.

When the inserted carbon is in a vacancy position, it is noted as substitution impurity. The formation energy is calculated as sWt = nX@+oo7q$Br sWtsB$BWB$Y#: -< n @+oo7q$Br lp@p#@i: >

where:

• nX@+oo7q$Br sWtsB$BWB$Y#: is the total energy of a supercell with C substitution ; • n is an atomic energy in a reference crystal (diamond structure); • @+oo7q$Br lp@p#@i ?: is the total energy of a supercell with a vacancy x in it.

If the formation energy sWt of a vacancy x is positive, it means nX@+oo7q$Br sWtsB$BWB$Y#: > n @+oo7q$Br lp@p#@i: , the C atom is more stable in the reference crystal than in the vacancy position.

However, the formation energy sWt is compared to m , if it is smaller than m , it can be said that for the C atom it is favorable to stay in the vacancy position than to be interstitial.

Sometimes there will be interstitial-vacancy interaction and substitution-vacancy interaction, at this time it better to use the binding energy to understand the interaction between the two point defects. The binding energy is calculated as below:

t$#k$#h = 7 n7$#B/sWt:X@+oo @+oo7q$Br lp@p#@i ?: : -7 n7$#B/sWt:X@+oo7q$Br lp@p#@i ?: @+oo :

(21)

n7$#B/sWt:X@+oo @+oo7q$Br lp@p#@i ?: n7$#B/sWt:X@+oo7q$Br lp@p#@i ?: @+oo

C V V C + - + ( ) (
)
where:

• n7$#B/sWt:X@+oo7q$Br lp@p#@i: is the total energy of a supercell with a C atom (interstitial or substitution) and a vacancy x; • n7$#B/sWt:X@+oo is the total energy of a supercell with a carbon atom in an interstitial position or vacancy position; • @+oo7q$Br lp@p#@i ?: is the total energy of a supercell with a vacancy x; • @+oo is the total energy of a pure supercell.

If t$#k$#h is positive, which means that n7$#B/sWt:X@+oo @+oo7q$Br lp@p#@i ?: > @+oo n7$#B/sWt:X@+oo7q$Br lp@p#@i ?: , the interaction between C atom and vacancy x is attractive, the C atom is more stable to be near vacancy than to be in the int/sub position far away from the vacancy. Other else, if the binding energy is negative, it means the interaction is repulsive, the C atom is more stable in the int/sub position.

XI Method to calculate the migration energy (NEB) and attempt frequency.

To study the carbon diffusion at atomic scale, it is necessary to characterize all the atomic jumps. For one carbon position, all the possible neighbor sites should be identified with reasonable distance and by using the Transition State Theory (TST) [E. P. Wigner 1937, H. Eyring 1935]. Once the potential landscape has been obtained for each path between neighbor sites, we can compute the transition rate between the initial and final sites 7u $Z : along this path.

It is a statistical method based on Boltzmann's distribution in phase space which leads to the migration rates through the energy barriers and attempt frequencies.

Diffusion by vacancy or interstitial mechanisms: the jump probability for diffusion is given by:

u = v $Z w 7- j x!y H * L : (22)
where:

• v $Z is the attempt frequency from the initial site to final site;

• "$h is the migration energy; • N . is the constant Boltzmann;

• is the simulation temperature.

v $Z and "$h depend on the local environment of the migrating atom. Thus these 2 quantities have to be determined for all the distinct atomic jumps.

"$h is calculated thanks to the NEB method and v $Z is calculated by the vibration frequencies of carbon atom at stable site and transition state.

XI.1 Migration energy calculation (NEB)

Using the precedent methodologies, the studied model can be constructed and optimized to find its stable state. To study the behavior of an element (C in our project), first, the possible interstitial sites will be investigated to study its atomic migrations. Then, by using Nudged Elastic Band (NEB) method (Ref H. Jonsson 1998), the Minimum Energy Path (MEP) between these sites will be determined (Figure 7). To use the NEB method, the initial state and final state should be stable. So firstly, the initial site and final site are optimized by VASP, and their equilibrium stable positions are used as the limit states. Then, several "images" are built between these two states. These images are specific C positions, "guessed" by locating the diffusing atom in a "thought to be reasonable" diffusion way. After that, each image finds its ionic position which has the lowest energy (or force, depends on the command used), while the images maintain equal spacing with each other. This equal-spacing-maintaining ability is realized by the NEB, who adds spring forces along the band, and projects out the component of the force due to the potential perpendicular to the band.

There are a few improvements to the NEB method which are not yet included in the current version of VASP. A climbing image method and a better tangent definition combine to allow for the more accurate finding of saddle points using the NEB with fewer images than the original method. The setup and operation of this implementation can be identical to what is described in the VASP manual under the elastic band section. The climbing image method can be turned by setting LCLIMB = .TRUE. in the INCAR file [G. Henkelman, 2000]. With Climbing Image NEB method, it is possible to calculate not only the MEP for the reaction but also the transition state configuration at the saddle point. The climbing image NEB method has one modification which makes the image with the highest energy up to the saddle point (Figure 8). This image does not see the spring forces along the band. Instead, the true force at this image along the tangent is inverted. In this way, the image tries to maximize its energy along the band, and minimize in all other directions. When this image converges, it will be at the exact saddle point.

Because the highest image is moved to the saddle point and it does not feel the spring forces, the spacing of images on either side of this image will be different. It can be important to do some minimization with the regular NEB method before this flag is turned on, both to have a good estimate of the reaction coordinate around the saddle point, and so that the highest image is close to the saddle point. If the image is initially very far from the saddle point, and the climbing image was used from the outset, the path could develop very different spacing on either side of the saddle point. The energy barrier for diffusion atom to cross during the migration can be calculated to be the energy difference of the initial site and transition site, using

"$h = B/p#s$B$Y# -$#$B$po (23)
where:

• "$h is the migration energy for atomic jump from the initial site to final site;

• B/p#s$B$Y#
is the total energy of the system when the diffusion atom is at its transition site; • $#$B$po is the total energy of the system when the diffusion atom is at its initial site.

In our calculation 7 or 8 images are considered for all migration barrier calculations. The simulations of CINEB are performed with k-points grid 2×2×2 on the supercell of 96 atoms.

XI.2 Attempt frequency calculation

When the atom is at the transition state, it is in the equilibrium position. And its potential is harmonic. For one atom, it has three proper modes. So for the N atoms in the bulk, there are 3N frequencies.

For the atom at the transition state, since it is not stable, it has an imaginary frequency, and all the other two frequencies are real. For the model used in the thesis, it contains an atom which does the diffusion in the supercell. Except the diffusion atom, all the other atoms are almost stayed at the same positions, so the change of their vibration frequencies at initial site and transition site can be neglected and only the attempt frequency of carbon atom will be considered. By the Vineyard theory [G. H. Vineyard 1957], the attempt frequency can be calculated using

v $Z = ∑ z ! !{!|!}• € ∑ z ! |•}{'!|!ƒ{ •€"… = z … z z • z … † z † (24)
where:

• v $Z is the attempt frequency;

• v $ are the C vibration frequencies at the initial state;

• v $ ‡ are the C real vibration frequencies at the transition state. Now, the attempt frequencies for the possible atomic jumps are investigated, so the jump probabilities of carbon atom can be calculated with Eq. 20. The characterizations of carbon migrations at atomic scale have been performed. In the next step, the Kinetic Monte Carlo method can be used in the macroscopic scale simulation to calculate the diffusion coefficient of carbon in bulk.

XII Kinetic Monte Carlo method (KMC)

Once calculating the atomic jumps with NEB method, the carbon diffusion simulation at the macroscopic scale can be performed with the kinetic Monte Carlo (KMC) method. It is a Monte Carlo method intended to simulate the time evolution of some processes occurring in nature. Typically, these are processes that occur with a given known rate. It is important to understand that these rates are inputs to the KMC algorithm, the method itself cannot predict them.

The main idea behind KMC is to use transition rates that depend on the energy barrier between the states, with time increments formulated so that they relate to the microscopic kinetics of the system.

For example,diffusion in the bulk with the activation energy in macroscopic scale is represented at atomic scale with the migration energy of the carbon atom which is energy barrier from one site to another neighbor one (Figure 9) and the rate constant for diffusion can be calculated using VTST theory, which is presented as jump probability in Part XI.

Before starting KMC simulation we have to make a list of all possible events that can be realized during the simulation and calculate rates for each event. When the rate constants of all processes are known, we can perform KMC simulation in the time domain. In the case of a single process, the reciprocal of the rate of the process determines the time required for the reaction occurs. This quantity can be set equal KMC time. In the case of a many time is less straightforward and several modifications of KMC exist.

First, the rates of all atomic jumps from the initial site to final states have to be (Figure 10). the rate constants of all processes are known, we can perform KMC simulation in the time domain. In the case of a single process, the reciprocal of the rate of the process determines the time required for the reaction occurs. This quantity can be set equal KMC time. In the case of a many-particle multi-process system, however, introduction of time is less straightforward and several modifications of KMC exist.

rates of all atomic jumps from the initial site to final states have to be the rate constants of all processes are known, we can perform KMC simulation in the time domain. In the case of a single process, the reciprocal of the rate of the process determines the time required for the reaction occurs. This quantity can be set equal to the process system, however, introduction of rates of all atomic jumps from the initial site to final states have to be evaluated For a single constant time step, the probability of the considered processes is between 0 and 1 (Figure 11).When the rate constants of all processes are known, we can perform KMC simulation in the time domain. Time increments are defined by the rates of all processes and are formulated so that they relate to the microscopic kinetics of the system. By time residence algorithm [Young 1966], the average residence time dt equals to the inverse of the sum of the probabilities.

Cˆ= ∑ ‰ ! { ! (25)
This method should be used when kinetics rather than equilibrium thermodynamics dominates the structural and/or compositional changes in the system.

XIII Diffusion coefficient calculation

In the previous part, it is known that the carbon does random jumps in the bulk, the displacement of the atom is noted as r. For the diffusion in three dimensions, there is relation between the mean-square of displacement (MSD) and the diffusion coefficient as below:

< > = 6OEˆ (26)
In Figure 12, the slope of the curve equals to the 6D. So the diffusion coefficient at temperature T could be calculated. The diffusion coefficient in solids at different temperatures is generally found to be well predicted by the Arrhenius equation:

OE = OE exp 7- j }•| H * L : (27)
where:

• OE is the diffusion coefficient;

• OE is the pre-factor for calculation of D;

• is the absolute temperature;

• p@B is the activation energy in the diffusion; • N . is the boltzman constant.

This equation can be transformed to:

' OE = ' OE - j }•| H * L (28)
The equation presents time dependence of diffusion coefficient. The linear relation between lnD and temperature T is presented in Figure 13. So by simulating the carbon's trajectories at different temperatures, the related diffusion coefficients at the temperatures can be calculated. By doing the Arrhenius plot for T and D, the prefactor D 0 and the activation energy E act can be calculated.

Conclusion

The carbon atom diffusion behavior is performed with multi-scale method simulations in this thesis. At atomic scale, the calculation model is built from the unit cell. The geometrical parameters and calculation parameters such as cut-off energy and pseudopotential set will be optimized by DFT calculations using VASP code. All the calculations are performed with the GGA functional PW91 [Perdew 1992]. Then to study the carbon atomic jumps, all the possible sites will be investigated in a large supercell (to avoid the volume variation caused by carbon insertion). For each two neighbor interstitial sites, the migration energy and attempt frequency will be calculated using VTST theory and the jump probabilities for all possible atomic jumps will be investigated. At macroscopic scale, using the atomic jump results of carbon from the atomic scale simulations, the carbon's trajectory at certain temperature will be performed with KMC approach. By repeating the simulation of carbon at different temperatures, the relation between diffusion coefficient and temperature can be calculated.

Finally diffusion coefficient of carbon at storage condition can be calculated.

All the multi-scale simulations will be performed in hcp zirconium and monoclinic zirconia.

The diffusion mechanism in this work contains two types of mechanism: diffusion by interstitial sites, and diffusion by vacancy sites.

Chapter 2

Carbon diffusion in bulk zirconium

I.1.B. Optimization of the

I.2 Macroscopic

Introduction

In the PWR reactor, the cladding is made of zircaloy-4, which is principally zirconium (98.4%) and contains oxygen atoms, which may become C14 under irradiation. So to study the carbon behavior and diffusion in the claddings, it is necessary to perform the carbon simulation in the zirconium.

The hexagonal close packed (hcp) structure is the crystallographic structure of zirconium alloys. So in this study, the hcp phase was chosen for the simulation in zirconium bulk. The simulations are performed to characterize, by the mean of a multi-scale approach, the carbon diffusion behavior in the Zr hexagonal close packed bulk to estimate the waste safety of the fuel rods' claddings in the storage condition.

Very few experimental studies of C in Zr have been published. In the G. Murch's experiments, the measurements at constant hydrostatic pressure rather than constant volume were made as implied in the lattice gas treatment. It was found that the compositional dependence of the tracer carbon diffusion coefficient was reasonably well represented using the partial molar energy of carbon as input to the equation for the diffusion coefficient [G. E. Murch 1979]. So it means that in theoretical study, an isolated carbon atom with its molar energy could be chosen to perform the carbon behavior which can describe the carbon diffusion coefficient very well.

For the C behavior in the zirconium bulk, there was only one old experimental study to measure the diffusion coefficients of C. The diffusion of carbon in zirconium was studied in the temperature range 873-1523 K, using the residual activity technique. [R. P. [START_REF] Agarwala | [END_REF]. It showed that the diffusion coefficient at 923 K was 5.1×10 -15 cm 2 /s.

However, regarding theoretical calculations, to our knowledge, no calculation of C behavior in bulk hcp zirconium have been performed.

Few simulation studies by first principle calculations have been performed with Zr. For example, the self-diffusion was performed by the first principal calculation in the reference in which the average value of the migration energy was found to be in a very good agreement with available experimental data. 2009] and experiments [M. Iwasawa 2008]. For the diffusion in the grain boundaries, it was proved by calculations that both the formation energies and the migration free energies were significantly lower than in the bulk which demonstrates the role of the grain boundary as a vacancy sink and a fast diffusion channel [J. R. Fernandez 2000].

In this chapter, the carbon diffusion coefficient in the zirconium has been investigated using a multi-scale approach: atomic scale using the periodic density functional theory (DFT) and macroscopic scale with Kinetic Monte Carlo (KMC) simulations.

In the first part, a Zr bulk model was optimized and the possible interstitial carbon sites were investigated. Then, for each interstitial carbon site, all the neighbor interstitial sites were identified so that all the migration possibilities between each two neighbor sites were characterized. For the possible jumps, the migration energies, the attempt frequencies and the migration rates were calculated.

In the second part, all the results of the previous calculations were used to simulate the carbon diffusion by using the Kinetic Monte Carlo (KMC) method. Once knowing the carbon diffusion behavior in the perfect hcp Zr, the vacancy effect was considered to calculate the carbon diffusion in a more realistic environment.

I: Carbon behavior inzirconium

In this work, the carbon diffusion simulation will be performed using multi-scale simulations. That means the simulation has two scales: atomic scale to study the atomic jump between two neighbor sites and macroscopic scale to study the macroscopic diffusion in bulk.

I.1 Atomic scale simulations

As the first step, in atomic scale simulations, calculation parameters were optimized by comparing calculated properties to experimental ones for zirconium bulk.

I.1.A Optimization of the calculation parameters

A-1 Zr model

The zirconium metal under ambient environment is a hexagonal close-packed (hcp) The crystal structure of Zr can be constructed by 2 atoms in one unit cell (Figure 2), a motif with coordinates (0, 0, 0) and ( pseudopotential has to be chosen to do the simulation.

The zirconium electron configuration is pseudopotentials are considered for Zr: Zr_4 and Zr_12. For the four electrons (4d 2 5s 2 ) are considered as valence electrons. For the Zr_12 pseudopotential, the last twelve electrons (4s 2 4p 6 4d was evaluated by comparing the tot the experimental values when the cut

The geometrical optimization was performed with different k energies. For Zr_12, the total energy conv k-points grid. But, for Zr_4, the total energy had already converged at 200 eV with a 11×11×11 k-points grid. With (Zr_4) and (Zr_12), the calculated bulk parameters of hexagonal zirconium are the same: a= 3 the experiment ones 3.23Å and 5.15 Å [C. Kittel 2005] electrons and has a higher convergence energy which needs more calculation time, (Zr_4) is thought to be more suitable to be used in the later calculations.

numbered 194[B. Olinger 1973], with parameters: a=3.232 Å, b=3.232 Å, c=5.147 Å, α=90°, The crystal structure of Zr can be constructed by 2 atoms in one unit cell (Figure 2), a motif with coordinates (0, 0, 0) and ( G , G , ). Once the initial structure was built, the corresponding pseudopotential has to be chosen to do the simulation.

The zirconium electron configuration is:1s 2 2s 2 2p 6 3s 2 3p 6 3d 10 4s 2 4p 6 4d 2 5s pseudopotentials are considered for Zr: Zr_4 and Zr_12. For the Zr_4 pseudopotential, the last ) are considered as valence electrons. For the Zr_12 pseudopotential, the 4d 2 5s 2 ) are considered as valence electrons. The pseudopotential was evaluated by comparing the total energy convergence and geometrical parameters with the experimental values when the cut-off energies and k-points change.

The geometrical optimization was performed with different k-points and different cut energies. For Zr_12, the total energy converges from 400 eV cut-off energy with a 11×11×11 points grid. But, for Zr_4, the total energy had already converged at 200 eV with a points grid. With (Zr_4) and (Zr_12), the calculated bulk parameters of hexagonal zirconium are the same: a= 3.23 Å, c= 5.17 Å which have a good agreement with the experiment ones 3.23Å and 5.15 Å [C. Kittel 2005]. Since Zr_12 contains more valence electrons and has a higher convergence energy which needs more calculation time, ( Zr_4) is e to be used in the later calculations. points grid. But, for Zr_4, the total energy had already converged at 200 eV with a points grid. With (Zr_4) and (Zr_12), the calculated bulk parameters of .23 Å, c= 5.17 Å which have a good agreement with . Since Zr_12 contains more valence electrons and has a higher convergence energy which needs more calculation time, (Zr_4) is

In order to compare the optimization result of bulk Zr with previous theoretical and experimental results, the cohesive energy for a Zr unit cell was also calculated, using the following formula:

@Yr+s$l+ = 7 BYBpo -pBY" * pBY" :/ pBY"

where:

• BYBpo : total energy of the Zr unit cell; • pBY" : energy of an isolated Zr atom; • pBY" : number of atoms in the Zr unit cell.

Since in the Eq (1), it is necessary to know the energy of a single Zr atom, a cell of 10 Å×10 Å×10 Å with a single Zr atom in its center was built to calculate it. Then the cohesive energy for a Zr unit-cell can be calculated. The Zr_4 pseudopotential gives accurate results comparing with the experimental and theoretical values. So Zr_4 can be used for this system. Now, in order to study the carbon diffusion, the carbon's pseudopotential has to be tested.

First, a pure diamond structure is considered to perform the pseudopotential evaluations.

A-2 Diamond C structure

The pure diamond under ambient environment is a diamond-cubic structure, atomic number 6, and atomic mass 12.011 a.m.u. Its space group is Fd3m, numbered 227 [K. Kobashi 2005] with parameters: a=3.567 Å, b=3.567 Å, c=3.567 Å, α= 90°, β= 90°, γ= 90°. In our simulation, the crystal structure of diamond was built with 8 atoms in one unit cell, a motif with coordinates (0, 0, 0) (

( G , G , ).
The electron configuration of carbon is: noted "C", "C_s" and "C_h" are tested in this

As indicated in their name, "C_s" is a "soft" pseudopotential, "C_h" is a "hard" pseudopotential, and "C" is a "normal" one. A "hard" pseudopotential means it has strongly oscillating pseudo-wave functions, By comparing the convergence of associated total energy and geometrical parameters depending on the cut-off energies and k found.

For "C", the total energy converges from k (550 eV). The geometrical parameters converge from 3×3×3 k 450 eV.

For C_s, the total energy converges from k Geometrical parameters converge from 3×3×3 k For C_h, the total energy and geometrical parameters converge from 3×3×3 K very slowly with the cut-off energy.

cell of C diamond structure

In our simulation, the crystal structure of diamond was built with 8 atoms in one unit cell, a with coordinates (0, 0, 0) ( , , ), (0, , ), ( G , , G ), ( , 0, ), ( ,

The electron configuration of carbon is:1s 2 2s 2 2p 2 . Three types of pseudopotentials of carbon noted "C", "C_s" and "C_h" are tested in this part. They all have 4 valence electrons (

As indicated in their name, "C_s" is a "soft" pseudopotential, "C_h" is a "hard" pseudopotential, and "C" is a "normal" one. A "hard" pseudopotential means it has strongly wave functions, and requires many pseudo-wave functions in calculations.

By comparing the convergence of associated total energy and geometrical parameters off energies and k-points mesh, the more suitable pseudopotential is energy converges from k-point 5×5×5, but slowly with cut (550 eV). The geometrical parameters converge from 3×3×3 k-point and a cut For C_s, the total energy converges from k-point 5×5×5 with a cut-off energy of 550 eV.

eometrical parameters converge from 3×3×3 k-point and a cut-off energy of 350 eV.

For C_h, the total energy and geometrical parameters converge from 3×3×3 K off energy.

53

In our simulation, the crystal structure of diamond was built with 8 atoms in one unit cell, a G , G ), ( , , 0), and

. Three types of pseudopotentials of carbon part. They all have 4 valence electrons (2s 2 2p 2 ).

As indicated in their name, "C_s" is a "soft" pseudopotential, "C_h" is a "hard" pseudopotential, and "C" is a "normal" one. A "hard" pseudopotential means it has strongly wave functions in calculations.

By comparing the convergence of associated total energy and geometrical parameters points mesh, the more suitable pseudopotential is point 5×5×5, but slowly with cut-off energy point and a cut-off energy of off energy of 550 eV.

off energy of 350 eV.

For C_h, the total energy and geometrical parameters converge from 3×3×3 K-point, but So C_h pseudopotential is not suitable for the future calculations because of its slow convergence of total energy and geometrical parameters. C and C_s are both good candidatures for future calculations.

Also, to compare our optimization result with the previous calculation results and experiment ones, the cohesion energy will be calculated. A cell of 10Å×10Å×10Å with a single C atom in its center was built to get the energy of one carbon atom pBY" .

So the cohesive energy could be calculated using Eq (1). The results were listed in In our simulation, we construct the crystal structure of ZrC using 4 Zr atoms and 4 C atoms in one unit cell, a motif with coordinates (0, 0, 0) (0, , ), ( , 0, ), ( , , 0) for Zr and ( , 0,0), (0,0, ), ( , , ), (0, , 0) for C was used. Two combinations of pseudopotentials are tested in this part. We study one ZrC structure using the first set of pseudopotentials (Zr_4/C), and the other one using the second set of pseudopotentials (Zr_4/C_s). By comparing their total energy convergence and geometrical parameters as a function of the cut-off energies and k-points, we can determine which set of pseudopotentials is more suitable for this system.

For the set (Zr_4/C), the total energy converges from 9×9×9 K-points with a cut of 500 eV, and the geometrical parameters converge with a cut-off energy of 450 eV and k

For the set (Zr_4/C_s), the total energy converges from 9×9×9 K-points with a cut The cohesive energy is also compared with experimental and theoretical values. Results

were listed in Total energies and geometrical parameters for the (Zr_4/C) set:

Total energies and geometrical parameters for the (Zr_4/C_s) set:

A-4 Summary

Based on the previous results, the pseudopotential set simulations.

I.1.B. Optimization of the Zr bulk super

In the previous parts of this study, pseudopotentials, cut optimized to simulate pure Zr and C structures as well as ZrC one. However, to simulate the behavior of one carbon atom in the Zr he optimize the size of a Zr model which should be big enough to correctly describe the real system when a carbon is inserted (which means a weak geometrical relaxation of the cell volume in order to perform later concentration).

In this part, the geometrical parameters are calculated for different supercell sizes of Zr before and after the insertion of one carbon atom. By comparing the cell volume differ (before and after insertion), an accurate supercell model can be optimized to study C behavior in the Zr hexagonal close-packed structure.

By reproducing the unit cell of zirconium (Figure 6) in three dimensions, different sizes of supercells were obtained. For example, if the unit cell is reproduced 2 times along each axe (a, b and c), it will be a 2×2×2 supercell. At the first step, the variation sizes of supercells are calculated, when one carbon atom is added in the structure in site 8 (Table 4). The following supercells were tested (Figure 7). Taking into account: acceptable calculation time; small variation of geometric parameters; convenience for future calculation of vacancy defect, the 4×4×3 supercell has been chosen for next calculations (Figure 8). At the second step, calculation will be performed at constant volume, the cut-off energies are optimized, using the K-points large enough. We find that the cut-off energy of 300 eV can be used in the future calculation.

E cut =450 eV n(atom) Size ∆(a) (Å) ∆(a) (%) ∆(c) (Å) ∆(c) (%) ∆(V) (Å 3 ) ∆(V) (
At the third step, we optimize the k-points grid, and find the k-points 3×3×3 can be used in the future calculation. Now in the next part, the investigation of possible carbon interstitial sites can be calculated.

I.1.C Investigation of the carbon interstitial sites

The supercell model used in this part was built from the 2 atomic unit cell extended by 4×4×3 in the three directions as presented in the last part. By the optimization results, k-points mesh was optimized to 3×3×3 with a cut-off energy of 300 eV and constant volume calculations will be performed with this model.For the frequency calculation the convergence criteria of the force on atoms is set to be 10 -5 eV/Å in order to get an accurate relaxed configuration which is required to obtain an imaginary frequency at the saddle point.

To calculate the carbon diffusion in the Zr bulk, at first, the possible C interstitial sites were identified as well as their relative stabilities. It showed that there are two possible interstitial sites in the Zr hexagonal compact structure: an octahedral site named as O site (Figure 9 In fact, from an initial BT position (Fig. 9 panel a) after the relaxation process, it was found that there are three stable symmetrical BT_d sites as a triangle (Figure 9 panel c). The distances between a starting BT position and BT_d sites are very small at about 0.4 Å.

The formation energy of a carbon in interstitial position (E I ) was calculated using Eq (2): m = nX)/ -7 n + )/ :

(2)

where:

• )/
is the total energy of a 96 atomic zirconium supercell; • n is a carbon atomic energy in a diamond structure; • nX)/ is the total energy of a 96 atomic zirconium supercell with a carbon atom in an interstitial position.

Calculations showed that the carbon in the O site has a formation energy of -1.70 eV for Zr_4 or -1.68 eV for Zr_12 and in the BT_d site the formation energy is -0.25 eV for Zr_4 and -0.14 eV for Zr_12. It means that the O site is the more favorable one (∆E=E I (BT_d)-E I (O)=1.45 eV for Zr_4 and 1.54 eV for Zr_12). Thus, most of the carbon would stay in the O site (cf Table 5).

Configuration (eV) Zr_4 Zr_12

O -1.70 -1.68 BT_d -0.25 -0.14 By comparing the formation energy, it was found that for the two pseudopotentials (Zr_4 and Zr_12), the formation energies are very similar for O site, but there is a small difference for the distorted BT site. Since a carbon can jump between different interstitial sites, the migration energies for all possible jumps of Zr_4 and Zr_12 will be calculated in the migration part. However, the calculations of diffusion coefficient were only performed with Zr_4, the results of Zr_12 are given for comparison.

I.1.D Migration between neighbor sites

To study the carbon's diffusion in zirconium bulk, it is first necessary to calculate the needed migration energy, E mig , to jump from one interstitial site to a neighbor one.

D-1 Migration energies (NEB)

The Nudged Elastic Band method (NEB) [H. Jonsson 1998] was used to calculate the minimum energy pathway between the possible jumps in order to identify transition states and migration energies (Chapter I.X).

Following each migration pathway between two neighbor sites, several (7 or 8) carbon positions (images) were considered. After optimization, the energetic profile was plotted according to the total energy of these different points (Figure 10).

First, we test the jump between two neighbor BT_d sites in the zirconium bulk and it was found that the migration path would always pass through the BT position with an energy barrier 0.40 eV which corresponds to the energy difference of the two types of site (Figure 10 panel a).

Then, the migration between BT_d site and O site was calculated. As it can be seen in

u = v $Z w 7- j x!y H * L : (3)
Since all the migration energies for the possible atomic jumps have been calculated, for certain temperature T, the corresponding attempt frequency v $Z should be calculated to get the jump probability.

The results of the calculated attempt frequencies corresponding to the possible migration were listed in the As the first conclusion, at low temperature, the interstitial carbon should preferentially diffuse through the path BT_d site and O site, which are connected. For all possible migration paths, the migration energies, the attempt frequencies and their jump probabilities were calculated in this part. The results were summarized in Table 8.

Possible jump

E mig νif Γif O site O site 3.51 eV 1.3×10 13 s -1 1.3 × 10 G š X w (- 3.51 N . )
O site BT_d site 2.13 eV 1.7×10 13 s -1 1.7 × 10 G š X w (-2.13 N . )

BT_d site O site 0.70 eV 1.6×10 13 s -1 1.6 × 10 G š X w (-0.70 N . )

BT_d site BT_d site 0.40 eV 1.3×10 13 s -1 1.3 × 10 G š X w (-0.40 N . )

Table 8: The results of all possible atomic jumps between each two neighbor sites

I.2 Macroscopic scale simulations

After doing the elementary mechanism calculations at the atomic scale, the macroscopic scale simulations were performed by using Kinetic Monte-Carlo simulations (KMC) [C. Kittel 2005, H. Jonsson 1998] to study the carbon diffusion in the zirconium bulk. The main idea behind KMC is to use migration rates that depend on the energy barrier between the states and the time step is given by the residence time algorithm [C. S. Becquart 2010]

I.2.A KMC simulations

The isolated carbon used in the KMC simulation does a random walk. From the carbon's trajectory of three dimensions at a given temperature (Figure 11), the mean square displacement of the carbon could be calculated. The theoretical relation between r

<

where:

• is the displacement of carbon;

• OE is the C diffusion coefficient;

• ˆis the diffusion time.

The diffusion coefficient, at a fix temperature, can thus be obtained from the slope of Eq. (4).

z(Å) lated carbon used in the KMC simulation does a random walk. From the carbon's trajectory of three dimensions at a given temperature (Figure 11), the mean square displacement of the carbon could be calculated.

Example of C's atomic trajectory at T=900K for 100000steps

Its linear relation with the residence time could be proved by drawing the graphic of as a function of time t.

The theoretical relation between r 2 and t could be described with the formula:

>= 6OEˆ
is the displacement of carbon;

is the C diffusion coefficient;

The diffusion coefficient, at a fix temperature, can thus be obtained from the slope of Eq. (4).

y(Å)

x(Å) lated carbon used in the KMC simulation does a random walk. From the carbon's trajectory of three dimensions at a given temperature (Figure 11), the mean square Its linear relation with the residence time could be proved by drawing the graphic of and t could be described with the formula:

(4)

The diffusion coefficient, at a fix temperature, can thus be obtained from the slope of Eq. (4).

Take an example at 900 K here in Figure 12, MSD is mean square displacement< >, the slope of the curve is 6D. So by the figure, the diffusion coefficient at 900 K is calculated to be 2.93×10 -14 cm 2 /s. After doing the same simulations for different temperatures from 700 K to 1100 K, the diffusion at different temperatures were calculated and summarized in Table 8. The relation between diffusion coefficient and the effective activation energy can be obtained using Eq. 5:

OE = OE exp 7- j }•| H * L : (5)
where:

• OE is the diffusion coefficient;

• OE is the prefactor of D;

• is the temperature; • p@B is the effective activation energy.

Indeed, Eq. 5 can be changed as Eq. 6 to find the linear relation between ' OE and 1/T.

' OE = ' OE - j }•| H * L (6)

I.2.B Determination of the Diffusion coefficient

By applying the formula on the calculation results in different temperatures (Figure 13), the prefactor of diffusion coefficient D 0 was calculated as 6.02x10 -3 cm 2 /s and the effective activation energy of carbon's diffusion turned out to be 2.02 eV, leading to: OE = 6.02 × 10 XG (žŸ /s)exp 7- The storage temperature condition for the used fuel claddings is at around 50°C, so the diffusion coefficient of a carbon by interstitial mechanism in perfect Zr crystal at this temperature is calculated to be ~10 -34 cm 2 /s.

I.3 Conclusion

The calculated diffusion coefficient is so small that we could conclude the carbon in perfect Zr bulk would be blocked in it since it is very hard for carbon to diffuse and so to exit from an ideal zirconium bulk structure.

Comparing our result with the experimental value, for C diffusion in Zr, the experimental diffusion coefficient at 923 K is about 5.1×10 -15 cm 2 /s [R. P. [START_REF] Agarwala | [END_REF] and in our calculation at 923 K the carbon diffusion coefficient is 5.6×10 -14 cm 2 /s which has a good agreement.

II: Carbon behavior in zirconium bulk with point defects

So far, the carbon behavior in the zirconium bulk had been studied. In the fuel claddings, there is not perfect zirconium but point defects are present especially under and after irradiation. So it is important to study the effect of vacancy on C behavior in the bulk. All the calculations of defects were performed on the 96 atomic zirconium supercell.

II.1. Zr vacancy in pure Zr bulk

The formation energy to form an atomic vacancy was calculated using the Eq. 8:

In this last part, the Zr vacancy effect on the formation energy of an interstitial C has been lp@p#@i: > (10)

is the total energy of 96 atomic supercell with a nearest Zr vacancy is the total energy of the 95 atomic zirconium supercell with a Zr The result of formation energies turned out to be -0.85 eV for C in O site and 0.46 eV for C in BT_d site when there is a nearby vacancy (Figure 14).

It showed that the carbon is relatively stable at the interstitial O site even when there is a nearby vacancy, comparing with the substitution formation energy 1.82 eV. So it was more stable for C to be in the interstitial position than in the substitution. However, considering the -1.70 eV formation energy for C in O site for an interstitial carbon without a vacancy, it is found that the carbon is getting less stable when it is near the vacancy.

II.4 Vacancy effect on carbon diffusion by calculation of binding energy

This conclusion can also be found by calculating the binding energy. The binding energy is the interaction energy between the vacancy and the interstitial carbon. It can be calculated using Eq. 11:

t$#k$#h = 7 @X)/ + @+oo7q$Br lp@p#@i: : -7 nX)/7q$Br lp@p#@i: @+oo :

where:

• nX)/7q$Br lp@p#@i: is the total energy of 96 atomic supercell with a nearest Zr vacancy from an interstitial carbon; • nX)/ is the total energy of a 96 atomic zirconium supercell with a carbon atom in an interstitial position; • @+oo7q$Br lp@p#@i: is the total energy of the 95 atomic zirconium supercell with a Zr vacancy; • @+oo is the total energy of the 96 atomic zirconium supercell.

The vacancy-carbon binding energy is -0.86 eV considering C in O as reference (Table 9).

The negative value means that it is less stable when the carbon interstitial has a vacancy nearby than the carbon interstitial is far away from the vacancy. Therefore, it can be said that the closer an inserted carbon is to a Zr vacancy, the less stable it is. So the carbon always prefers to stay away from zirconium vacancies. In the hand of diffusion, we can say that the zirconium vacancy wouldn't decrease the carbon diffusion.

Conclusion

We have calculated the carbon diffusion coefficient in pure Zr bulk to model the carbon behavior in used fuel claddings by coupling atomic and macroscopic scale simulations. For the DFT atomic scale simulation, a 96 atomic hcp zirconium supercell was considered. Two types of stable interstitial sites for carbon were found: BT_d site which is less stable than the O site by 1.45 eV for Zr_4. The migration energies calculated from the DFT-NEB approach showed that the easiest migration pathway for diffusion should be the successive diffusion between BT_d site and O site. Then, the attempt frequencies for each possible jump were calculated in order to get, from the transition state theory, the migration rates between each two neighbor sites.

In a second step, a macroscopic simulation by the Kinetic Monte-Carlo approach (KMC) was performed using jump probabilities. From the KMC simulations at different temperatures, the effective activation energy as well as the diffusion coefficient as a function on temperature Zr crystal, at the maximum storage temperature (50°C), is calculated to be very small (~10 -34 cm 2 /s). Therefore, in this condition, it will be very difficult for a carbon atom to diffuse in a perfect Zr crystal. However, even though there arevacancies in Zr, it was calculated that the C stability decrease when it gets closer to a Zr vacancy. Therefore, the zirconium with vacancies wouldn't increase the carbon diffusion by the interstitial diffusion mechanism. Carbon would be blocked in the bulk zirconium.

Chapter 3

Carbon Diffusion in the bulk monoclinic zirconia

Introduction

In the previous chapter, the carbon diffusion coefficient in hcp zirconium has been calculated. Since the external part of zirconium cladding contacts with water, it will form an oxide layer composed of monoclinic zirconia and tetragonal zirconia. At the storage temperature, monoclinic zirconia is the most stable phase in the oxide layer, so this phase is chosen to perform the carbon diffusion behavior. As the same as zirconium, the simulations are performed in bulk in this chapter.

On the one hand, some experiments have proved that the monoclinic structure is the most stable one for zirconia [A.S. [START_REF] Foster | [END_REF].

On the other hand, some DFT calculations have been performed for zirconia, either monoclinic structure or tetragonal structure[G. Stapper 1999, G. M. Hood 1993, H. Jiang 2010]. The static simulations were used to confirm the good geometrical parametersand structures of zirconia. The theoretical results are in a good agreement with some experiments, such as the parameters of zirconia. However, theyarenotin agreement with other experimental results, such as the activation energy. In this work, the theoretical results will be compared to the experimental results and other theoretical ones.

In this chapter, the carbon diffusion calculations were performed with multi-scale approach in bulk zirconia. For the pure bulk, the diffusion by interstitial mechanism will be studied.

Different interstitial sites for carbon will be investigated and the atomic jumps between interstitial sites will be calculated through the migration energies and transition rates. The macroscopic scale simulations using KMC method will be performed to evaluate the diffusion coefficient at storage temperature. However, there are vacancies in the oxide layer after irradiation. So, the diffusion by vacancy mechanism will also beinvestigated. First, the selfdiffusion of O will be calculated. Calculations will also be performed to investigate the possibilities of carbon diffusion by vacancy mechanism.

I Carbon behavior in perfect zirconia

I.1 Atomic scale simulations

The study of carbon behavior in monoclinic zirconia is performed as the one in bulk zirconium, by using the multi-scale simulations: the atomic jumps with DFT simulations at the atomic scale, and the diffusion with KMC approach at the macroscopic scale. So at first, a unit cell of monoclinic zirconia will bebuilt to study carbon diffusion.

I.1.A Optimization of zirconia calculation parameters

To do the simulation of carbon behaviour in zirconia, a 12 atomic cell will be built which is the smallest cell to model the monoclinic zirconia. To decide which pseudopotential set can better describe the structure of monoclinic zirconia, different pseudopotentials were used for Zr, O and finally ZrO 2 .

The optimization of bulk zirconium (which is also hexagonal close packed phase here) was The cell parameters still do (Zr_12/O_s), the total energy and the geometrical parameters of the = 350 eV. There is no significant difference of the total . Here are the figures udopotential set Zr_4/O_s (Figure 3 and ). But Zr_12 which contains more valence electrons need more Since the set of pseudopotentials Zr chosen to optimize the ZrO 2 cubic and tetragonal structure.

I.1.B Study of relative stability of the three phases of ZrO

ZrO 2 has 3 types of structures: monoclinic, tetragonal, cubic. Now, the optimizations of parameters for structure cubic and tetragonal are also performed to study their stabilities. Since the set of pseudopotentials Zr_4 and O_s seems to be a goodchoice in this part, cubic and tetragonal structure.

stability of the three phases of ZrO 2 (tetragonal, and has 3 types of structures: monoclinic, tetragonal, cubic. Now, the optimizations of parameters for structure cubic and tetragonal are also performed to study their stabilities.

12 atomic supercell is composed with 4 motifs of 3 atoms with space C. J. Howard 1988].

: Cubic Zirconia structure

The influence of the different k-points grids in the Brillouin zones and different cut energy of the 12 atomic supercell from 3×3×3 to 9×9×

After relaxation, the value chosen for parameter calculations of cubic zirconia is k point grid

The optimized parameters are listed in Table the other DFT calculation results and experiment data, they have good correlation.

84 to be a goodchoice in this part, it was and cubic) has 3 types of structures: monoclinic, tetragonal, cubic. Now, the optimizations of parameters for structure cubic and tetragonal are also performed to study their stabilities.

is composed with 4 motifs of 3 atoms with space in the Brillouin zones and different cut-off 9×9×9 are showed in After relaxation, the value chosen for parameter calculations of cubic zirconia is k point grid in Table 2. Comparing with the other DFT calculation results and experiment data, they have good correlation. For ZrO 2 tetragonal, the 6 atomic tetragonal bulk is composed of 2 motifs of 3 atoms with space group P42/nmc, a=b≠c (Figure 8). The oxygen atoms were displaced by pair along the z direction To study the relative stability of the three ZrO phases were also calculated and compared by using the formula:

where:

• @Yr+s$l+ is the cohesive energy of supercell;

•
)/ is the atom numbers of Zr in the supercell;

. The report of c/a is set to be 1.414 before optimization : Tetragonal Zirconia structure in the case of bulk cubic ZrO 2 , the convergence of the total energy was points meshes and different cut-off energies. study the relative stability of the three ZrO 2 structures, the cohesive calculated and compared by using the formula: So as showed in Table 4, the relative stability of the three phases are in the flowing order:

@Yr+s$l+ = sW¢+/@+oo -< )/ )/
monoclinic> tetragonal > cubic
The resultsarein good agreement with the experimental and theoretical ones, showing that our calculations can correctly describe the ZrO 2 structures .

As discussed above, ZrO 2 could adopt various phases in different ranges of temperature and pressure. In particular, the three low pressure phases: cubic, tetragonal, and monoclinic structures have been extensively studied for both their scientific and technological importance 1996]. For the cubic structure, it is only thermodynamically stable at temperatures higher than 2650 K [O. Ruff 1929]. All the transition results are ina good agreement with the stability evaluation by using Eq.1. Therefore, the phase transition mechanisms have been well

understood through modeling studies and experimental observations.

I.1.CMonoclinic zirconia supercell construction

Since there is low concentration of carbon in monoclinic zirconia, to simulate a carbon atom behavior, it is necessary to construct a larger size of zirconia supercell. Moreover, to perform calculation in constant volume, the relaxation of the supercell has to be refigured when one carbon atom is inserted.

To compare the results with other literature calculation ones [A. [START_REF] Foster | [END_REF], G. Jomard 1998], the most used supercell in the references contains 96 atoms and only 3.27% volume relaxation was chosen to perform the calculations (Figure 9). It was extendedfrom the 12 atomic supercell by two in each directions.

I.1.D Investigation of carboninterstitials in the supercell of zirconia

Once building a supercell of monoclinic zirconia, a carbon atom isinserted to investigate its possible interstitial sites to study its diffusion behavior. In the Chapter II, different pseudopotential sets were used for the ZrC structure, and Zr_4/C_s pseudopotential set was found to correctly describe the Zr-C bond. So in the simulation of carbon in bulk monoclinic zirconia, the same pseudopotential of C is used: C_s.

In order to know the stability of a carbon in the model of ZrO 2 monoclinic, the possible positions to insert it should be investigated.

At first, a carbon is inserted in the center of different available volume cages and optimized to find its minimum energy state. But it was found that for one cage, since the anisotropic character of monoclinic structure through the three directions x, y and z, there are different interstitial sites in one volume cage. By this method, the possible interstitial sites could not be completely found. So the investigation of interstitial sites should be performed by the way of a code.

Therefore, a code was developed and performed on the 12 atomic supercell in the 96 atomic supercell to find all the possible interstitial sites. The 96 atomicsupercell was equally separated by 20 times in each direction. Since the supercell is made by a 12 atomic cell extended by 2 in each direction, the 12 atomic cell is separated by 10 (Figure 10), and each grid is a cubic cage with dimension about 0.5 Å×0.5 Å×0.5 Å. First, a carbon was put in a random position of each small cubiccell with reasonable distance C-O (>1,3 Å) and C-Zr(>1.9 Å), the cut-off energy is set to be 300 eV and k points mesh is 1×1×1, then it was optimized to get the total energy convergence and geometrical parameters convergence.

Then at the second step, the optimized carbon was displaced the direction x, y, and z. The carbon in the new direction will be optimized again. If it gets back to the optimized position in the first step, it would be kept to verify its stability in the next step. If it goes to another posi

The third step is the last but also the most complicated step for carbon to verify if i possible interstitial position sites.

frequencies were calculated. If interstitial site position for a carbon So by all the three steps performed by the found. All their total energies and struct 96 atomic supercellfor carbon in the 78 positions are presented.

Figure 11: The founded positions with their total energies

There are 10 different system energies in the Figure 12. The three lowest ones have closed energies, similar structures, and they are closed within distance smaller than 0.2 Å. So to make the simulation simpler, only the most stable sites were kept, and th interstitial sites were considered to be the same position as the most stable ones (red color in Figure 11). After neglecting the two types of interstitial positions, the types of interstitial sites are still as many as eight. For each type of interstitial carbon is analyzed in detail.

Then at the second step, the optimized carbon was displaced by about 0.2 Å separately by the direction x, y, and z. The carbon in the new direction will be optimized again. If it gets back to the optimized position in the first step, it would be kept to verify its stability in the next step. If it goes to another position, then it would be deleted.

The third step is the last but also the most complicated step for carbon to verify if i possible interstitial position sites. For all the relaxed carbon position

. If their values are all real (without imaginary value) carbon atom.

So by all the three steps performed by the code, 78 positions in the 12 atomic found. All their total energies and structures were analyzed. In Figure 11, the carbon in the 78 positions are presented.

: The founded positions with their total energies

There are 10 different system energies in the Figure 12. The three lowest ones have closed energies, similar structures, and they are closed within distance smaller than 0.2 Å. So to make the simulation simpler, only the most stable sites were kept, and th interstitial sites were considered to be the same position as the most stable ones (red color in Figure 11). After neglecting the two types of interstitial positions, the types of interstitial sites are still as many as eight. For each type of interstitial position, the structure around the interstitial carbon is analyzed in detail. 90 about 0.2 Å separately by the direction x, y, and z. The carbon in the new direction will be optimized again. If it gets back to the optimized position in the first step, it would be kept to verify its stability in the The third step is the last but also the most complicated step for carbon to verify if it is a positions, the vibration (without imaginary value), it will bean positions in the 12 atomic supercell are , the total energies of There are 10 different system energies in the Figure 12. The three lowest ones have closed energies, similar structures, and they are closed within distance smaller than 0.2 Å. So to make the simulation simpler, only the most stable sites were kept, and the other two interstitial sites were considered to be the same position as the most stable ones (red color in Figure 11). After neglecting the two types of interstitial positions, the types of interstitial sites interstitial position, the structure around the Eight distinct possible sites have been found in the monoclinic structure (Figure 12) and for each site there are four symmetric positions in the 12 atomicsupercell which have the same structure around carbon but in different orientations. Here as example, the four symmetrical positions of site 2 are presented in Table 5.

Estimation of the formation energy of an interstitial carbon atom was calculated by the formula: (2)

The results of interstitial formations energies for the eight types of sites were listed in Table 6. The insertion formation energies are positive, which means E C-ZrO2 > (E C +E ZrO2 ). So considering diamond carbon as ground state, the insertion of an interstitial carbon is not favorable in the perfect ZrO 2 monoclinic bulk.

I.1.EMigration of C between two neighbor sites in the monoclinic zirconia without defects

Now, all the possible interstitial sites were investigated. In the monoclinic zirconia without defects, to calculate the carbon diffusion coefficient in bulk, the diffusion by interstitial mechanism will be determined. Since the diffusion at atomic scale is atomic jumps, here in this part, the migration of carbon between different neighbor sites will be calculated. Starting from the interstitial sites investigated in the last part, the first step is to identify all the neighbor sites for each interstitial site. Since there are four symmetrical positions for one type of sites, for each interstitial site position, there are many neighbor sites. In the atomic jump, to make the migration simple to calculate, for one interstitial position, among all its neighbor sites, the atomic jump study will concentrate only the nearest neighbor position of each type, which means there are eight neighbor positions from the eight types of sites. So within the 12 atomic supercell, all the neighbor sites are identified with their corresponding migration vectors (Figure 13). Then for each two neighbor sites, the NEB calculations are performed to test the possibility of migrations.

During the NEB calculation, the migration energies (E mig ) from one site to its neighbor site are firstly calculated with the migration vectors. For each migration, 8 imagepositions were built from the relaxed initial to thefinal equilibrium position. After the optimization the saddle points have been determined. It is found that for two neighbor sites, sometimes there is no direct migration path (the carbon would prefer to go to another neighbor site first), so it is not possible to do a direct atomic jump. For all atomic jumps, to study their jump probabilities, the attempt frequencies were also calculated, as in Chapter II. First, the KMC calculation wasperformed from 300 K, it was found that carbon could not diffuse in three dimensions until 2500 K. So to study the carbon diffusion, the carbon's to 3500 K were calculated. The diffusion coefficients at each simulation temperature were calculated and listed in Table 9.

D (cm 2 /s) lnD (Å 2 /s) 6.38×10 -6 24.88 9.61×10 -6 25.29 1.60×10 -5 25.80 3.09×10 -5 26.46

ifferent diffusion coefficient at the simulation temperatures of different diffusion coefficients as function of temperatures in the Figure 14.

he diffusion coefficient of C in the perfectzirconia 95 attempt frequencies macroscopic scale, the Kinetic using the residence time algorithm can be used to monoclinic zirconia.

performed from 300 K, it was found that carbon could not diffuse in three dimensions until 2500 K. So to study the carbon diffusion, the carbon's

The diffusion coefficients at each as function of temperatures for C in pure So it corresponds to the formula, the prefactor D 0 equals to 1.6×10 -3 cm 2 /s and the effective activation energy E act is 1.19 eV, leading to:

OE = 1.6 × 10 XG (žŸ /š) X ….…¥ 2¦ § * ¨(3)
At the real storage condition, the higher final temperature is around 50 ᵒC, the carbon diffusion in the perfect zirconia bulk at the storage condition was calculatedas ~10 -22 cm 2 /s.

As example, in the claddings where the zirconia layer is about 50 µmof thickness, if a carbon atom is in the center of the oxide layer, it will take to the carbon 10 11 years to go out from the zirconia layer. It provides that the carbon does not diffuse in the perfect ZrO 2 bulk at the storage temperature.

II Carbon behavior in zirconia with point defects

The carbon diffusion in perfect zirconia was calculated in the last part. However, in the real claddings, in the oxidation layer, there are point defects. So it is interesting to complete this study to analyze the vacancy effect on the carbon diffusion in zirconia.

II.1Vacancyand carbon substitution in the zirconia supercell

II.1.A Vacancy of O1, O2and Zr in monocliniczirconia

There are two types of oxygen in the ZrO 2 monoclinic: O1 which is three coordinated with zirconium and O2 which is four coordinated with zirconium (Figure 2). So a vacancy of O1, O2, and Zr would be first considered in the supercell of 96 atoms. The formation energies to form a vacancy werecalculated by the formula: lp@p#@i? = @+oo7q$Br lp@p#@i?: ? -@+oo (4)

where:

• lp@p#@i? : the formation energy for the supercell to form a vacancy x; • @+oo7q$Br lp@p#@i?: :the total energy of supercell with a vacancy x; • ? : the atom energy of x, where x is the atom of vacancy; • @+oo : the total energy of the pure supercell.

For the three kinds of carbon substitutions inO1, O2, and Zr positions, the formation energies are 3.06 eV, 3.33 eV, -7.02 eV respectively (Table 11).

Substitution energy sub Zr -7.02 eV sub O1

3.06 eV sub O2

3.33 eV So the carbon atom insertion in the position of a zirconium vacancy is very favorable. After comparing these formation energies with the interstitial energies (6.05 eV to 7.76 eV), it turned out that a carbon atom always preferred to be in a vacancy position than in an It can be seen that when there was no substitution, a Zr has 6 neighbor oxygen atoms. After it was substituted by a carbon atom, the carbon atom has only 2 neighbor oxygen atoms with formation energy -7.02 eV.

Since it is proved that carbon is stable as substitution of zirconium and there is more space in the substitution site, it should be tried to insert more than one carbon atoms in the site of a zirconium atom.The formation energies calculated by the Eq. ( 5) from one to four substitution carbon atoms were listed in Table 12. Since it is proved that carbon is stable as substitution of zirconium and there is more space in the substitution site, it should be tried to insert more than one carbon atoms in the site of a zirconium atom.The formation energies calculated by the Eq. ( 5) from one to four substitution carbon atoms were listed in Table 12.

N br C 1 2 3 4 E sub (eV) -7.02 -5.53 -1.84 0.89 As we can see in Table 12, the formation energy is negative until the fourth carbon atoms was inserted in Zr vacancy site, which means the supercell was more stable when the carbons were in the vacancy site. All the carbon substitutions in the Zr vacancy site formed a framework in the space (Figure 16). Comparing with the formation energy of carbon in interstitial position, the inserted carbon atoms still preferred to be in Zr vacancy site than be in interstitial position. So if there are Zr vacancies in the zirconia, the carbon atoms will be trapped in the vacancy sites and it won't diffuse.

II.1.C Migration by vacancy mechanism

In the previous part, the carbon is found to be more stable in the vacancy position than in an interstitial position. When there is a Zr vacancy, the carbon atom will be trapped in the vacancy. When there is an O vacancy, the carbon atom is more stable in the vacancy position than being in the interstitial sites. So it is necessary to test the possibility of carbon diffusion by oxygen vacancy mechanism in bulk zirconia. The distances between the O1 vacancy and neighbor O1 atoms were noted in red color in atoms, the two pairs of symmetrical atoms are much easier to migrate to the O1 vacancy site.

2)Migration to a vacancy O2

As example, if the atom No.48 was taken away to form an O2 vacancy, there were 9neighbor atomsin the radius of 3.5 Å: five neighbor O1 atoms with distances noted in red color, and four neighbor O2 atoms with distances noted in green colorwhich includes two symmetrical atoms No.45 and No.77 (Figure 18, Table 14). To compare with the migration of O without carbon, still the O1 vacancy No.35 was considered with all its neighbor positions. The binding energies for C in the neighbor positions were calculated and listed in Table 15with their corresponding distances to the O1 vacancy.

All the binding energies have positive values, meaning that the carbon is more stable in neighbor O positions of the O1vacancy site than being in the position far away from the vacancy.

In a second time, the migrations of a neighbor carbon atom from the different neighbor positions were calculated. The energy barrier for carbon in some neighbor positions does not exit since there is no direct migration path ( The binding energy decreaseswhen the distanceincreases. Since the binding energy calculated here is always larger than 2 eV, the DOS calculations were performed to test the electron structure effect on the binding energy.

II.2 Electronic structure: Density of states

In the last part, the different migrations mechanisms of O or C have been determined. Since the competition is complicated and also the position of C have a large effect on the energy barrier of migrations, to better understand the total energy of bulk with different defects, the density of states are performed.

At the first step, the DOS of the 96 atomic supercell will be calculated as reference and also will help to understand the composition of the density of states.

Then, at a second step, three DOS of supercells with different defects were calculated. First a supercell with a carbon atom in interstitial position was considered to study the effect of inserted carbon on the electron density. Then the DOS of a supercell with an O vacancy was calculated to see the vacancy effect on the electron density. Since the carbon atom prefers to go to vacancy position than the interstitial site, the DOS of supercell with a carbon substitution in O vacancy was calculated to evaluate the electronic density.

At last, the DOS of supercell which contains a C su calculated.The binding energy between neighbor C substitution and the O vacancy will be discussed.

II.2.A DOS of the perfect 96 atomic

The DOS were calculated in the structure of 96 atoms without carbons or defects. After cell optimization, the file POSCAR was replaced by CONTCAR, and a larger k The number of electrons in the material is the integral of the density of states (Figure 26).

As presented in Figure 26, for the supercell which contains 96 atoms (32 Zr and 64 O), the electron number is calculated to be 32x4+64x6=512.

II.2.B DOS of the 96 atomic supercell of monoclinic zirconia with defects

All defect calculations were made using a 96 atomic supercell, which wasgenerated by extending the 12 atomic monoclinicsupercell by two in all three dimensions. For the optimized supercell, the four types of models which are: 1) supercell with an interstitial carbon atom, 2)supercell with an O vacancy, 3) supercell with a carbon substitution,4)

supercell with a carbon substitution and an O vacancy will be calculated. During defect calculations the lattice vectors of the supercell were frozen which means the volume wasconstant and all atoms were allowed to relax until atomic forces were converged with a variation smaller than 0.001 eV/Å. The change of electrons numbers caused by the defects is very small comparing to the total electron numbers of the supercell (512 electrons for 96 atomic supercell), so it is sometimes very difficult to get the DOS peaksof defects. At this time, it is necessary to do the derivation of the integral density to get a more accurate DOS.

As explication in Chapter I, the derived of the integral density of states is defined as below:

= k+/$l+k = k' 2 7j: kj (6)
By derivation of the integral DOS, any change of electron numbers will be shown..

B-1 DOS of 96 atomic supercell with a carbon interstitial atom

DOS for the 96 atomic supercell with a carbon in interstitial position was calculated to see the change of electron density caused by a carbon atom in interstitial position. Here is DOS plotted for the supercell when there is a C atom in site 8 which is the most stable interstitial position (Figure 27). At this time, the total DOS plot did not present the peaks caused by defects while the derived integral DOS plot showed the new peaks.For static calculations, the file PROCAR contains the spd and site projected wave function character of each band.Comparing with the total DOS of perfect supercell, the compositions of the different peaks can be analyzed from the PROCAR (Table 19). In all the analyses by PROCAR, the spd site projected wave function character of the band located at the peak value at Gamma point were listed in the To identify the peaks, the partial DOS of Zr and O are plotted (Figure 29). It was found that the two electrons presented in Figure 28 are Zr_4d offered from the neighbor Zr atom of the vacancy.

When there is an O vacancy, the valence band was below the Fermi level. The Fermi level corresponds to the last level which was occupied by the electrons of the system. A neutral oxygen vacancy leads to the appearance of a DOS peak in the gap because in zirconia, the chemical bonds have ionic character, and the O atom was charged by two electrons. So the created vacancy was charged and the zirconia contains plus two electrons.

B-3 DOS of 96 atomic supercell with a carbon substitution in O2 position

Since a C atom is more stable in vacancy position than in an interstitial site, the DOS of supercell when there is a C in O2substitution position is calculated(Figure 30). To understand the total DOS, the partial DOS for a neighbor Zr and a neighbor O of the carbon substitution was alsoperformed(Figure 31). The compositions of different peaks caused by defects were analyzed in

B-4 DOS of 96 atomic supercell

To study the vacancy effect on the substitution, an oxygen vacancy was made with a distance 3 Å, 5 Å and 7 Å from the Å are performed in the 96 atomic supercell. A larger supercell of 324 atoms wa perform the simulation of distance 7 Å (the 12 atomic supercell each direction x, y and z). The total DOS DOS of the supercell in which the (Figure 32). The analyzed results of the different peaks using PROCAR To study the vacancy effect on the substitution, an oxygen vacancy was made with a distance 3 Å, 5 Å and 7 Å from the substituted carbon. The calculations for distance 3 Å and Å are performed in the 96 atomic supercell. A larger supercell of 324 atoms wa perform the simulation of distance 7 Å (the 12 atomic supercell was expanded three times in each direction x, y and z). The total DOS for the three supercell models were DOS of the supercell in which the vacancy is 5 Å away from the substitution : The total DOS of zirconia supercell with a carbon substitution and away from the substitution(derived integral DOS plot)

The analyzed results of the different peaks using PROCARand DOSCAR To study the vacancy effect on the substitution, an oxygen vacancy was made with a carbon. The calculations for distance 3 Å and 5 Å are performed in the 96 atomic supercell. A larger supercell of 324 atoms was built to For the supercell with distance 5 Å between carbon in substitution and the O2 vacancy, the binding energy was calculated to be 2.21 eV by Eq. 9. For distance at 3 Å and 7 Å, the binding energies calculated were 3.25 eV and 1.99 eV (cf. Table 22).

The binding energy decreased when the distance increased. However, there should not bea binding energy between vacancy and substitution when the distance is as far as 7 Å. But the calculated binding energy presented here is as large as 1.99 eV.

So to understand where 1.99 eV comes from, the binding energy was also calculated with charged vacancy. Indeed, for the 324 atomic supercell, when an oxygen vacancy is neutral, the binding energy is 1.99 eV. But if an O vacancy is charged in the calculations which means the supercell with an O vacancy loses two electrons(losing 2 electrons since one O 2-vacancy) the binding energy calculated for distance 7 Å is now 0.04 eV. So that means 1.99 eV comes from the neutral oxygen vacancy. There is no binding energy when the carbon atom is 7 Å far away from the vacancy.

With the same method, for distance 5 Å between the C substitution and the O vacancy, if an O vacancy is considered to be charged, the binding energy is 0.12 eV. The corrected binding energy proves that there is no vacancy effect on the carbon atom in substitution position when the vacancy is5 Å away from the substitution. So when the carbon is further than 5 Å to the vacancy, the vacancy effect on the migration of C can be neglected. For distance 3 Å, the corrected binding energy with the same method is 1.19 eV, so there is still interaction between the carbon substitution and the O2 vacancy. Since for anO vacancy, the neighbor positions were in radius of 3.5 Å, the energy difference 2.06 eV caused by neutral vacancy for distance 3 Å is chosen to correct all the binding energies of the neighbor atoms to the vacancy.

The binding energies in Part II.1.C-2 are correctedand listed in Table 23 By correcting binding energies presented in these two tables, the binding energies between neighbor C and O vacancy get smaller. When the binding energy is very small to be neglected, the carbon substitution has no interaction with the O vacancy, so it is more difficult to do the migration to the O vacancy.

Conclusions

Since there is oxidation reaction between zirconium and water, an oxide layer which is principally zirconia is formed (about 50 µm at the external layer). Zirconia at the storage temperature exists in monoclinic phase, so the monoclinic zirconia is chosen to do the simulations in order to study carbon behavior in it.

Density Functional Theory, in static approach, was used to evaluate the carbon behavior in the zirconia bulk. Cell models which are the smallest cells to reproduce the structure of three different structures for zirconia were built, such as cubic, tetragonal, and monoclinic. It is proved that monoclinic phase was the most stable one of the three phases, which corresponds to the phase transitions of zirconia at different temperature.

For the monoclinic phases, the possible stable sites for carbon atoms were characterized in the perfect bulk without any vacancies. By doing the Monte Carlo simulations for C diffusion by interstitial mechanism, it was found that the carbon atom does not diffuse in the perfect monoclinic zirconia at storage condition.

To simulate more realistic conditions, the formation energies of a vacancy, a substitution, and an interstitial carbon with a vacancy were studied. According the results, it was shown that it is more favorable for carbon atoms to go to the position of vacancies as substitution than in the interstitial sites.

The migration of O, the migration of C and their competition by vacancy mechanism have been performed. Comparing with the migration of C by interstitial mechanism (without vacancies), the migration energies of the carbon atoms are comparable. However, it was shown that a neighbor carbon atom of an O vacancy increases the migration of neighbor O atomsto the vacancy, leading to a stronger competition between the oxygen migration and the carbon migration to the same vacancy.Therefore, we can say that the carbon diffusion by vacancy mechanism is more difficult than the carbon diffusion by interstitial mechanism. As at the storage temperature, it was shown that there is no carbon diffusion by interstitial mechanism, we can say that there would be no carbon diffusion by vacancy mechanism.

To study the carbon behavior in the used fuel claddings, a multi-scale approach has been developed and performed for carbon atom in bulk zirconium and bulk monoclinic zirconia.

Behavior of a carbon atom in bulk zirconium

For the carbon atom in zirconium, first, with the DFT atomic scale simulation, a 96 atomic hcp zirconium supercell was optimized to perform the diffusion simulation. Two types of possible interstitial sites for the carbon atom were characterized: distorted Basal Tetragonal site noted as BT_d site and Octahedron site noted as O site. BT_d site is less stable than the O one by 1.45 eV. To study the atomic jump between the different neighbor interstitial sites, all the migration energies and the attempt frequencies were calculated in order to get the corresponding jump probabilities from the transition state theory for all possible atomic jumps.

As a preliminary result, calculation from the DFT-NEB approach showed that the easiest migration pathway for carbon diffusion should be the successive jumps between BT_d site and O site.

In a second step, a macroscopic simulation by the Kinetic Monte-Carlo approach (KMC) was performed using jump probabilities. From the KMC simulations at different temperatures, the carbon's trajectory as well as the diffusion coefficient (as a function of temperature) were The carbon diffusion coefficient in perfect Zr crystal, at the storage temperature (around 50°C), is calculated to be very small (~10 -34 cm 2 /s). Therefore, these first investigations show it will be very difficult for a carbon atom to diffuse in a perfect Zr crystal in the range of the storage condition.

In order to investigate the vacancies effect on the carbon diffusion, DFT simulations were also performed with vacancies. It was proved that the carbon atom stability decreases when it gets closer to a Zr vacancy. Therefore, bulk zirconium with vacancies won't increase the carbon diffusion relative to interstitial diffusion mechanism without vacancies.

To summarize, this first part has pointed out that the carbon will be blocked in zirconium bulk with and without Zr vacancies.

Behavior of a carbon atom in the zirconia

For the oxide layer of the fuel claddings (about 50 µm at the external layer), it is zirconia in monoclinic phase that principally exists at the storage temperature, so it was chosen to do the simulations in order to study carbon behavior in it.

At atomic scale, using Density Functional Theory, first, an optimized bulk model was built for the three different structures for zirconia, such as cubic, tetragonal, and monoclinic. It was proved that monoclinic phase was the most stable one of the three phases which corresponded to the phase transitions of zirconia at different temperatures.

For the pure monoclinic phase, the possible interstitial sites for carbon atoms were investigated. There are eight types of different interstitial sites for carbon atoms and each type has four symmetric positions. For each interstitial site position, all the neighbor positions were identified. By studying the possible atomic jumps, it is identified that there are 28 possible atomic jumps between the different neighbor interstitial positions.

The migration energy and attempt frequency were calculated for each atomic jump to get the jump probabilities. Then, using the atomic results to do the macroscopic simulations by the Kinetic Monte Carlo approach, it was found that it is very difficult for a carbon atom to diffuse by interstitial mechanism in monoclinic zirconia structure with effective activation energy of 1.19 eV. As example, a carbon atom in the center of the oxide layer of 50 µm will need 10 11 years to leave it.

At last, the vacancy effect on the carbon diffusion was also considered. Firstly, the binding energy between carbon atom in an interstitial position, and an O vacancy were calculated which indicated that the closer they were, the larger the binding energy was. It means that the carbon atom prefers to go near a vacancy position. So the study of the carbon diffusion by vacancy mechanism was performed. We have demonstrated that:

1) For the self-migration of oxygen atoms, the migration energies of different vacancy positions are between 1.3 eV and 3 eV, which have low jump probabilities.

2) The binding energies when a carbon atom is in different neighbor positions of an O vacancy were calculated. It was found that a carbon atom in substitution was more stable in neighbor position of an O vacancy site. So the possibility of carbon diffusion by vacancy mechanism was calculated. The migration energies of C were calculated for all the possible atomic jumps which turned out between 0.17 eV and 2.03 eV.

3) Since there is a positive binding energy between an O vacancy and the substituted carbon atom near to the vacancy, and the substituted carbon changes the geometrical structure around the vacancy, the migration energies of neighbor O atoms would be affected. All the possible atomic migrations of O to the vacancy were recalculated with presence of a neighbor C atom. It is found that the migration energies get smaller when there is a carbon atom nearby the vacancy. The migration energies were calculated between 0.27 eV and 3.20 eV. Since the migration energies are comparable to the migration energies of C (between 0.17

  .. I.2 Macroscopic scale simulations ..................................................................................... II Carbon behavior in zirconia with point defects ............................................................ II.1 Vacancy and carbon substitution in the zirconia supercell .......................................

  .. 1) Migration to a vacancy O1 .............................................................................................. 2) Migration to a vacancy O2 .............................................................................................. C-2 Migration of C by vacancy mechanism .......................................................................... 1) C substitution in a neighbor O position of an O1 vacancy .............................................. 2) C substitution in a neighbor O position of an O2 vacancy .............................................. C-3 Migration of O with a neighbor carbon to a vacancy: competition effects .................... 1) Carbon in a neighbor O1 position of the O1 vacancy ..................................................... 2) Carbon in a neighbor O2 position of the O1 vacancy ..................................................... 3) Carbon in a neighbor O1 or O2 position of an O2 vacancy ............................................

  structure: Density of states ..................................................................... II.2.A DOS of the perfect 96 atomic supercell of monoclinic zirconia ......................... II.2.B DOS of the 96 atomic supercell of monoclinic zirconia with defects ................. B-1 DOS of 96 atomic supercell with a carbon interstitial atom ......................................... B-2 DOS of 96 atomic supercell with an O2 vacancy ........................................................... B-3 DOS of 96 atomic supercell with a carbon substitution in O2 position ........................ B-4 DOS of 96 atomic supercell with a carbon substitution in O2 position with an O vacancy

Figure 1 :

 1 Figure 1: Schematic view of PWR fuel assembly (Mitsubishi Nuclear Fuel)

Figure 2 :

 2 Figure 2: The zirconium layer and oxide layer (ZrO

Figure 2 :

 2 Figure 2: The zirconium layer and oxide layer (ZrO 2 ) of thefuel cladding

  Figure 3: The process of nuclear waste treatment (a) cut; (b) compacted; (c) introduced in container CSD-C

Figure 4 :

 4 Figure 4: View of a concrete package to enhanced retention capacity storage of CSD-C containers

Figure 1 :

 1 Figure 1: Pseudopotential of atoms

  off energy which corresponds to the different off radius r is used to separate the valence electrons from the core electrons. The VASP principles pseudopotentials are derived from an atomic reference state, requiring that the pseudo and all-electron valence eigenstates have the same energies and amplitude (and thus density) outside a chosen core e is simplified by Ψ pseudo to calculate , the exact wavefunction is Many different basis sets have been proposed in the literature historically, and the two most -Spherical Bessel-functions (Rappe-Rabe-Kaxiras-Joannopoulos (RRKJ) method (Ref A.

Figure 2 :Figure 3 :

 23 Figure 2: The calculation process of VASP

Figure 4 :

 4 Figure 4: The reciprocal lattices (dots) and corresponding first Brillouin zones of (a) square lattice and (b) hexagonal lattice.

4:

  The reciprocal lattices (dots) and corresponding first Brillouin zones of (a) square In order to define the Brillouin zone we need to define first the reciprocal lattice. The set of that yield plane waves with the periodicity of a given Bravais lattice is belongs to the reciprocal lattice of a Bravais , and all R in Bravais (11) Factoring out, we can characterize the reciprocal lattice as the set of wave vectors k (12)

Figure 5 :

 5 Figure 5: The different point defects in crystal

Figure 6 :

 6 Figure 6: Binding energy calculation

Figure 7 :

 7 Figure 7: Principle of NEB method [D. Sheppard 2008]

Figure 8 :

 8 Figure 8: The comparison of NEB and Climbing Image NEB (CINEB) [G.Henkelman, 2000]

Figure 9 :

 9 Figure 9: The migration energy barrier a) E migration from site 2 to site 1

Figure 10 :Figure 9 :

 109 Figure 10: The atomic possible jumps from the initial state to the final states

Figure 10 :

 10 Figure 10: The atomic possible jumps from the initial state to the final states

Figure 12 :

 12 Figure 12: The mean square displacement of carbon atom at temperature T

Figure 13 :

 13 Figure 13: The linear relation of lnD and temperature T

  Introduction ............................................................................................................................ I: Carbon behavior inzirconium ......................................................................................... I.1 Atomic scale simulations ............................................................................................. I.1.A Optimization of the calculation parameters .

  Zr bulk super-cell size ............................................................. I.1.C Investigation of the carbon interstitial sites ............................................................ I.1.D Migration between neighbor sites ......................................................................... D-1 Migration energies (NEB) .

  scale simulations ....................................................................................I.2.A KMC simulations .................................................................................................... I.2.B Determination of the Diffusion coefficient ............................................................. I.3 Conclusion .................................................................................................................... II: Carbon behavior in zirconium bulk with point defects ............................................... II.1. Zr vacancy in pure Zr bulk ......................................................................................... II.2. Substitution of Zr by carbon ..................................................................................... II.3 Interstitial carbon with a nearest Zr vacancy ............................................................ II.4 Vacancy effect on carbon diffusion by calculation of binding energy ..................... Conclusion ............................................................................................................................... References ...............................................................................................................................

Figure 1 :

 1 Figure 1: Zirconium HCP structure Figure 2: Unit

Figure 1 :

 1 Figure 1: Zirconium HCP structure Figure 2: Unit-cell of Zirconium HCP structure

  51numbered 194[B.Olinger 1973], with parameters: a=3.232 Å, b=3.232 Å, c=5.147 Å, α=90°, cell of Zirconium HCP structureThe crystal structure of Zr can be constructed by 2 atoms in one unit cell (Figure2), a motif initial structure was built, the corresponding 5s 2 . In this case, two Zr_4 pseudopotential, the last ) are considered as valence electrons. For the Zr_12 pseudopotential, the ) are considered as valence electrons. The pseudopotential al energy convergence and geometrical parameters with points and different cut-off off energy with a 11×11×11

Figure 3 :

 3 Figure 3: Unit-cell of C diamond structure

Figure 4 :

 4 Figure 4: Unit-cell of ZrC structure

  energy of 450 eV, and the geometrical parameters converge with a cut-off energy of 350 eV By comparing the convergences of total energy and geometrical parameters as a function of points grids, similar results are found. Therefore, (Zr_4/C_s) can be off of 350 eV and a 9×9×9 k-points grid in the later calculations since its lower calculation time cost and good description of the structure. The parameters a, b, nd c are 4.725 Å, which have a good agreement with experimental results 4.68 Å [A. C. Lawson 2007] and theoretical values 4.72 Å [F. Vines 2005]. 55 Two combinations of pseudopotentials are tested in this part. We study one ZrC structure using the second set of pseudopotentials (Zr_4/C_s). By comparing their total energy convergence and geometrical points, we can determine which set of points with a cut-off energy off energy of 450 eV and kpoints with a cut-off off energy of 350 eV By comparing the convergences of total energy and geometrical parameters as a function of s grids, similar results are found. Therefore, (Zr_4/C_s) can be points grid in the later calculations since its lower calculation time cost and good description of the structure. The parameters a, b, nd c are 4.725 Å, which have a good agreement with experimental results 4.68 Å [A. C.

Figure 5 :

 5 Figure 5: ZrC total energies and geometrical parameters with different values of k-points and cut-off energies, for the (Zr_4/C) and (Zr_4/C_s) pseudopotentials

Figure 6 :Figure 7 :

 67 Figure 6: Zr hexagonal close-

Figure 8 :

 8 Figure 8: Supercell 4×4×3 of Zr with one interstitial carbon atom in (O) site (grey atom)

  panel b), and a distorted basal tetragonal one named as BT_d site (Figure 9 panel c). For the BT_d site, the neighbor Zr atoms were displaced to get a stable state of the structure. The distances between Zr and C were changed relative to an ideal BT site (Fig. 9 panel a).

Figure 9 :Figure 9 :

 99 Figure 9: The interstitial sites for C in the Zr model composed of 96 Zr atoms. a: basal tetragonal site (BT); b: octahedral site (O); c: BT_d sites one BT site)

Figure 10

 10 Figure 10 panel b, if C migrates from BT_d site to O site, it needs 0.70 eV to cross the migration barrier. However, it is much harder to migrate from O site to BT_d site, which needs 2.13 eV to pass the barrier. The migration path between neighbor BT_d site and O site is always following the path: BT_d site BT_d site (n times) O site BT_d site …

Figure 10 :

 10 Figure 10: Migration between two neighbor sites for C in Zr bulk (Zr_4). a: migration path from BT_d1 site to BT_d2 site; b: migration path between BT_d site and O site; c: migration path between two neighbor O sites

Figure 11 :

 11 Figure 11: Example of C's atomic trajectory at T=900K for 100000steps

Figure 12 :

 12 Figure 12: The mean square displacement of carbon atom at 900 K.

Figure 13 :

 13 Figure 13: linear relation between lnD and 1/T from KMC simulations

  the carbon trajectory, the carbon diffusion pathway is through: O site BT_d site BT_d site BT_d site (n times) O site….

Figure 14 :

 14 Figure 14: A C in interstitial site (gray color) with a neighbor Zr vacancy (blue cycle) in the 95 zirconium supercell. a: C in O site with a neighbor Zr vacancy; b: C in BT_d

  was calculated. It was found that 2.02 eV effective activation energy corresponded to the lowest migration energy path:

  done in Chapter II, Part I.1. Both Zr_4 and Zr_12 can well describe the hcp structure of zirconium except Zr_12 which demands higher cut-off energy to convergebothforce and energy. So the next step is to choose the pseudopotential set of Zr/O to better describe the monoclinic zirconia. For pseudopotential of O, two kinds of pseudopotentials were tested: O_s and O_h which both contains 6 valence electrons. Since O_h demands an higher cut-off energy to arrive the convergence of geometrical parameters, O_s is better to describe the O 2 molecule. For ZrO 2 , different pseudopotentials were compared separately for Zr and O. So now it is important to estimate the different pseudopotential set to verify which one can better describe the Zr-O bond. Monoclinic zirconia has a structure with space group P21/c (Figure 1)[G. M Rignanese 2005, C. J. Howard 1988].The model of 12 atoms was built with atomic coordinates in Wyckoff notation (x, y, z), (-x, y+1/2, 1/2-z),-(x, y, z),-(-x, y+1/2, 1/2-z).In the bulk structure, there are two types of O, the first has three bonds with zirconium, noted as O1, and the other one has four bonds with zirconium noted as O2. The zirconium is six coordinated by the oxygen atoms (Figure2).The convergence of the total energy were tested with different k-points and different cut-off energies for the 4 sets of pseudopotentials, (Zr_4/O_s), (Zr_12/O_s), (Zr_4/O_h), and (Zr_12/O_h).Different total energies with different geometrical parameters a, b, c were calculated, to evaluate the proper values of E cut-off and cell parameters which can well describe the monoclinic zirconia.

Figure

  Figure 1: The 12 atomic super

  /), the situations are similar. The cell parameters still do off energy is already as high as 700eV.(Zr_12/O_s), the total energy and the geometrical parameters of the cell begin to converge at E cut-off = 350 eV. There is no significant difference of the total points grids. The 5×5×5 k-points grid is chosen. Here are the figur of total energy and geometrical parameter for the pseudopotential set Zr_4/O_s (Fig 82

Figure

  Figure 3: The total energy of the zirconia unit (Zr_4/O_s)

Figure 5 :

 5 Figure 5: Cubic Zirconia structure

Figure

  Figure 6: The energy of a cubic structure unit

parameters of the cubic zirconia 6

 6 atomic cell which is the smallest structure to reproduce the is composed of 2 motifs of 3 atoms with space group P42/nmc, ). The oxygen atoms were displaced by pair along the z direction is the smallest structure to reproduce the is composed of 2 motifs of 3 atoms with space group P42/nmc, parameters ). The oxygen atoms were displaced by pair along the z direction which the distance is noted as dz. The report of c/a is 1988].

Figure 8 :

 8 Figure 8: Tetragonal Zirconia structure

[

  F. F. Lange 1992, G. Teufer 1962]. At ambient pressure, ZrO 2 has a monoclinic structure up to a temperature of 1450K [J. D. Mccullough 1959]. At temperature higher than 1450K, the monoclinic to tetragonal transition occurs with a decrease in unit cell volume since the 7-fold coordination in monoclinic phase of Zr changes to 8-fold in tetragonal phase[M. Wilson

Figure 9 :

 9 Figure 9: 96 atomic supercell of monoclinic zirconia extended from the 12 atomic cell

Figure 10 :

 10 Figure 10: The small grids for the monoclinic supercell (10 times of separation for the 12 atomic supercell)

Figure 12 :

 12 Figure 12: Eight types of interstitial sites in monoclinic zirconia

Figure 13 :

 13 Figure 13: Example of site 1's neighbor sites

  At the atomic scale, the migration vectors, migration energies and the have already calculated to study atomic jumps. Then Monte-Carlo method [H. Jon determine the diffusion path and diffusion coefficient in First, the KMC calculation was diffuse in three dimensions until 2500 K. So to study the carbon diffusion, the carbon's trajectories from 2500 K to 3500 K were calculated simulation temperature were calculated and listed in TableT

Figure 14 :

 14 Figure 14: The diffusion coefficient of C in the

  interstitial one. When there is an O vacancy, the carbon atom might diffuse by vacancy mechanism. The structure surrounding the carbon in the case of the substitution of type O1, O2 were almost like the structure before the insertion of carbon substitution.The substitution by C was three-coordinated or four-coordinated (Figure 15 panel a and Figure 15 panel b) and the relaxed distances are given. The C -Zr distances are around 0.1 Å larger than O -Zr distances. The configuration of Zr substitution is particular. After optimization, the carbon has a strong relaxation, and is not in the same position as the Zr vacancy but in the zone of the vacancy. The structure before and after the substitution are displayed (Figure 15 panel c).

Figure 15 :

 15 Figure 15: Structures before and after the carbon substitution

  before the substitution Structure around Cafter the substitution Structure around O2 before the substitution Structure around Cafter the substitution Structure around Zr before the substitution Structure around Cafter the substitution

Figure 16 :

 16 Figure 16: The framework of carbons in the position of a Zr vacancy (a) one C substitution; (b) two C substitutions; (c) three C substitutions; (d) four C substitutions.

  vacancy effect more clear, a model of an O vacancywith neighbor atoms were simulated (Because in the real condition, there are mostly O vacancies since the zirconia came from the oxidation reaction of zirconium.).C-1Self-diffusion of O by vacancy mechanism without carbonAt first, in the 96 atomic supercell, an O vacancy was included to study the O migration by vacancy mechanism. Since there are two kinds of oxygen atoms in bulk monoclinic zirconia: one is three coordinated with neighbor zirconium noted as O1 and the other is four coordinated with neighbor zirconium noted as O2, the oxygen migration to an O vacancy wereperformed and summarized in two parts: O1 vacancy and O2 vacancy. All the migration barriers for the neighbor oxygen atoms were calculated. 1) Migration to a vacancy O1 As example, if atom No.35 was taken away to form an O1 vacancy, there are 10 neighbor O atoms (Figure 17)in the radius of 3.5 Å. There are five O1 atoms includingtwo symmetrical neighbor atoms No.34 and No.66, two symmetrical neighbor atoms No.44 and No.60, and one neighbor atom No.41.

Figure 17 :

 17 Figure 17: The neighbor O atoms around an O1 vacancy (blue position)

Figure 18 :

 18 Figure 18: The neighbor O atoms around an O2 vacancy (blue position)

Figure 19 :

 19 Figure 19: the O1 vacancy with all the neighbor O atoms

Figure 19 :

 19 Figure 19: the O1 vacancy with all the neighbor O atoms (red) and the C substitution

Figure 20 :

 20 Figure 20: The migration energies of neighbor O atoms to an O1 vacancy when C is in different neighbor O1 positions (a) C in No.41; (b) C in No.66; (c) C in No.44 (no possible atomic jump for C)

(

  5×5×5) was used to do the DOS calculations. Only the electrons with energies in t 20 eV, 20 eV) weretaken into account in the calculations. The Fermi state was shifted at the 0 eV point, the total density of states (Figure an O atom were plotted (Figure 2 The DOS for the monoclinic phase is very similar to the DOS obtained for the cubic and tetragonal phases and has three clear bands. The total DOS is composed of two valence bands, where one is centered around -5 eV corresponding to O_2p principally and Zr_4d. They are separated conduction band by a gap of values (3.4 eV/3.42 eV) of the literatures [J. X. Zheng 2007, the value 5.83 eV found by Ultraviolet absorption experiments eV found by electron-energy-loss spectroscopy the energy gap calculated by DFT may display discrepancies with experimental data since this theory is strictly valid for the ground state properties.

Figure 24 :

 24 Figure 24: The total DOS of monoclinic ZrO

Figure 25 :Figure 26 :::

 2526 Figure 25: The partial DOS of one Zr atom (panel a) and one O atom (panel b)

Figure 27 :

 27 Figure 27: The total DOS of zirconia supercell with a carbon interstitial DOS plot)

Figure

  Figure 28: The total DOS of zirconia supercell with

Figure 30 :::

 30 Figure 30: The total DOS of zirconia supercell with a carbon

Figure 32 :

 32 Figure 32: The total DOS of zirconia supercell with a carbon substitution and 5 Åaway from the substitution

::

  The total DOS of zirconia supercell with a carbon substitution and an O vacancy at DOSCAR is presented in 0.25 eV are both composed by electrons from C_2p and Zr_4d. For the peak at Fermi level energy (0 eV), it is principally from the electrons of The compositions of different peaks in DOS of supercell with a carbon atom in O2

  calculated. It is found that the carbon diffusion in zirconium has 2.02 eV effective activation energy, which corresponds to the lowest migration energy path:

  eV and 2.03 eV), there is a competition between the migration of O and the migration of C.The simulations for the migration of O, the migration of C and their competition by vacancy mechanism have been performed. Comparing the carbon diffusion by two different mechanisms, the migration energies of C by vacancy mechanism (between 0.17 eV and 2.03 eV) are comparable with the activation energy of C diffusion by interstitial mechanism (1.19 eV). The neighbor O atoms of an O vacancy have competition with the carbon atom to migrate to the vacancy, which makes the C diffusion by vacancy mechanism more difficult. In the two ways, the carbon diffusion by vacancy mechanism is proved to be also extremely limited, the same as for the carbon diffusion by interstitial mechanism. As at the storage temperature, based on our calculations, and considering a half life of 5600 years for 14 C, the consequences are the following: After 100.000 years, C should have diffused on 10 -7 µm in Zr and 0.1 µm in ZrO 2 , and the activity should be reduced by 10 -6 in the meantime.After 200.000 years, C should have diffused on 10 -7 µm in Zr and 0.2 µm in ZrO 2 , and the activity should be reduced by 10 -11 in the meantime. Thus, after 100.000 or 200.000 years, and a diffusion distance on 0.1 to 0.2 µm, C should remain within the grain of zirconia (whose size should be around 1 µm size) and very few amounts may reach grain boundaries. Even if grain boundaries and micro-cracks, which were not studied in this work, are short circuits for diffusion (e.g. with very high Carbon diffusion coefficients), their impact on C release should remain small.

Table 1 : Composition of zircaloy-4 in the fuel cladding (cladding of ASTM B353)

 1

  Simeone 2002] D. Simeone, D. Gosset, J. L. Bechade and A. Chevarier, Journal of

	Nuclear Materials, 300, 27 (2002)
	[E. C. Subarao 1974] E. C. Subarao, H. S. Maiti and K. K. Srivastava, Physica Status Solidi.,
	9, (1974)
	[D. Wallace 1977] D. Wallace, "Carbon-14 Production in nuclear reactors", Oak Ridge
	National Laboratory report (1977)
	[M. C. Wittels 1956] M. C. Wittels and F. A. Sherrill, Journal Applied Physics, 27, 643 (1956)
	[M. C. Wittels 1959] M. C. Wittels and F. A.

  Table1and we used the good set for the GGA-

	PW91 calculations.
	atoms	Valence electrons in VASP	Notations in the work
	Zr	Zr_4d 2 5s 2 Zr_4s 2 4p 6 4d 2 5s 2	Zr_4 Zr_12
	O	O_2s 2 2p 4 O_2s 2 2p 4	O_h O_s
		C_2s 2 2p 2	C_h
	C	C_2s 2 2p 2	C_s
		C_2s 2 2p 2	C

Table 1 : The different pseudopotentials for atom Zr, O, and C from VASP library

 1

Table 1 : Cohesive energies and lattice parameters for a hcp Zr unit -cell, using Zr_4

 1 The results were listed in Table1.

		E cohesive	Parameter	Parameter
		(eV):	a, b (Å):	c (Å):
	VASP This work	6.73	3.23	5.14
	VASP [R. Li 2013]	6.74	3.23	5.17
	VASP [N. Poletz 2006]	6.36	3.25	5.15
	Experimental [C. Kittel 2005]	6.49	3.23	5.15

Table 2

 2

	.

Table 2 : Diamond carbon cohesive energy and lattice parameter

 2

	The zirconium carbide under ambient environment is a cubic close packed structure. Its
	space group is Fm3m, numbered 225, with parameters: a=4.698 Å, b=4.698 Å, c=4.698 Å, α=
	90°, β= 90°, γ= 90° [D. R. Lide 1998].

So it is provided that both C and C_s give accurate results comparing with experimental values. In this case, both can be used for this system.

A-3 ZrC structure

In the previous parts, it is found that the Zr_4, C and C_s pseudopotential are good candidates for pure Zr and C structure. In this part, to verify if the pseudopotentialscan accurately describe the Zr-C bonds, the ZrC structure was calculated.

  Table 3, which prove a good agreement with experiment data.

		E cohesive (eV):	Parameter a, b, c (Å):
	VASP This work	16.22	4.725
	VASP PW91[F. Vines 2005]	16.08	4.72
	Experiment [V. P. Zhukov 1985]	15.86±0.19	4.68

Table 3 : ZrC cohesive energy and lattice parameter

 3

	-77.3892				-77.5
						-77.55
	-77.7 -77.5446 Total Energy (eV)	Etot(7*7*7) Etot(9*9*9) Etot(11*11*11) Etot(13*13*13)	Total Energy (eV)	-77.7 -77.65 -77.6	Etot(7*7*7) Etot(9*9*9) Etot(11*11*11) Etot(13*13*13)
		300 400 500 600 700 800			300 400 500 600 700 800
		Cut-off Energy (eV)			Cut-off Energy (eV)
	Geometrical Parametre (Å)	4.7 4.71 4.72 4.73	Ecut=400 eV Ecut=450 eV Ecut=500 eV Ecut=550 eV Ecut=600 eV	Geometrical Parametre (Å)	4.73 4.7 4.71 4.72	Ecut=350 eV Ecut=400 eV Ecut=450 eV
		Kpoint				Kpoint

Table 4 : Volume variations before and after the insertion of one carbon in site 8 for different supercell sizes

 4

					3c
	2c		2c
	2b	2a	3b	3a	3b	3a
	3c		4c		5c
	4b	4a	4b	4a	5b	5a

Table 5 : Formation energy (in eV) for the different C interstitial positions

 5

Table 6 : Migration energies (in eV) of the possible migration pathway in zirconium bulk for Zr_4 and (Zr_12) in parenthesis

 6 To study the possible jumps of carbon through interstitial sites, the possible neighbor sites were identified. A carbon in BT_d site has 2 O sites and 2 BT_d sites as neighbor sites and for a carbon in O site, it has 6 neighbor BT_d sites and 2 neighbor O sites. So for a carbon in the O site, it can do 8 possible jumps, and for a carbon in the BT_d site, it can do 4 possible jumps. As summarize, there are four kinds of atomic jumps: O site BT_d site; BT_d site

	D-2 Attempt frequencies and jump rates

O site; O site O site; BT_d site BT_d site… In order to calculate the transition rate associated to each possible jump, the TST theory presented in Chapter 1 is used

Table 7 .

 7 The chosen displacement of the C atom at the stable sites and transition states in the 3 directions is 0.005 Å.

	final site
		O site	BT_d site
	initial site
	O site	1.3×10 13 s -1	1.7×10 13 s -1
	BT_d site	1.6×10 13 s -1	1.3×10 13 s -1

Table 7 :

 7 The attempt frequencies " "• of

the possible migration pathway in zirconium bulk for Zr_4

Table 8 : Different diffusion temperatures of carbon atom at different temperatures in the monoclinic zirconia

 8

	Temperature (K)	Diffusion coefficient (cm 2 /s)
	700	1.72×10 -17
	800	1.13×10 -15
	900	2.93×10 -14
	1000	3.97×10 -13
	1100	3.34×10 -12

Table 9 : Vacancy -carbon formation energy (E for ) and binding energy (E b ) for substitution and interstitial with a nearby vacancy (in eV)

 9

	Energy	Zr_4	Zr_4	Zr_12	Zr_12
	(eV)	E for	E b	E for	E b
	Subs	1.82	-3.50	1.90	-3.56
	BT_d	0.46	-2.16	0.41	-2.09
	O	-0.85	-0.86	-0.86	-0.83

  Monoclinic zirconia supercell construction ........................................................... I.1.D Investigation of carbon interstitials in the supercell of zirconia ........................... I.1.E Migration of C between two neighbor sites in the monoclinic zirconia without defects ............................................................................................................................... -4 DOS of 96 atomic supercell with a carbon substitution in O2 position with an O vacancy.................................................................................................................................................

	I.1.B Study of relative stability of the three phases of ZrO 2 (tetragonal, and cubic) ......
	I.1.C II.2.A DOS of the perfect 96 atomic supercell of monoclinic zirconia .........................
	II.2.B DOS of the 96 atomic supercell of monoclinic zirconia with defects .................

Introduction ............................................................................................................................ I Carbon behavior in perfect zirconia ................................................................................ I.1 Atomic scale simulations .............................................................................................. I.1.A Optimization of zirconia calculation parameters................................................... I.2 Macroscopic scale simulations ..................................................................................... II Carbon behavior in zirconia with point defects ............................................................ II.1 Vacancy and carbon substitution in the zirconia supercell ....................................... II.1.A Vacancy of O1, O2 and Zr in monoclinic zirconia ............................................... II.1.B Substitution of O and Zr by carbon ...................................................................... II.1.C Migration by vacancy mechanism ...................................................................... C-1 Self-diffusion of O by vacancy mechanism without carbon ............................................ 1) Migration to a vacancy O1 .............................................................................................. 2) Migration to a vacancy O2 .............................................................................................. C-2 Migration of C by vacancy mechanism .......................................................................... 1) C substitution in a neighbor O position of an O1 vacancy .............................................. 2) C substitution in a neighbor O position of an O2 vacancy .............................................. C-3 Migration of O with a neighbor carbon to a vacancy: competition effects .................... 1) Carbon in a neighbor O1 position of the O1 vacancy ..................................................... 2) Carbon in a neighbor O2 position of the O1 vacancy ..................................................... 3) Carbon in a neighbor O1 or O2 position of an O2 vacancy ............................................ C-4 Conclusion ...................................................................................................................... II.2 Electronic structure: Density of states ..................................................................... B-1 DOS of 96 atomic supercell with a carbon interstitial atom ......................................... B-2 DOS of 96 atomic supercell with an O2 vacancy ........................................................... B-3 DOS of 96 atomic supercell with a carbon substitution in O2 position ........................ BConclusions ........................................................................................................................... References .............................................................................................................................

Table 1 : Geometrical parameters of the

 1

	At 350 eV, the Zr with 4 e -
	the experimental data (Table 1
	calculation time.

: The

total energy of the zirconia unit-cell with different k-points and : Lattice parameter of the monoclinic structure with different k-

			Exp [Ref G. M.
	Zr_4 O_s	Zr_12 O_s	Rignanese 2005, C. J.
			Howard 1988
	1.01	1.01	1.01
	1.03	1.03	1.03
	5.20	5.21	5.15

parameters of the monoclinic zirconia for different pseudopotential sets orthe

Zr with 12 e -with O_s have similar results, comparing to 1). But Zr_12 which contains more valence electrons need more 83 points and E

cut-off for -points and E cut-off for

	Ref G. M.
	gnanese 2005, C. J.
	Howard 1988]
	1.01
	1.03
	5.15
	monoclinic zirconia for different pseudopotential sets
	results, comparing to

Table 2 : Geometrical parameters of the cubic

 2

	• Tetragonal ZrO 2

: The energy of a cubic structure unit-cell ZrO 2 with different k-points and Ecut parameter of cubic structure ZrO 2 with different k

		DFT [A. S.	Experimental
	calculated	Foster 2001, G.	[Ref C. J.
		Jomard 1998]	Howard 1988]
	5.12	5.09 / 5.15	5.09
	33.60	32.97 / 34.06	32.97

Table 3 .

 3 The calculated values are in agreement with the the literature.

	Cell parameters	calculated
	a (Å)	3.63
	c/a	1.44
	Dz	0.05
	Motif volume(Å 3 )	33.00

Table 3 : structural parameters of the tetragonal

 3

Table 4 : Cohesive energy (in eV) for the three phases of zirconia

 4 XtWoH is the atomic energy of Zr in hcpzirconium; • £ is the atom number of O in the supercell; • £ is the atomic energy of O in molecular O 2.

	geometrical parameters
	parameters were compared

XtWoH is the cohesive energy of supercell; is the atom numbers of Zr in the supercell; 86 1.414 before optimization [C. J. Howard , the convergence of the total energy was It was found that the

Table 5 : Four symmetrical positions of site 2 with their structures around carbon

 5

	Site	1	2	3	4	5	6	7	8
	N(site)	4	4	4	4	4	4	4	4
	E I (eV) 7.23 7.50 7.76 7.11 8.13 6.67 6.57 6.05

Table 6 : E I is the formation energy for the eight different interstitial sites of carbon in monoclinic zirconia. N is the number of equivalent site in the 12 atom supercell.

 6

Table 11 : The different C substitution energies sWt? in the bulk monoclinic zirconia

 11

Table 12 : The formation energy for different numbers of carbon in the Zr vacancy site

 12

Table 13 .

 13 There are also five neighbor O2 atoms:No.48, No.45, No.37, No.40 and No.39. The distances between the O1 vacancy position and neighbor O2 atoms were noted in green color in Table13.

	Possible jumps	34/66 (O1)	44/60 (O1)	41 (O1)	48 (O2)	45 (O2)	37 (O2)	40 (O2)	39 (O2)
	E mig (eV) i O1	1.29	1.80	2.55 2.31 2.26 2.11 2.64 2.35
	Distance(Å)	2.78	2.84	2.58 2.97 3.02 2.53 3.25 2.97

Table 13 : Migration energies for all the neighbor O atoms migrate to the O1 vacancy

 13 It is found that O2 atom No. 37 and O1 atom No.41 are the nearest O2 atom and the nearest O1 atom. Also for atom No.37, it needs the smallest migration energy of the five neighbor O2 atoms, but for the nearest O1 atom No.41, it needs the largest migration energy of the five neighbor O1 atoms. So there is no correlation for the migrations of the neighbor O atoms. It is concluded that for one O1 vacancy, the energy barrier of a neighbor O2 atom's migration is relevant about the distance between the vacancy and the neighbor atom. For the neighbor O1

Table 14 : Migration energies for all the neighbor O atoms migrate to the O2 vacancy

 14

	Possible jumps	35 (O1)	43 (O1)	66 (O1)	74 (O1)	60 (O1)	45/77 (O2)	78 (O2)	94 (O2)
	E mig (eV) i O2	2.38	2.71 2.33 2.18 2.42 2.74 2.34 3.04
	Distance(Å)	2.97	3.25 3.02 2.53 2.97 2.66 2.59 2.73

Table 15

 15

	).

Table 15 : The binding energies (in eV) between C in different neighbor O positions and the O1 vacancy and the migration energy of C with the corresponding distance(X: no direct migration path)

 15 It was found that for an O1 vacancy, the carbon in the pair of symmetrical positions has the smallest migration energy 0.83 eV. Comparing with the O migration in the C-1 part, it was also smaller than the migration energies of O self-diffusion in bulk without carbon. So for O1 vacancy, it should be calculated the migration of O with a carbon atom in a neighbor O position to verify the carbon effect on the migration energies of O, to understand what is competition between migration of O and migration of C to the same vacancy. To compare with the migration of O without carbon, still the O2 vacancy No.48 was considered with all its neighbor positions. The binding energies between carbon atomin different neighbor O positions and the O2 vacancy site were calculated and listed in Table16.The calculated binding energies arealways positive, which means that that the carbon is more stable in a neighbor O position of the O2 vacancy than in a substitution O2 position far away from the vacancy. The carbon atom prefers to go near the O2 vacancy. All the energy barriers of C atom in different positions migrate to the O2 vacancies were listed in Table16.

	2) C substitution in a neighbor O position of an O2 vacancy
	When the C is substitution in a neighbor O position of an O2 vacancy, there are 9 different
	O positions for C, since as seen previously, there are 9 neighbor O atom positions for an O2
	vacancy. sites	35 (O1)	43 (O1)	66 (O1)	74 (O1)	60 (O1)	45/77 (O2)	78 (O2)	94 (O2)
	E binding	2.34	3.44	2.54	3.34	3.13	3.45	3.62	3.30
	E mig (eV) i O2	0.39	X	0.78	0.17	0.48	1.59	1.09	1.62
	Distance(Å) 2.97	3.25	3.02	2.53	2.97	2.66	2.59	2.73

Table 16 :The binding energies (in eV) between C in different neighbor O positions and the O2 vacancy and the migration energies of C with corresponding distances(X: no direct migration path)

 16 It is found that,for the carbon in the nearest position of O1 atom No.74,it is much easier to do the migration since the energy barrier (0.17 eV) was smaller than that one of O self-

	diffusion in bulk without carbon calculated in C-1 part. So for O2 vacancy, as the same as for
	O1 vacancy, it should be calculated the O migration with a carbon substitution nearby to
	check what is the competition.
	As a conclusion for the migration of C by vacancy mechanism (C-1 and C-2), if a carbon
	atom diffuses by migration path of vacancy O2-O1-O2, the smallest energy barrier path for C

Table 19 : The compositions of different peaks in DOS of supercell with atom

 19

	There are two new electron p
	DOSCAR and PROCAR, the two
	and C_2p Zr_4d respectively.

total DOS of zirconia supercell with a carbon interstitial atom

	10.2	-7.3	-1.4	0
	0.856	0.826	0.819	0.838
	0.314	0.132	0.235	0.289
	0.207	0.198	0.097	0.056
	0.043	0.043	0.060	0.216
	: The compositions of different peaks in DOS of supercell with a carbon interstitial a carbon interstitial
	here are two new electron peaks centered at -1.4 eV and 0 eV. By analyzing th eV. By analyzing the
	DOSCAR and PROCAR, the two peaks were identified as C_2p from the inserted carbon s were identified as C_2p from the inserted carbon,
	and C_2p Zr_4d respectively.
				atom (derived integral
	DOS of 96 atomic supercell with an O2 vacancy
	and derived DOS of the supercell with an oxygen vacancy w and derived DOS of the supercell with an oxygen vacancy were presented in
	. They are perfectly matched which means the derived integral integral DOS can well
	present the density of states for the supercell.
				117

: The

  Table 20from the PROCAR file. The peaks at -0.15 eV and 0 eVare both

	composed by C_2p and Zr_4d.
	Peak (eV)	-0.15	0	0.06
	total	0.926	0.859	0.842
	C	0.333	0.333	0.327
	O	0.041	0.036	0.052
	Zr	0.134	0.148	0.126

Table 20 : The compositions of different peaks in DOS of supercell with a carbon atom insubstitution O2 position

 20

Table 21

 21

	. The peaks at -0.54 eV and
	Zr_4d. For the peak at Fermi level energy (0 eV), it is principally from the electrons of
	from the substituted carbon.
	Peak (eV)	-0.54
	total	0.874
	C	0.316
	O	0.022
	Zr	0.162

Table 21 : The compositions of different peaks in DOS of supercell with a carbon atom in O2 substitution position

 21

	supercell with a carbon substitution in O2 position with an O

compositions of different peaks in DOS of supercell with a carbon atom in O2 substitution position

			an O vacancy
	0.54 eV and -0.25 eV are both composed by electrons from C_2p and
	Zr_4d. For the peak at Fermi level energy (0 eV), it is principally from the electrons of
	0.54	-0.25	0
	0.874	0.839	0.759
	0.316	0.323	0.279
	0.022	0.039	0.033
	0.162	0.154	0.088
			121

: The

Table 22 : The binding energies and corrected binding energies between C in different neighbor O positions and the O2 vacancy

 22

  and Table 24.

	Sites	34/66 (O1)	44/60 (O1)	41 (O1)	48 (O2)	45 (O2)	37 (O2)	40 (O2)
	E binding (eV)	3.15	3.27	3.33	3.73	3.35	4.55	3.68
	E binding corrected (eV)	1.09	1.21	1.27	1.67	1.29	2.49	1.62

Table 22 : The binding energies and corrected binding energies between C in different neighbor O positions and the O1 vacancy

 22

	Sites	35 (O1)	43 (O1)	66 (O1)	74 (O1)	45 (O2)	77 (O2)	78 (O2)	94 (O2)
	E binding (eV)	2.34	3.44	2.54	3.34	3.45	3.45	3.62	3.30
	E binding corrected (eV)	0.28	1.38	0.48	1.28	1.39	1.39	1.56	1.24

Table 23 :

 23 The binding

energies and corrected binding energies between C in different neighbor O positions and the O2 vacancy

lp@p#@i = @+oo7q$Br lp@p#@i:

)/ -@+oo (8)

where:

• @+oo7q$Br lp@p#@i: is the total energy of the 95 atomic zirconium supercell with a Zr vacancy in it;

)/ is a zirconium atom's energy in the 96 atomic zirconium supercell; • @+oo is the total energy of the 96 atomic zirconium supercell.

The formation energy to form a zirconium vacancy is 1.87 eV. The value has a good agreement with previous theoretical ones 1.90 eV and 1.86 eV which were obtained by using DFT simulations [P. Gasca 2010, C. Domain 2005].

II.2. Substitution of Zr by carbon

In the case of substitution in zirconium, a carbon was placed in the position of a Zr vacancy.

The total energy after geometrical optimization was noted as E C-Zr(with substitution) .The formation energies of C in a Zr vacancy were calculated using the Eq. 9:

sWt = nX)/7q$Br sWtsB$BWB$Y#: -< n @+oo7q$Br lp@p#@i: >

where:

• nX)/7q$Br sWtsB$BWB$Y#: is the total energy of 96 atomic supercell with a Zr substitution by one carbon; • n is the reference carbon atom's energy in a diamond crystal; • @+oo7q$Br lp@p#@i: is the total energy of the 95 atomic zirconium supercell with a Zr vacancy.

The formation energy is 1.82 eV for a C substitution in a Zr vacancy position. So, the carbon atom is less stable in the position of a zirconium vacancy than in an interstitial position (by 3.52 eV relative to O site). A carbon atom always preferred to be in an interstitial position than to be in a substitution one.

II.3 Interstitial carbon with a nearest Zr vacancy

In this last part, the Zr vacancy effect on the formation energy of an interstitial C has been investigated using Eq. 10:

$#BXlp@p#@i = where:

• nX)/7q$Br lp@p#@i: is the total energy of 96 atomic supercell with a nearest Zr vacancy from an interstitial carbon; • n is a carbon atomic energy in a diamond crystal; • @+oo7q$Br lp@p#@i: is the total energy of the 95 atomic zirconium supercell with a Zr vacancy.

II.3 Interstitial carbon with a nearest Zr vacancy

In this last part, the Zr vacancy effect on the formation energy of an interstitial C has been = nX)/7q$Br lp@p#@i: -< n @+oo7q$Br lp@p#@i is the total energy of 96 atomic supercell with a nearest Zr vacancy from an interstitial carbon;

energy in a diamond crystal;

is the total energy of the 95 atomic zirconium supercell with a Zr All the migration barriers were calculated for all the possible atomic jumps and the results for energy barriers were listed in Table 7. For each possible atomic jump, the attempt frequencies were listed in For the formation energy results, the value is positive, which means the monoclinic zirconia is less stable when there is a vacancy. For Zr vacancy, the formation energy is larger than the one of an O vacancy, which provided a Zr vacancy less stable than an O vacancy. It is easier to form an O vacancy than a Zr one.

II.1.BSubstitution of O and Zr by carbon

Substitution calculationswere performed to verify whether a carbon prefers to stay in a vacancy or not. A carbon was positioned in the position of a vacancy O1, O2 or Zr. The total energy after geometrical optimization is noted as E C-ZrO2(with substitution) .The energies to form the substitution were calculated by the formula:

where:

• sWt? : the formation energy for the supercell to form a vacancy x, where x is the atom of vacancy; • nX)/£ 7q$Br sWtsB$BWB$Y#: : the total energy of the supercell with a carbon substitution; • n : the atom energy of carbon; • @+oo7q$Br lp@p#@i?: :the total energy of supercell with a vacancy x. Comparing the migration of carbon atom presented here with the migration of carbon atom by interstitial mechanism, the migration barriers are found to be comparable. Since there exits competition from the migration of neighbor O atoms, it makes the carbon migration more difficult.

3) Carbon in a neighbor O1 or O2 position of an O2 vacancy

The O2 vacancy No.48 was still used in this part of calculations. For carbon atom in one of the nine neighbor positions, the migrations of all the neighbor O atoms were performed. For carbon atom in neighbor O1 position, the different migration energies of neighbor O atoms were calculated and presented in Figure 22.

Conclusion

Consequently, when C will diffuse out of the zirconia, after several thousands or millions of years it will no more be active. Thus the release of 14 C should be expected to be really very small, and it should have a very limited contribution to the release of radionuclides outside.

This extremely low 14 C release should allow to improve the models of the container evolution in storage conditions, and contribute to reduce the cost of storage.

Perspectives

In the chapter 3, it is found that the charged O vacancy will change the binding energies between an O vacancy and neighbor C/O atoms. To see charge effect on the migration energies, the migration by charged vacancy mechanism can be simulated.

In this first step toward an understanding of the carbon behavior in the fuel claddings, it was proved that there is no carbon diffusion in bulk zirconia with and without vacancies. However, if there are some carbon atoms in the ground boundaries or cracks of the bulk, it should be interesting to simulate them.

Once the diffusion of carbon in ZrO 2 is well understood, it should be interesting to investigate the carbon form at the ZrO 2 / water interface. . Currently stored at La Hague, the waste will be stored in deep geological environment. The MA-VL waste contains RN including carbon-14, which comes from the neutron activation of nitrogen-14 and oxygen-17 present in the oxide Zircaloy. The objective of this thesis, which is taken in collaboration with EDF and AREVA, is to answer the question of the carbon-14 behavior in ZrO 2 and Zr metal. This thesis is to perform numerical simulations, using the Theory of Density Functional (DFT) to model the behavior of carbon-14 in the CSD -C. The simulations are performed with multi-scale approach: 1) At atomic scale, optimization of atomic models to represent the different phases of ZrO 2 ; identification of the insertion sites for the carbon atom in interstitial sites and substitutions of an oxygen or zirconium atom; modeling different minimum energy path for the migration of a carbon atom from one site to another by interstitial and vacancye mechanisms.

ANNEXE

2) At macroscopic scale, determination of diffusion coefficients in pure bulk monoclinic ZrO 2 and pure bulk Zr by the Monte Carlo method.