Keywords: dynamiques des dislocations, dynamiques moléculaires, plasticité cristalline, simulation dislocation dynamics, molecular dynamics, crystal plasticity, simulation

Ce travail vise à reproduire les interactions individuelles entre les dislocations vis et les boucles induites par l'irradiation en utilisant les simulations de dynamique des dislocations en accord avec les simulations de dynamique moléculaire. Un tel accord se caractérise par la reproduction de la réaction et avoir un valeur des contraintes critiques résolues pour franchir les obstacles. Cette approche fournit le moyen de calibrer notre code de dynamique des dislocations avec les paramètres des simulations de dynamique moléculaire. Par conséquent, il permet d'effectuer des simulations massives à l'échelle mésoscopique. Dans ce cadre, ce travail se compose de deux parties, une identification du modèle énergétique et une identification des mécanismes élémentaires. Dans la première partie, nous proposons une procédure de calibrage de la tension ligne basée sur le mécanisme d'Orowan en utilisant une étude de sensibilité. Dans la deuxième partie, nous avons identifié les mécanismes de glissement dévié et le maclage/antimaclage comme étant essentiels pour reproduire les interactions individuelles de dislocation-boucle. Les simulations de la dynamique des dislocations sont réalisées à l'aide d'un code nodal 3D appelé NUMODIS, où les développements récents dans ce code sont présentés. Un des caractéristiques de ce code est sa capacité à gérer et contrôler les collisions entre les segments des dislocations. Cela se fait au moyen en utilisant un ensemble d'algorithmes génériques avec un minimum de règles locales.

Twinning and anti-twinning directions are indicated for each plane . . . . .

List of tables

List of acronyms

General introduction

Industrial context

Most of the world's currently operating nuclear power plants have a life span of 30-40 years that could be possibly extended if safety conditions allow. During this period, several structural components suffer from a degradation in the mechanical properties. Such degradation is due to corrosion, mechanical stresses, neutron irradiation and other factors.

The reactor pressure vessel (RPV) is considered as the main critical structural component in nuclear power reactors, since it contains the core and it is connected to steam generators and other structural elements. In addition, the RPV serves as the second barrier against irradiation, the first being the fuel cladding, therefore it should maintain its integrity during normal conditions and in case of accidents [1]. A schematic diagram of a nuclear power plant and the pressure vessel is shown in Figure (1).

Figure 1 -Schematic diagram of the major components of a nuclear power plant of a pressurized water reactor. The RPV is situated inside the confinement, where fuel and control rods are indicated [2].

In pressurized water reactors (PWR), the RPV is designed to operate under a high pressure of 155 bar and a temperature around 300 • C [3]. During its lifetime, it is subjected to microstructural changes due to the presence of highly energetic fast neutrons. The microstructural damage occurs when such neutrons collide with atoms of the RPV. Atoms collided with neutrons are then displaced from their lattice site position in a relatively short time of few ps. This process leads to a continuous generation of radiation-induced defects such as self-interstitial atoms (SIA), clusters of solute atoms and voids [4]. Displacement per atom (dpa) is the unit used to describe the irradiation damage in the material. It represents the average number of lattice position change per atom and it was first proposed by Norgett et al. in [5].

Plastic deformation in crystalline materials, such as the steel used in RPV, is mainly related to the movement of dislocations. Strain hardening in metals can be attributed to the progressive introduction of barriers to the free movement of dislocations. The existence of radiation-induced defects also hinders the movement of dislocations. This leads to an increased hardness and embrittlement of the material [6]. One of the major design basis accidents in a nuclear reactor is the loss of coolant accident (LOCA) from the primary cooling circuit. When the DBT temperature decreases due to irradiation, in case of a LOCA for example, the temperature of the RPV might drop to values below the DBT temperature. It leads therefore to the initiation of a brittle fracture in the RPV. Hence, it is important to provide the end-users (operators) and decision makers with accurate models that describes the microstructural evolution of the RPV. Such models are used to predict the residual lifetime of the existing nuclear reactors to ensure safe operation in conformity with national and international regulations.

Irradiation effects on the RPV and internal steels was the core topic of successive European projects such as LONGLIFE, PERFECT and PERFORM 60. This study is part of a sequel project called SOTERIA. It stands for "Safe long-term operation of light water reactors based on improved understanding of radiation effects" [9]. This project is focused on the radiation-induced defects and microstructural heterogeneities using a combination of experimental and modeling work. In order to achieve the goals of SOTERIA, a materials multi-scale modeling (MMM), from atomic to macroscopic scale, is adopted. MMM is necessary to provide the means to understand the different phenomena that occur in the material due to irradiation.

Objectives

As a part of this MMM approach, the global objective of this work is to build in mesoscale simulations, namely dislocation dynamics (DD) simulations, as many inputs from atomistic simulations regarding the interaction between dislocations and radiation-induced defects. More specifically, the objective of this PhD work, is to reproduce the individual interactions of screw dislocations with different radiation-induced defects.

There are two types of dislocations in materials depending on the direction of their Burgers vector relative to the dislocation line, screw or edge dislocations, where a dislocation of mixed character can also exist [10]. In body-centered cubic (BCC) metals, yield stress is affected by the thermally activated motion of screw dislocations. The temperature above which the lattice resistance vanishes is identified as the athermal transition temperature (T a ) and it is found to be 340 K in α-Fe [11][12][13]. Because of this phenomenon, screw dislocations have a lower mobility at low and moderate temperatures and the deformation process is mainly controlled by their movement rather than edge dislocations.

In the case of interactions between dislocations and radiation-induced loops, molecular dynamics (MD) results show that the nature of the interactions does not depend only on the character of the dislocation (edge or screw), but also on the type of the defect, its size and the strain rate [14]. Interactions with screw dislocations have a high critical yield stress compared to those with edge dislocations. In this work, we concentrate on studying the interactions of screw dislocations with radiation-induced loops for the reasons mentioned before.

Irradiation defects are generated through rapid cascades in the RPV. These cascades take place in short length and time scales and they have consequences on the long-term at the macroscopic scale [15]. Displacement cascades have been extensively studied using atomistic simulation methods (ab-initio or molecular dynamics) [16]. Nevertheless, these methods lack the capacity to perform simulations at representative volumes relatively close to the macroscopic scale. Therefore, a materials multi-scale modeling is adopted to provide the link between different physical phenomena and processes at different time and space scales. A typical scheme of multi-scale modeling is shown in Figure (3). Materials multi-scale modeling aims to provide the link between different physical phenomena occurring at different time and space scales, which ultimately result in the macroscopic mechanical properties. The range of multi-scale modeling spans from atomistic to the macroscopic scale. More specifically, time scale spans from 10 -15 seconds to few minutes at the laboratory scale, while the space or length scale ranges from 10 -12 m 3 to few cm 3 . At each scale, there are different experimental techniques that provide an insight of the plasticity. The extracted information from each scale can be passed to the superior scale as input parameters.

Among the different simulation techniques is the ab-initio method1 . It provides a precise evaluation of the structure and energy of defects such as their formation and migration energies. Such values are not only used in some analytical models at the macroscopic scale, but also in Monte-Carlo simulations (cf. Fu et al. in [18]) and cluster dynamics (cf. Barouh et al. in [19]). In addition, ab-initio method helps in the development of semi-empirical interatomic potentials that can be used in molecular dynamics simulations.

At a larger scale of few nanometers, classical molecular dynamics provides an insight on the elementary deformation processes in materials by solving the equations of motion at the atomic scale2 . These calculations are based on semi-empirical interatomic potentials, such as EAM potentials (ex: [20][21][22]).

At the mesoscopic scale comes the dislocation dynamics, it represents an intermediate level linking the atomic scale to the level of continuum mechanics of mono/polycrystals [23]. It is based on linear elastic theory, in which dislocations are moving in an elastic continuum medium [24]. DD technique is a powerful tool to gain an insight into crystal plasticity through the complex evolution of dislocations. DD calculations are not only performed at the nanometric scale but also at the size of grains (few µm) to provide and enrich crystal plasticity laws. The latter describe the mechanical behavior of statistically representative population of dislocations based on their density only.

At a scale higher comes the finite-elements method (FEM). This method aims at solving the set of partial differential equations describing the continuum motion and/or equilibrium of an object. In this MMM approach, it is mainly used to simulate the time evolution of the microstructure at space and time scale in the range of few cm 3 and few seconds, respectively. Crystal plasticity laws are used in FEM. The results of FEM help in deriving polycrystalline plasticity laws that describe the deformation at the macroscopic scale [25].

This work was performed using a 3D nodal DD code developed at CEA, CNRS and INRIA, called NUMODIS which stands for NUmerical MOdeling of DISlocations [26]. The originality of this work lies in the attempt to reproduce the individual interactions of "screw dislocations" with radiation-induced defects in DD in a good agreement with prior MD studies. This approach ensures the validity of our results when performing massive simulations. The latter is used to predict plasticity at a larger scale using complex geometries of junctions and radiation-induced defects, showing the power of DD to handle such calculations.

Organization of manuscript

This thesis is organized in the following way:

Chapter 1 is a literature review of the most important topics related to the scope of this work. It consists of a general description of the steel used in French RPV concerning its chemical composition, metallurgical and mechanical properties. A description of radiationinduced defects and their effect on the mechanical properties of RPV steel is then made. The importance of screw dislocations in controlling the deformation mechanism in RPV steel is also underlined. Finally, results from available molecular dynamics studies are presented for later comparison with dislocation dynamics.

Chapter 2 is an overview of our methodology to compare dislocation dynamics to molecular dynamics and how to transfer information through different scales. This chapter provides a detailed description of the code NUMODIS used for this work, with particular focus on one of the most important features; the nodal split algorithm.

Chapter 3 is dedicated to the description of the modifications and models implemented in the code to reproduce molecular dynamics results using our dislocation dynamics code.

Chapter 4 represents the core of this thesis. In this chapter we present the work done to reproduce each individual interaction between a screw dislocation and radiation-induced loops. In particular, interaction with 1 0 0 square SIA loops of different orientations and 1/2[1 1 1] hexagonal SIA loops of different size. We present the details of each interaction and how to reproduce it correctly with DD in NUMODIS.

Steel in the nuclear industry

Introduction

This chapter is essentially a bibliographic review of the mechanical behavior of the French reactor pressure vessel steel. We mainly focus on the radiation-induced defects and study their influence on the mechanical properties of the steel at the microscopic scale.

In this chapter, the first section is a briefing on steels used in the nuclear industry. Then we focus on the low-alloy steel used in RPV, its chemical composition, metallurgical and microstructural properties. In section (1.3), we concentrate on the mechanical properties of the RPV steel and the difference between irradiated and unirradiated states at the macroscopic scale. The origin of the mechanical degradation of the RPV steel is then related to the presence of radiation-induced defects. We focus mainly on self-interstitial loops and precipitates for the sake of this study. Comparison of the interactions between radiation-induced defects and dislocations using either molecular dynamics or dislocation dynamics is finally presented.

Steel in the nuclear industry

Nuclear energy proved itself as a reliable and efficient technology that provides a stable source of energy. Many countries adopted the nuclear technology not only for electricity generation, but also for research purposes. The total capacity of nuclear energy is around 383 GW(e) and delivers around 2410 TW(e).h of electricity supply in 31 countries around the world [1]. There are several types of nuclear power reactors, in this study we concentrate on the pressurized water reactors. PWR is the most predominant design and it represents around two-thirds of world's nuclear reactors [1].

The selection process of nuclear reactors materials is of a specific importance. Due to irradiation effect, nuclear power plants are subjected to a unique operation environment that differs from any other application. Hence, certain properties are required for the selected materials to ensure a safe long term operation of the nuclear reactors [2]. These requirements are listed below:

• Mechanical properties; strength, ductility, toughness and structural integrity.

• Production; machinability, fabricability, compatibility, availability and cost.

• Thermal properties; heat transfer, thermal stability.

• Chemical properties; stability, corrosion resistance.

• Neutronic properties; irradiation stability. Different grades of steels are widely used in nuclear power plants. Steels are mostly used for structural materials such as the steam generators, piping, reactor pressure vessel and internals. They have a great resistance against internal and external stresses while maintaining the integrity of the structure. Regarding their thermal properties, they show a good stability against deformation at high temperatures. In terms of production, iron exists in immense amount and it has a competitive cost and availability compared to other metals. Steels manufacturers have wide experience with steel fabrication and machining. In addition to the above mentioned properties, the need for a material able to withstand the harsh irradiation environment made steels a viable option in the nuclear industry [2].

Table (1.1) shows the different types and grades of steel used in nuclear reactors. The choice of each type of steel depends mainly on the function of the component. Carbon steels are mainly used for components with low irradiation exposure, e.g. piping in the secondary circuit. The latter are connected to the RPV from the outer side, the irradiation dose is therefore low [2].

Components

Secondary piping

Vessel Pressurizer

RPV liner piping, Pump

The nuclear reactor coolant system transfers the generated heat in the reactor core to the turbines. Its main function is to maintain the temperature of the fuel within limits to avoid any meltdown. Light water is used as a coolant and a moderator in PWRs. Water circulates under a pressure of around 155 bars (15.5 MPa). The existing pressure prevents water from boiling before it reaches the steam generators.

Ferritic or low alloyed steels are commonly used to produce reactor pressure vessels and the pressurizer [5]. Although these components are less exposed to irradiation damage than the fuel, their mechanical behavior remains crucial during the lifetime of the nuclear reactor. One of the major problems that encounters nuclear engineers is the ductile-brittle transition phenomenon. Below this specific temperature, the material undergoes a brittle fracture instead of a ductile one [3]. This phenomenon might be detrimental in case of accidents such as loss of coolant accidents [5].

Austenitic stainless steels are a good choice for components in direct contact with water due to their corrosion resistance properties. Furthermore, they are mainly used in components of moderate and high exposure to irradiation e.g. primary piping, pumps and steel liner inside the RPV [2]. Generally, austenitic steels are like other face-centered cubic (FCC) crystals, they do not exhibit a ductile to brittle transition. Nevertheless, Mullner et al. reported that DBT is observed in some types of austenitic steels of specific compositions at certain temperatures [6].

Generalities on the RPV steel

In summary, low-alloy steels are used in the reactor pressure vessel for their mechanical properties, corrosion and irradiation resistance. In the next section, we concentrate on the RPV steel through a general overview of the bainitic structure.

Generalities on the RPV steel

In this section we focus on the low-alloy steel used in RPV. Mainly, its chemical composition and the heat treatment to obtain the desired microstructure. We also present metallurgical and microstructural observations.

RPV as a primary component in nuclear reactors

The RPV is considered as the primary component of the reactor coolant system. The RPV consists of a thick-walled cylinders stacked on top of each other and two hemispheres all welded together. One hemisphere is in the bottom and the other on the top. The latter contains openings for the control rod driving system. The RPV is coated inside with a stainless steel liner of small thickness for purposes of anti-corrosion protection. The vessel contains the reactor core which is composed of fuel assemblies, control rods and supports internal structures. It is shown as a cut-away in Figure (1.1). The size of the pressure vessel is inversely related to the power density of the reactor. For example, a pressurized water reactor has a smaller vessel than a boiling water reactor [7]. Several aspects must be taken into account in the design process of a PWR pressure vessel [2]:

• geometry, specify its height, external and internal diameter, top and bottom geometry, internal configuration of structural supports and the location of penetrations of inlets and outlets.

• operating conditions, including maximum temperature and pressure in steady-state and in case of accidents or transient-state.

• mechanical degradation due to irradiation conditions in terms of swelling and creep.

• ease of access for refueling, inspection and maintenance.

Chemical composition and heat treatment

The RPV is made with industrial low alloyed steel 16MND5 of bainitic structure. The French RPV is similar to the A508CL3 (forged) pressure vessels used in US reactors. The chemical composition of French RPV steel as specified by the French code of design and construction (RCC-M) is indicated in Table (1.2) [8]. Each alloying element has a specific effect on the microstructure and the final properties of the pressure vessel. As mentioned before, the RPV is composed of forged rings, flanges and nozzles welded together. The forging process of the initial ingot takes place in several steps at different temperatures between 1000 and 1300 • C [9]. In order to obtain the desired bainitic-ferritic microstructure, three distinct phases of heat treatments are applied to the initial forged pieces as follow [10,11]:

• Homogenization phase:

-Slow cooling in the furnace after the forging process until 350 • C.

-Austenization at 900 -950 • C during at least 6 hours.

-Quenching in water down to 350 • C during 20 -42 minutes.

-Annealing at 650 • C during at least 6 hours for stress relief purposes.

-Cooling inside or outside the furnace.

• Enhancement of the mechanical properties:

-Austenization at 865 -890 • C during 3 -6 hours, to obtain a fine grain size.

-Quenching in water during 58 -62 minutes, this gives the desired bainitic structure.

-Annealing at 635 -655 • C during 5 hours and half to remove internal stresses.

-Cooling inside or outside the furnace.

• Tempering phase:

by reheating the vessel up to 615 • C and maintaining the temperature for 10 hours followed by slow cooling at 20 • C/h. This phase is necessary to increase the toughness of the steel and for stress relief in the welding zones.

The initial temperature of the ingot and the cooling rate during quenching from the austenitic phase is the main factor to obtain the bainitic-ferritic microstructure. The resulting microstructure is described using the continuous cooling transformation diagram (CCT) as in Figure (1.3). A unique CCT diagram exists for each alloy. The bainitic microstructure is formed using a moderate cooling rate, slower than that used to form martensite and faster than that to form pearlite. Bainite is formed of two distinct categories, upper and lower bainite [3]. Since the RPV has a considerable thickness of more than 20 cm, the inner and outer sides are subjected to different cooling rates. This has an an effect on the resulting microstructure of the RPV, which is discussed in the following section.

Microstructure of the RPV

The microstructure of the RPV steel is not homogeneous along its thickness. Since the RPV has a thick-walled geometry, it is subjected to different cooling rates at the inner and outer surfaces. The outer shell cools down at a rate of ∼ 1-2 • C/s which leads to a totally bainitic microstructure. On the other side, the core cools down at a slow rate of ∼ 0.1-0.2 • C/s. The bainitic structure, in general, consists of a dislocation-rich ferrite with carbides grouped in carbon-enriched clusters. Upper bainite is comprised of lath bundles, while lower bainite is in the form of individual plates. [13]. Observations of the bainitic microstructure of 16MND5 RPV steel is done by optical microscope and scanning electron microscope (SEM). Observations using these techniques are shown in Figure (1.4) [14]. As can be seen from Figure (1.4), the microstructure of bainitic steels looks like the martensitic steel. It is therefore hard to be distinguished only by images, even when high resolution microscopy techniques are used. The similarities in the morphology are due to the existence of laths of different orientations.

During the previously described heat treatment of the RPV steel, it is possible to form carbon rich zones in form of cementite (Fe x C), where x ranges between 2-3. More precisely, it is formed by a chemical segregation during the solidification process [15]. Observation of the carbides in the RPV steel using SEM in the inter-granular and intra-granular sites is shown in Figure (1.5). Carbide segregation phases contribute to the degradation of the mechanical properties of the RPV steels, because they are considered as susceptible zones for crack initiation [17].

In conclusion, we presented the bainitic structure of the RPV steel and how to obtain such microstructure. In the next section, we present the mechanical properties of the bainitic steel at the macroscopic scale.

Mechanical behavior of RPV

In this section, we present the mechanical properties of the RPV. We also show the difference between the irradiated and unirradiated states as reported in literature.

The mechanical behavior of 16MND5 RPV steels has been widely studied in the literature at microscopic and macroscopic scales [18][19][20]. The RPV steel has a bodycentered cubic crystallographic structure, the stress-strain relation is therefore affected by temperature. Such relation is shown in Figure (1.6) for 16MND5 RPV steel under a strain rate of 5 × 10 -4 s -1 [18]. As can be seen in Figure (1.6), the yield stress increases when temperature decreases. Furthermore, the fracture stress decreases when temperature increases due to the ductile to brittle transition. The transition temperature between the thermal and athermal regimes was reported to be around 300 K at a strain rate of 10 -4 s -1 . DBT temperature is closely related to the strain rate [15,18] Embrittlement due to irradiation can be characterized using the Charpy test. In French nuclear reactors, the reference temperature for DBT temperature is taken when the average rupture energy equals 56 joules [21]. As mentioned before, a degradation in the mechanical properties occurs due neutron-induced embrittlement. The yield stress is observed to increase when the irradiation dose increases. The increase in the yield stress leads to a substantial hardening of the material. Stress-strain curves can be fitted using Ramberg-Osgood correlation [24]:

ε = σ E + 0.002( σ σ 0.2 ) n (1.1)
where E is the Young's modulus, σ 0.2 is the equivalent yield stress at 0.2% elongation and n is the hardening exponent. The latter is a fitted parameter that describes the knee sharpness in the stress-strain curve with respect to irradiation.

In summary, the ductile to brittle transition is observed in RPV steel. The DBT temperature increases with irradiation, the material therefore might be fragile at high temperatures. This could be detrimental in case of accidents in nuclear reactors. In sections (1.4, 1.5), we present a review of dislocations, their relation to mechanical deformation and the origin of the irradiation effect on the mechanical properties at microscopic scale.

Dislocations and plastic deformation

This PhD thesis, among many other dissertations written by students such as [19,20,25,26], discusses the role of dislocations in the mechanical deformation of materials. Moreover, basic information about the theory of dislocations can be found in general textbooks such as [3,[27][28][29]. We dedicate this section to remind of the most important basic concepts of dislocations that are discussed in the following chapters.

1.4. Dislocations and plastic deformation

Dislocations in general

In the early 20s of the bygone century, the difference between theoretical and experimental values of applied shear stress to plastically deform a single crystal was explained by the existence of dislocations. Plastic deformation occurs by atomic planes sliding over each other when shear stress is applied. A dislocation is simply a line defect in crystalline material that allows for such plastic shear without moving all the atoms in the shear planes. There are two main types of dislocations; edge and screw. An edge dislocation is formed by the insertion of an extra half-plane of atoms in the crystal. It has a Burgers vector normal to the dislocation line. A screw dislocation is created by displacing part of the crystal in opposite directions relative to each other. It has a Burgers vector parallel to the dislocation line. In addition, a dislocation of mixed character can also exist. Such dislocation is by definition a partial superposition of dislocation with edge and screw character. The formation process of edge or screw dislocations is explained in Figure (1.10).

Cross-slip mechanism

The cross-slip is known phenomenon for screw dislocations where the dislocation may deviate in a cross-slip plane. The dislocation Burgers vector is along the intersection of the primary and the cross slip planes. A simple illustration of such process is depicted in Figure (1.11).

Twinning and anti-twinning directions

Experiments performed on single crystals made of pure iron by Spitzig and Keh in 1970 [30] and Molybdenum by Him, Lau and Dorn in 1972 [31] showed asymmetrical deformation processes around {1 1 2} planes in BCC metals. This feature, called twinning/anti-twinning, can be related to the atomic asymmetry around these planes as illustrated by Figure (1.12) (e.g. Chaussidon 2007). If we imagine cutting the crystal in two parts along the (1 1 2) plane, this asymmetry makes shifting the upper-part with respect to the lower part along the b vector more difficult that in the opposite direction. This argument can be extended to the motion of 1 1 1 screw dislocations in {1 1 2} planes, which is clearly asymmetrical along 1 1 0 directions. Figure 1.12 -Illustration of the origin of twinning / anti-twinning directions for 1 1 1 screw dislocations: (left) Atomic structure of BCC as seen from a [1 1 1] direction. Along this direction, the BCC structure can be seen as the superposition of three atom strings as highlighted by three colors (green, blue, red), shifted by b/3 along the [1 1 1] direction. (right) Atomic structure as seen from the side. The atomic arrangement is not symmetric {1 1 2} planes (after Chaussidon 2007) [32].

Figure (1.13) shows the possible twinning and anti-twinning {1 1 2} planes associated to a [1 1 1] screw dislocation with both a line vector and Burgers vector pointing outwards. This phenomenon has also a direct influence on the ability of screw dislocations to cross-slip between {1 1 0} planes (dashed green line in Figure 1.13) under an applied resolved shear stress (e.g. Chaussidon 2007). In a nutshell, these screw dislocations will effectively avoid cross-slipping in the [-30 • ; +30 • ] angular sectors around {1 1 2} anti-twinning planes, while cross-slip will be favored between neighboring {1 1 0} planes if the maximum resolved shear stress direction lies in a twinning angular sector.

Radiation-induced defects

In this section, we present the different types of radiation-induced defects. We mainly focus on the self-interstitial atom loops and precipitates, because of their high density and obstacle strength. Moreover, we present different molecular dynamics studies of the interaction between dislocations and radiation-induced defects. A comparison between MD and DD is also done for some interactions.

Microstructural defects are continuously generated in RPV through rapid cascades (few ps) at a fine scale (few nm) when subjected to fast neutron flux [33]. When an energetic particle collides with an atom, the latter might be knocked out from its lattice position. The first displaced atom is then called a primary knocked atom (PKA). The PKA might travel to an interstitial site creating a vacancy behind. If the PKA has a high energy, it might launch a cascade and create a large agglomeration of knocked atoms. Such nuclear reactions between lattice atoms and neutron irradiation occur in the RPV steel. After certain reactions, atoms transmutation can also take place which produce a new metal atom, or gas such as Hydrogen or Helium [33].

Radiation damage is generally described in terms of displacements per atom from their lattice positions as proposed by Norgett et al. [34]. Atoms of the RPV in light-water reactors (LWR) normally has a displacement rate of 0.03 dpa/yr. In addition, dpa in research reactors is much higher than that for LWR [4]. Generally, a small fraction of the displaced atoms actually induce defects. This is because most of the displaced atoms recombine with existing vacancies through an annihilation process. Atoms that do not recombine, migrate to sinks such as grain boundaries, dislocations and existing cavities [35].

There are several types of radiation-induced defects in nuclear reactors, such as Cu-rich precipitates, cavities, solute cluster (e.g. Mn, Cr, etc.) and loops. Their density and volume fractions as function of neutron dose were reported in Meslin et al. In general, transmission electron microscopy (TEM) is used to observe the radiationinduced loops and dislocations [28]. Other techniques can be used to characterize different radiation-induced defects. Atom probe tomography (APT), scanning electron microscopy and positron annihilation spectroscopy (PAS) are examples of such characterization techniques.

The characterization of radiation-induced defects provides information on the size and concentration. In the following, we present the effect of each type of defects on the mechanical properties. Lambrecht et al. in [37] investigated the size and concentration of radiation-induced defects in pure iron and in different alloys including the RPV alloy. The contribution of each type of radiation-induced defects to the hardening process is shown in Figures (1.15, 1.16). In case of low irradiation dose, radiation-induced loops are not observed in RPV steel. However, Lambrecht et al. suggested that there are invisible loops in the RPV steel. Such invisible loops have a high obstacle strength factor (as presented in Equation 1.2) and contribute to the hardening process following the equation [37].

∆σ = αM µb √ N.d (1.2)
where α represents the obstacle strength, M is Taylor factor = 3.06, µ is the shear modulus = 71.8 GPa and b is the Burgers vector = 0.249 nm. N is the number density of radiation defects and d is the mean size. The values of N and d are obtained from the experimental methods.

The contribution of each defect to the hardening process is a function of the obstacle strength and their density as shown in Figure (1.16). SIA loops have a high obstacle coefficient α in Equation (1.2) compared to other defects. Although SIA loops have a high obstacle strength, their contribution to hardening is lower than the precipitates. This is because the latter have a much higher density than the SIA loops. Radiation-induced defects have a substantial influence on the mechanical properties of the RPV steel [38,39]. For example, embrittlement and fracture, hardening, swelling and creep. At the operation temperature of a RPV, the embrittlement is attributed to the presence of solute cluster and point defects cluster. The latter consists of dislocation loops or cavities [38,[40][41][42].

The four types of radiation-induced defects are discussed in the following.

Copper-rich precipitates

Although Copper is present only in minor concentrations in the composition of the RPV steel, it tends to form nanometric defects under the influence of irradiation. Other alloying elements form precipitates as well e.g. Ni, Mn, Si and Cr. However, Buswell et al. show that Cu precipitates are found to be pure [43]. Other experimental (cf. Pareige et al. [44]) and modeling (cf. Odette et al. [45]) studies show that other alloying elements are concentrated at the Cu-Fe matrix interface. The morphology of such precipitates is revealed using atom probe tomography. Figure (1.17) shows the pure precipitates in irradiated binary iron-copper alloy (Fe-0.7%Cu, Fe-1.4%Cu) at 288 • C. Odette et al. reported that Cu precipitates have a size of about 1.5 nm and a concentration of around 10 24 /m 3 when steel is irradiated at 10 19 n.cm -2 and 290 • C [45]. The crystallographic structure depends on the size of the precipitate. The precipitates might undergo a martensitic transformation from BCC to FCC when the size increases [46]. The previous finding is coherent with prior studies obtained with atomistic simulations (cf. Bacon et Osetsky in [47]). The martensitic transformation is also reported by Harry et Bacon, when a screw dislocation interacts with Cu precipitates [48].

Copper precipitates are considered as a pinning source to the dislocation motion. Bacon and Osetsky 2003 and 2004 studied the contribution of Cu precipitates to the mechanical strengthening at the atomic level. The mechanism of dislocation-precipitate interaction depends on the size of the Cu precipitate. For small diameters, no climb or martensitic transformation is observed. For large diameters, the observed reactions are non-conservative and the production of vacancy and interstitial point defects is observed. Such interaction is found to cause the precipitate to change its crystallographic structure and lead to possible climb of the dislocation [47,49]. The relation between the critical shear stress (CSS) and the size of the precipitates and temperature is given in Figure

(1.18a) (a) (b)

Cavities

Vacancies are known to migrate to sinks in the crystal at high temperature and irradiation dose. They agglomerate in a three dimensional manner and form cavities of large geometry [50][START_REF] Eldrup | Dose dependence of defect accumulation in neutron irradiated copper and iron[END_REF][START_REF] Stoller | Primary damage formation in bcc iron[END_REF]. Cavities can be divided into two types, bubbles and voids, depending on gas content and pressure. Bubbles are Helium gas atoms generated by (n,α) reaction with boron atoms in the steel. Commonly, the nanometric cavities have a nearly spherical morphology. It is suggested that their shape is influenced by the presence of alloying elements [START_REF] Meslin | Mécanismes de fragilisation sous irradiation aux neutrons d'alliages modèles ferritiques et d'un acier de cuve : amas de défauts[END_REF][START_REF] Miller | Embrittlement of RPV steels: An atom probe tomography perspective[END_REF].

Cavities are generally considered as strong obstacles to the motion of dislocations. Thus, they have a higher contribution in the hardening process [40,[START_REF] Terentyev | Interaction of an edge dislocation with voids in α-iron modelled with different interatomic potentials[END_REF]. Several studied of cavities using MD and DD simulations are found in literature [START_REF] Monnet | Mechanical and energetical analysis of molecular dynamics simulations of dislocation-defect interactions[END_REF][START_REF] Monnet | Mesoscale thermodynamic analysis of atomic-scale dislocation-obstacle interactions simulated by molecular dynamics[END_REF][START_REF] Queyreau | Modeling the dislocation-void interaction in a dislocation dynamics simulation[END_REF]. Osetsky et al. studied the effect of temperature on the critical shear stress of an edge dislocation to overcome a 2 nm cavity. The latter study was performed using atomistic simulations. The result is shown in Figure (1.19). As the temperature increases, the critical shear stress continuously decreases. Terentyev et al. in [START_REF] Terentyev | Interaction of an edge dislocation with voids in α-iron modelled with different interatomic potentials[END_REF] also studied the temperature effect on CSS of an edge dislocation to overcome a periodic array of voids in α-iron using molecular dynamics. The configuration is shown in Figure (1.20). Two different interatomic potentials [START_REF] Ackland | Computer simulation of point defect properties in dilute Fe-Cu alloy using a many-body interatomic potential[END_REF][START_REF] Ackland | Development of an interatomic potential for phosphorus impurities in α-iron[END_REF] are used for comparison. In a static simulation study at 0 K, Osetsky and Bacon investigated the critical shear stress dependence on void's size and regular mean spacing as shown in Figure (1.22). It is clear that the critical shear stress is directly related to the size of the void. However, it is inversely proportional to the center-to-center spacing for a periodic row of voids in a crystal. Such dependence is formally predicted with Equation (1.3) [START_REF] Bacon | The effect of dislocation selfinteraction on the orowan stress[END_REF]:

τ c = Gb 2πL [ln(D -1 + L -1 ) -1 + B] (1.3)

Solute Cluster

Solute atoms might be alloying elements intentionally added to adjust certain properties of materials. They can also be impurities from the elaboration process, or they can be created through nuclear reactions during the lifetime of the RPV. Typical solute atoms in the RPV steel are Cu, Ni, Mn, P and Si [START_REF] Jumel | Interactions between dislocations and irradiation-induced defects in light water reactor pressure vessel steels[END_REF]. They are characterized using atom probe tomography or by small angle neutron scattering (SANS). These techniques revealed the existence of clusters smaller than 1.5 nm with a density of around 10 24 m -3 when steel is irradiated at 10 19 n.cm -2 neutron fluence and 290 • C [START_REF] Akamatsu | Experimental evidence of several contributions to the radiation damage in ferritic alloys[END_REF]. Solute cluster are known to have a dilute morphology compared to copper-rich precipitates, as they might contain up to 85% Fe atoms [START_REF] Jumel | Interactions between dislocations and irradiation-induced defects in light water reactor pressure vessel steels[END_REF].

Radiation-induced segregation is affected by the presence of point defects since they might act as sinks trapping solute atoms or cause them to repel. An extensive description of the formation process of solute clusters and their microstructural evolution under irradiation is reported in Vincent 2006 [START_REF] Vincent | Simulations numériques à l'échelle atomique de l'évolution microstructurale sous irradiation d'alliages ferritiques[END_REF].

Becquart et Domain reported that the interaction between solute atoms and point defect (SIAs or vacancies) in iron has elastic and chemical parts. The size and the electronic structure of solute clusters determine the predominant effect. The elastic effect is found to be predominant for large clusters. The chemical effect is more important in case of different electronic structures between solute clusters and the matrix [42].

Pascuet et al. studied the segregation of solute atoms around an edge dislocation dipole using MD. They also investigated the effect of solute atoms on the pinning of dislocation lines [START_REF] Pascuet | Solute effects on edge dislocation pinning in complex alpha-Fe alloys[END_REF]. A visualization of such segregation is shown in Figure (1.23). The segregation of solute atoms at the dislocation core plays an important role in the mobility of dislocations. The latter is reduced due to the presence of solute clusters that act as strong obstacles for edge dislocations [START_REF] Pascuet | Solute effects on edge dislocation pinning in complex alpha-Fe alloys[END_REF].

Radiation-induced loops

There are two types of dislocation loops based on their formation mechanism. Loops are formed either due to a collapse of atomic planes over an agglomeration of vacancies, or by a set of self-interstitial atoms. Only the latter are observed in the RPV and they are discussed in this section.

Prismatic dislocation loops are characterized by their Burgers vector, which is perpendicular to their plane. They can only glide in the direction of their Burgers vector. Therefore, they cannot expand or shrink unless a thermally activated climb occurs, and by an absorption process of other SIAs created by subsequent damage [28,[START_REF] Olander | Fundamental aspects of nuclear reactor fuel elements[END_REF].

As mentioned in the previous section, TEM is used to observe dislocation loops in irradiated materials. Zinkle et Singh showed that dislocation loops are only visible using TEM at high irradiation dose ∼ 0.2 dpa [50]. Figure (1.24) gives the microstructural evolution of irradiated steel at different doses. The majority of observed dislocation loops are interstitials and have a Burgers vector equal to 1/2[1 1 1] in {1 1 1} planes. MD simulations show that the lowest energy configuration of dislocation loops is for Burgers vector 1/2 1 1 1 with segments in 1 1 2 directions [START_REF] Bacon | Modelling dislocation-obstacle interactions in metals exposed to an irradiation environment[END_REF]. Other studies showed the existence of square dislocation loops of Burgers vector equal to 1 0 0 in the 1 0 0 or 1 1 0 directions [50,[START_REF] Horton | A TEM study of neutron-irradiated iron[END_REF][START_REF] Robertson | Low-dose neutron-irradiation damage in α-iron[END_REF][START_REF] Nicol | Matrix Damage in Iron[END_REF]. This configuration is most observed at high temperatures [START_REF] Dudarev | Effect of the α -γ phase transition on the stability of dislocation loops in bcc iron[END_REF]. The evolution of density and size with radiation damage in dpa is depicted in Figure (1.25a) [50]. Molecular dynamics is used to investigate the effect of loop's size and orientation on the interaction mechanism with screw of edge dislocations. In the following, three main MD studies are discussed. The first study is focused on edge dislocations, while the second and third studies are focused on screw dislocations.

All of these studies used an interatomic potential for pure iron by Ackland et al. [START_REF] Ackland | Development of an interatomic potential for phosphorus impurities in α-iron[END_REF]. This interatomic potential is fitted to ab-initio simulations using the embedded atom method (EAM). There are many other potentials used in the literature depending on the chemical composition of the studied material [START_REF] Ackland | Computer simulation of point defect properties in dilute Fe-Cu alloy using a many-body interatomic potential[END_REF][START_REF] Henriksson | Atomistic simulations of stainless steels: a many-body potential for the fe-cr-c system[END_REF][START_REF] Mendelev | Development of new interatomic potentials appropriate for crystalline and liquid iron[END_REF]. However, the interatomic potential of Ackland et al. provides a better prediction of the dislocation core structure in good agreement with ab-initio simulations [41,[START_REF] Domain | Simulation of screw dislocation motion in Iron by molecular dynamics simulations[END_REF][START_REF] Chaussidon | The glide of screw dislocations in bcc Fe: atomistic static and dynamic simulations[END_REF].

Interaction between edge dislocations and SIA loops

There are several factors that determine the nature of the interaction between a dislocation and a loop, such as the strain rate, the character of the dislocation (screw or edge), size and orientation of the loop [START_REF] Bacon | Dislocation-obstacle interactions at atomic level in irradiated metals[END_REF]. Interactions between edge dislocations and radiation-induced defects are widely studied in literature compared to that of screw dislocations. A summary of such interactions in iron are found in Granberg 2016 and in Bacon et al. 2009 [38, 78].

A study by Bacon et al. shows the interaction between an edge dislocation with a hexagonal loop with b =1/2[1 1 1] and 1/2[1 1 1] of 37 and 331 SIAs. The EAM interatomic potential used in this study is Ackland et al. in [START_REF] Ackland | Development of an interatomic potential for phosphorus impurities in α-iron[END_REF]. They performed static and dynamic simulations at T = 0 K and at different temperatures in the range between 100-450 K. The applied strain rate in this study is 10 7 s -1 . The initial configuration is shown in Figure (1.26) [START_REF] Bacon | Computer simulation of reactions between an edge dislocation and glissile self-interstitial clusters in iron[END_REF]. The primary results of static simulations at T = 0 K shows that the loop of 37 SIAs is transformed into a superjog of the same Burgers vector as the dislocation. The formation of the superjog in this case is done through the change of the crowdion axis of the loop from [1 1 1] to the [1 1 1] direction. A visualization of such interaction mechanism is also shown in Figure (1.27). In this case, the dislocation-superjog combination is glissile and can glide along the plane of the original dislocation.

Radiation-induced defects

The same mechanism of superjog formation is observed when dynamic simulations are performed at 300 K with a large loop of 331 SIAs. The dislocation is first attracted to the loop, then a [0 1 0] segment is formed. The formed [0 1 0] segment is sessile. The dislocation is therefore blocked and then glides and bows in the [1 1 1] direction. With the increased stress level, the sessile segment is transformed into a [1 1 1] segment. After this reaction, a superjog is immediately formed. The overall process is illustrated in Figure (1.28). A summary of the critical stress for different temperatures for small and large loops of 37 and 331 SIA, respectively, is shown in Figure (1.29). Two main results can be noted from this figure; i) temperature has nearly no effect on the critical stress for small loops of 37 SIAs, ii) for large loops of 331 SIAs, the critical stress decreases when temperature increases, iii) the critical stress of large loops is higher than that for small loops regardless of the temperature. Comparison between the critical stress to overcome the obstacle calculated with the two simulation methods is shown in Figure (1.34). The difference between the critical stress found in MD and DD in the three cases was 14%, 18% and 19% for C2, C4U and C5 cases, respectively. Although there is a gap between MD and DD results, it is completely justified due to the computational accuracy and approximations in the physical models in DD.

Interaction between screw dislocations and SIA loops

Screw dislocations in BCC metals are known to have a thermally activated glide controlled by the double-kink mechanism [START_REF] Domain | Simulation of screw dislocation motion in Iron by molecular dynamics simulations[END_REF][START_REF] Chaussidon | The glide of screw dislocations in bcc Fe: atomistic static and dynamic simulations[END_REF]. The difference in the core structure between screw and edge dislocations mainly leads to different reaction mechanisms.

Terentyev et al. studied the interaction between a screw dislocation and 1 0 0 radiation-induced loops of 128 SIAs in α-iron. The authors performed dynamic simulations at different temperatures of 100 K, 300 K and 600 K. The applied strain rate equals 1 × 10 7 s -1 that corresponds to a steady state velocity of 27 m/s. The initial configuration, the simulation box and the different orientations of SIA loops are depicted in Figure (1.35). The interactions at 300 K (ambient temperature) are discussed below, notably A3 and B3. In the case of A3, the loop has a Burgers vector of [0 0 1] and after the interaction, the loop was restored with the same initial Burgers vector. First, the dislocation cross-slips to meet the lower part of the loop at point C. Dislocation arms then convert segments BC and CD to 1/2[1 1 1] by the following reaction in Equation (1.4).

[0 0 1] -1/2[1 1 1] = 1/2[1 1 1] (1.4)
When dislocation arms bow out, the two sides of the mobile dislocation line recombine. The segments of the loop are then restored as [0 0 1] when the dislocation breaks away. In the case of B3, the original loop is absorbed as a helical turn. The first interaction is reported to be repulsive. As the level of stress is increased, the dislocation cross-slips and encounters the SIA loop at point D. Segment AD is transformed to 1/2[1 1 1] through the following reaction.

1/2[1 1 1] + 1/2[1 1 1] = [0 1 0] (1.5)
The formation of the helical turn then takes place. The dislocation is liberated as in Hirsch's mechanism, similar to the interaction mechanism with loop of 37 SIAs in Liu et Biner in [START_REF] Liu | Molecular dynamics simulations of the interactions between screw dislocations and self-interstitial clusters in body-centered cubic Fe[END_REF]. A summary of the normalized critical stress of each reaction in [START_REF] Terentyev | Reactions between a 1/2<111> screw dislocation and <100> interstitial dislocation loops in α-iron modelled at atomic scale[END_REF] is shown in Figure (1.39). As can be seen in this figure, the critical stress of the mechanisms that result in a helical turn is always higher than the other mechanisms. Reactions that restore the original loop have the lowest critical stress. In another study, Liu et Biner studied the interaction of screw dislocation with 1/2 1 1 1 hexagonal SIA loops of different sizes in α-iron. The diameter of each loop was determined by the number of SIAs contained in the loop. The latter quantity is given by a(2N ) 1/2 /3. The studied loops had 37, 127 and 271 SIAs which is equivalent of cluster size of 1.6, 3.0 and 4.4 nm respectively. The length of the screw dislocation was equal to 30 nm extending along the length of the simulation box. Periodic boundary conditions were applied along the directions of the dislocation line and free surfaces in the glide direction. Furthermore, free boundaries were applied on the upper and lower surfaces. These simulations were performed at 300 K and at a strain rate 8 × 10 6 s -1 that corresponds to a steady state velocity of 30 m/s [START_REF] Liu | Molecular dynamics simulations of the interactions between screw dislocations and self-interstitial clusters in body-centered cubic Fe[END_REF].

Two interaction mechanisms are observed for this configuration depend on the size of the SIA loop. For the smaller loop of 37 SIAs a helical turn is observed (see Figure 1.40). At the beginning of the simulation, the screw dislocation glides in the 1 1 0 plane toward the loop. The loop then moves along its crowdion axis below the glide plane then moves again toward the dislocation under the effect of the applied stress field and the attraction by the dislocation stress field. A sessile junction of type [0 0 1] is formed when the dislocation meets the loop. The sessile junction blocks the movement of the dislocation. Under increased stress, the junction is dissociated and converted to a 1/2[1 1 1] segment. The rest of the loop is transformed into a helical turn through successive steps of cross-slip during a relatively short time. The helical turn has six segments. It can only glide along the direction of its Burgers vector which is parallel to the direction of the dislocation. This mechanism was previously described by Hirsch and therefore it is often called as Hirsch's mechanism [START_REF] Hirsch | Point defect cluster hardening[END_REF]. In this mechanism, in order for the dislocation to be liberated it has to close the helical turn. It leaves a loop of the same Burgers vector as the initial screw dislocation, in this case it is 1/2[1 1 1].

Regarding loops of large size, a planarization mechanism is observed for loops of 127 and 271 SIAs. In this mechanism, the formation of the first junction is similar to that of 37 SIAs. However, the formation of the helical turn is not favorable due to the large size of the loop. The MD results of this interaction are shown in Figure (1.41). The critical stress for the three cases was measured at the junction breakup shear stress minus the flow stress needed to move the dislocation alone. A summary of these results is shown in Table (1.4). The strength coefficient is calculated using Orowan's relationship as ∆σ = αµb/L, where µ is the 1 1 1 shear modulus (71 GPa), b is the length of Burgers vector. L is the difference between the length of the straight screw dislocation and the loop diameter d. The values of α ranged between 0.4-0.6 for 300 K. As can be seen from Table (1.4), the formation of a helical turn leads to a critical stress 25% higher than the planarization mechanism. Such calculations are in good agreement with previous experiments made with neutron irradiation [START_REF] Nicol | Matrix Damage in Iron[END_REF][START_REF] Hashimoto | Deformation mechanisms in ferritic/martensitic steels irradiated in HFIR[END_REF].

Conclusion

In this chapter we show the following:

• Irradiation has a major influence on the degradation of the mechanical properties of the RPV. Such influence is mainly characterized by an increase in the yield stress and a shift in the DBT temperature. • Mechanical properties of the RPV steel are affected by irradiation due to the presence of radiation-induced defects. There are main four types of defects; precipitates, SIA loops, cavities and solute clusters. • Defect's concentration and size are the main two factors that determine the effect of each defect to the embrittlement. In this study we concentrate on the effect of precipitates and SIA loops since they have the highest obstacle strength. In case of SIA loops, Burgers vector of the loop and character of the dislocation determine the nature of their interaction. • Interactions of defects with screw dislocations are reported to have a higher critical stress than that with edge dislocations. This is because of the nature of the observed mechanisms, also the thermally activated glide of screw dislocation.

Introduction

Dislocation dynamics provides the link between the atomistic and macroscopic scale. It is based on the linear elastic theory. In this theory, dislocations glide in an elastic continuum medium [1]. In case of a single dislocation, a mathematical solution for the stresses can be found analytically. However, in the presence of N dislocations, the problem becomes chaotic and no analytical solution can be found. This is because each dislocation segment generates a stress field that is affected by the stress field of all other segments in the simulation box. The main purpose of DD simulations is to provide crystal plasticity laws that depend on initial parameters (ex, crystallographic properties, temperature and the corresponding elastic constants, deformation regime, etc.).

The advantage of DD is that it can handle massive simulations of large population of dislocations and radiation-induced defects. The present PhD is mainly focused on this scale and make use of a DD code that will be described later in this chapter. This work is part of a multi-scale materials modeling approach done within the framework of the European project SOTERIA [2]. SOTERIA aims to study the plasticity in nuclear reactor pressure vessel and internals. It focuses on the effect of radiationinduced defects and microstructural heterogeneities using a combination of experimental and modeling work.

In this chapter, we present the history and categories of DD codes. In section (2.1), we present the general algorithm on the DD code, NUMODIS, used in this study. In section (2.2.2) we put emphasis on one particular aspect of this code, which is the set of the local rules and core reactions used to handle interactions of dislocations. My contributions to the developments of NUMODIS code are illustrated in this chapter.

DD codes: history and present

The history of DD goes back to the 60's of the bygone century. The very first 2D DD codes were developed by Foreman in 1967 and Bacon in 1967 after the introduction of self-stresses concept by Brown in 1964 [3][4][5]. The aim was to evaluate the critical shear stress of certain mechanisms. Among these mechanisms is the activation of a Frank-Read source, or Orowan's mechanism (a dislocation bypassing an obstacle). In the late 80's other 2D DD codes were developed to simulate the interaction between several infinite straight dislocations [6,7].

The beginning of three dimensional DD codes goes back to 1992 [8]. The aim was to develop a code capable of handling individual or collective 3D dislocation behavior in a more realistic manner. Different mechanisms were studied such as climb, cross-slip or junction formation. Later on, more codes were developed and each showed an outstanding capability to model various physical problems. In particular, mechanical behavior of ferrite laths, fatigue studies, strain hardening induced by forest dislocations or other crystalline defects like precipitates or grain boundaries. Detailed description of the history of 3D DD codes can be found in Kubin 2013 [9]. Currently, a limited number of 3D DD codes are available because of the efforts and time needed for the development of such codes. Dislocation lines are discretized into segments of constant line character [8]. The nodes of Lattice-based DD codes are found at the corners of a three dimensional lattice. Hence, nodes are authorized to move on lattice positions only. Consequently, the degrees of freedom (dof ) of the nodes are reduced to a finite number. Tridis (SIMaP, Grenoble INP) or microMegas (LEM, CNRS-Onera) are two examples of lattice-base DD codes. In the former, dislocations are discretized in a succession of Edge-Screw segments, where in the latter, dislocations are discretized in a succession of Edge-Mix-Screw segments.

On the other hand, nodal codes provide a better representation of the initial dislocation line than lattice-based codes. More specifically, segments orientation and nodal position vary in a continuous way. Therefore, complex geometries and reactions can be described with more precision in nodal codes. NUMODIS (CEA, France) is an example of such an approach. Paradis (LLNL, USA) is another example of nodal DD codes. For these reasons, the work of this PhD is done using NUMODIS code.

Most of the existing DD codes share the same concept of discretization, forces and velocity calculations. The major difference lies in the analysis of core reactions. For example, annihilation reactions, junction formation, cross-slip and other reactions. Detailed description of the general algorithm of the NUMODIS code and some other features are presented in the next two sections; (2.1) and (2.2.2).

NUMODIS: a nodal DD code

NUMODIS is a nodal dislocation dynamics code. It is jointly developed by the SRMA-CEA (Commissariat à l'Energie Atomique), the SIMaP (Grenoble INP) and the INRIA (Bordeaux I). This code is written using an Object-oriented C++ language and runs on Linux platforms. The general algorithm used in this code is illustrated in Figure 2. Stress calculation at the nodes of the dislocation configuration.

3. Force calculation on each dislocation segment with the help of the Peach-Koehler form. The complexity of the interaction is of order O(N 2 ). This is because each segment in the N-body system interacts with all the other segments.

4. Velocity calculation at the nodes based on the force calculation made in step (3). Velocity depends on the form of mobility laws. In NUMODIS, in the simplest case of a linear mobility law, a set of linear equations with the form KV = F are solved to find the nodes velocity, where K is the viscosity matrix, V and F are the nodal velocity and force vectors respectively. Each node is then displaced to its new position in the simulated volume.

5. Discrete events are then treated. Such events include interaction between dislocations or collisions of dislocations with different obstacles (other dislocations, precipitates, grain boundaries, etc.).

At each time step, the concept of power dissipation, following the concept of virtual work (see [10,11]), is used. When a dislocation-dislocation reaction takes place, the selected reaction pace is that with the maximum power dissipation. NUMODIS 1.0 is a sequential code that runs on one processor of a computer. There is another parallel version of the code called NUMODIS 2.0. It runs on several processors using hybrid OpenMP / MPI implementation. As a result, it can simulate larger spatial and temporal scales within substantially shorter simulation time. The parallel version (NUMODIS 2.0) is however still under development and does not include all the physical models today developed in NUMODIS 1.0. This is why, in this PhD, NUMODIS 1.0 was preferred. Some details of the general algorithm are presented in sections (2.1.1-2.1.7).

Discretization of dislocation line

Dislocation lines in NUMODIS are discretized into nodes connected by segments. Two types of nodes, physical and topological nodes are considered. Physical nodes are those connecting more than two segments, or at the intersection of two segments in different gliding planes. These nodes have a physical meaning. On the other hand, topological nodes are introduced to discretize the geometry and curvature of dislocation lines between physical nodes. Each topological node is therefore connected to two segments sharing the same Burgers vector and glide plane.

Line discretization is made as follows. Straight segments are connected through a set of nodes. Each segment is assigned a Burgers vector b, a unit vector tangent to its length ξ and a normal to the glide plane n. These vectors are essential to determine the force acting on the corresponding segment as well as the velocity of the nodes. To ensure a non-redundant discretization, a pair of nodes is connected only by one segment. Furthermore, the summation of Burgers vector at physical nodes must equal zero in agreement with the Frank's Rule [13]. In general, the Burgers vector of segments is kept unchanged during the simulation, but it can be modified in certain cases when the segment undergoes specific reactions. This is explained in detail in section 2.2.2.

A fine discretization of a dislocation line is achieved by decreasing the distance between nodes and therefore by increasing the number of topological nodes. As the CPU calculation time increases with the number of segments, it is of importance to optimize the line discretization in simulations. In particular, regions with high stress gradients need more nodes, as certain reactions may be hindered if the dislocation curvature is not well captured while using coarse segments.

Nodal force calculation

Nodal forces are calculated to find the response of a dislocation to an external applied stress or due to its own self-stress. Nodal forces are expressed as the negative gradient of the total energy of the dislocation with respect to the segment position (cf. Equation 2.1).

f i = - ∂E tot (r i , b ij ) ∂r i (2.1)
It is convenient to split the total energy into two separate terms. The first term represents the long range elastic field and the other is the core energy.

E total = E elastic + E core (2.2)
The elastic strain energy has the dominant contribution in comparison with the core energy. In case of isotropic elasticity, the self-energy per unit length of a dislocation of an infinite straight dislocation takes the form (2.3):

E elastic (θ) = µb 2 4π(1 -ν) (1 -ν cos 2 θ) ln( R 0 r 0 ) (2.3)
where µ is the isotropic shear modulus, ν is the Poisson ration, b is the Burgers vector, R 0 and r 0 are the outer and inner radii of the dislocation line elastic energy, respectively. One can notice that this equation is orientation dependent on the character of dislocations.

For reasons of simplicity, core energy in NUMODIS is assumed to be a small fraction (percentage α core ) of the elastic energy of an infinite straight dislocation [9,13].

E core (θ) = α core * E elastic (θ) = α core * µb 2 4π(1 -ν) (1 -ν cos 2 θ) ln( R 0 r 0 ) (2.4)
For more convenience, we introduce at this point a term called the core energy parameter ζ core . This constant quantity take the form:

ζ core = α core ln( R 0 r 0 ) (2.5)
Another approach to express forces on dislocation line is done by line tension approximation models. Generally, a curved dislocation line tends to straighten out to minimize its strain energy and length in a process analogous to an elastic string. This phenomenon is known as the line tension. Line tension is defined as the change in the elastic self energy of a dislocation with respect to its length. The concept of line tension is applied in two ways in NUMODIS, either an isotropic or anisotropic model. The latter is derived using the following equation [14]:

F core (θ) = E core (θ) * T + dE core (θ) dθ * N (2.6)
where T and N are two vector tangential and normal to the dislocation line, respectively. Once elastic and core energies are determined, nodal forces can be expressed in terms of the negative gradient of the two dislocation line energy contributions as:

f i = f elastic i + f core i (2.7)
The elastic force can be calculated with the well-known formula proposed by Peach-Koehler in 1950 [16]. It should be noted that Equation (2.1) is a formulation to calculate the nodal force which is not directly related to the Peach-Koehler force. The latter takes the form (2.8):

f P K elastic (x) = (σ(x) • b) × ξ(x) (2.8)
where σ(x) is the stress at position x. In this equation, the stress is calculated as the sum of external applied stress and internal stress. In the case of curved dislocations or dislocation loops, internal stresses are calculated as follows [1].

σ αβ (x) = µ 8π ∂ i ∂ p ∂ p R[b m ε imα dx β + b m ε imβ dx α ] + µ 4π(1 -ν) b m ε imk (∂ i ∂ α ∂ β R -δ αβ ∂ i ∂ p ∂ p R)dx k (2.9)
where µ is the isotropic shear modulus, ν is the Poisson ration and ijk is the Levi-Civita notation; a permutation term used for the asymmetric properties.

Mathematical singularity appears in the dislocation elastic energy when R goes to zero, i.e. at the dislocation line. In order to avoid such problem, a method based on the non-singular model proposed by Cai et al. [17] is used. This standard method is now used in many simulation codes like the ParaDis code [18]. Cai et al. distributed the Burgers vector of the dislocation through a regularization function. R is therefore replaced by R a = R 2 + a 2 0 . In this formulation, even if R equals zero, R a is not null. The singularity is therefore eliminated near the dislocation core. The strength of this method is that it requires only one parameter (a 0 ), which is commonly defined as a core regularization parameter. This makes it convenient for implementation in DD simulations and this method is used in this PhD work.

There are several simplified models in the literature of Equation (2.9). All the above mentioned models reduce the triple integral into a double integral. This last integral can be simplified to a line integral over a closed path using Green's function as follows:

σint = σ ij (x) = C ijkl dislocation ε lnh C pqmn G kp,q (x -x )b m dx h (2.10) G ij (x -x ) = 1 8πµ [δ ij ∂ p ∂ p R - 1 2(1 -ν) ∂ i ∂ j R] (2.11)

Nodal mobility laws

In this section, calculation of the nodal responses to the applied forces is presented. This process is governed by predefined nodal mobility laws, whose purpose is to make a link between the force acting on a dislocation segment and its velocity. This relation is an important input of dislocation dynamics simulations and depends on many physical parameters especially when thermally activated mechanisms are involved [19].

The most general form of mobility function used in DD simulations is the Newtonian motion of dislocations. In such case, the system reaches a steady state condition in a relatively short time. This solution does not apply under extreme strain rates (e.g. shock propagation) where the driving force is not constant [20]. This point is considered as a limitation of DD codes to simulate deformations under high strain rates. A general form of mobility function is Equation (2.12).

v i = M(f j ) (2.12)
where v i is the velocity, M(f j ) is a mobility function. As mentioned before, mobility laws depend on many factors that are not explicitly defined in equation (2.12). The simplest form of mobility law is v = τ b/B, where v is the nodal velocity, τ is the shear stress, b is the Burgers vector and B is a drag coefficient. Assuming that, we can define a viscosity matrix on the nodes, the nodal velocity can then be calculated with Equation (2.13).

KV = F (2.13)
where K is the viscosity matrix, V is the corresponding nodal velocity vector and F is the force vector exerted on nodes. In NUMODIS, the viscosity matrix K has large dimensions. The computational time needed to solve this system is high. To speed up this calculation, NUMODIS divides the matrix K into blocks of sub-matrices dedicated to topological or physical nodes. Such solution takes advantage of the fact that the topological sub-matrix is a sparse band matrix.

In this work, we considered mobility laws used for BCC materials, such as (Fe-α). This is because it represents the main class of materials used in nuclear reactor pressure vessel. A simple mobility law is defined for edge and screw dislocations (see Bulatov et Cai [20]). An edge dislocation has two possible moves in the crystal, glide and climb. The mobility law for a pure edge dislocation is then expressed as:

B(ξ) = B eg (m ⊗ m) + B ec (n ⊗ n),
when ξ⊥b, (2.14) where B eg is the drag coefficient in the glide direction and B ec is the drag coefficient in the climb direction. Because climb is controlled by diffusion processes [13], the climb mobility used in the simulation is commonly phenomenological and does not account for all the details of the real physical process. The unit vectors are defined as n ≡ b × ξ/ b × ξ and m ≡ n × ξ. n is a unit vector perpendicular to the glide plane, m is the unit vector in the glide direction and ξ is the tangential unit vector of the segment as indicated in

Figure (2.3).

There is an equation to specify the drag coefficient tensor for all segments orientations in [20]. This drag tensor is singular and cannot be inverted. The problem is illustrated by Arsenlis et al. [18]. The reason is the alignment of dislocation segments. If two segments connected by a discretization (topological) node become collinear, the nodal drag tensor becomes ill-conditioned. No drag coefficient is indeed explicitly specified for the configuration forces on a discretization node along the segment direction. This problem is mitigated when a specific contribution is added to B. One can take the following form:

B(l) = B g (l)(m ⊗ m) + B c (l)(n ⊗ n) + B I t (ξ ⊗ ξ) (2.15)
where B t is a tangential drag viscosity associated with the discretization nodes. B t should be much less than B eg , B ec and B s . This ensures that the mechanical behavior of the system is not affected.

The tangential viscosity also prevents the rotation of the dislocation around itself. This implies that this term was introduced only for numerical stabilization and does not correspond to any physical process.

Time integration

The nodal equation of motion is based on the forces and mobility laws calculated at the nodes (the function g i ({r j })) in the form:

v i = dr i dt (2.16)
There are several numerical techniques to solve ordinary differential equations. Each method has its advantages and disadvantages in terms of accuracy, stability and computational time. One of the simplest methods is the Euler forward expressed in Equation (2.17).

r i (t + ∆t) = r i (t) + v i ∆t (2.17)
where r i (t) is the position of the node i and at time t. ∆t is the elementary time step. In case of stiff equations, a very small time step is needed to ensure stability and accuracy of results. Therefore, time steps must be carefully chosen in explicit Euler forward method. One of the advantages of Euler forward method that it is computationally inexpensive. In addition, it is easy to implement in DD codes. Therefore, this method is used in most DD code (including NUMODIS) as the numerical time integrator.

The Trapezoid method is a combination between Euler-forward and Euler-backward. It is a good example of implicit methods. Other methods could be used in DD codes (not used in NUMODIS for the moment) such as Newton-Raphson (matrix inversion) or Euler-Trapezoid (explicit method but it is computationally expensive) [18,[20][21][22].

Topological changes and operations

In DD simulation, dislocation segments interact with one another or with other radiation-induced defects. During this reactions, the length and curvature of dislocation lines evolve. Consequently to ensure a better representation of the dislocation line, the number of discretization (topological) nodes need to be modified on-the-fly. 1The key point in topological changes is the conservation of the Burgers vector. It is also important to ensure that each pair of nodes is connected by only one segment. Topological operations in NUMODIS, as in other nodal DD codes, can be summarized into two main algorithms, merge and split algorithm. The former is used to delete an existing node. The two neighboring segments are then merged. In this process, the number of degrees of freedom is reduced. The latter is responsible for adding a node in the middle of a segment. Thus, the number of degrees of freedom is increased. All other complex topological changes (i.e., in case of junction formation) are discussed in section (2.2.2).

To ensure an adequate representation of the dislocation lines in DD simulations, discretization nodes are added or deleted on-the-fly. Introduction and elimination of nodes are subjected to a simple and robust criterion. Such criterion depends on the length of the segment. Whenever the segment is longer than (L max ), a node is added along the line in the middle of the segment. The initial segment is therefore cut into two segments. On the other hand, if it is shorter than (L min ) a node is deleted. If L min and L max are not carefully chosen, repetitive discretization may happen during a simulation. This is why, it is recommended to set the values of the lower and upper boundaries of the segment's length as L max >2L min . The two discretization operations are illustrated in Figure (2.5).

Periodic boundary conditions

Computational methods at the microscopic or mesoscopic scales are inherently limited in terms of simulation cell size due to elevated time cost. The concept of periodic boundary condition (PBC) is a commonly used strategy to reduce the simulation volume. In this approach, the simulated volume is assumed periodic and we consider only a reference volume surrounded by a set of identical replicas.

In the present work, the elastic interaction of each dislocation segment with all the other segments of the simulation cell is performed under the concept of minimum image convention [20]. This interaction is actually computed with respect to the nearest replica of each segment. This concept is illustrated in Figure (2.6). In the NUMODIS code, one can actually go beyond the minimum image convention and choose the minimum number of replicas considered in all three axial directions. This choice nevertheless comes at higher computational cost, with no guaranteed convergence unless specific conditions are met [20]. In the following, we use the minimum image convention unless specified elsewhere.

It is important to avoid a cubic simulation box when PBC is used in massive simulations [24]. This is done to eliminate the auto-annihilation of dislocation lines with a dislocation image of opposite sign from the other side. This process is illustrated in Figure (2.7). Solution to optimize the shape of the simulated volume in relation to this problem can be found in [25].

Control modes in NUMODIS

There are several ways to apply a mechanical deformation on the simulation box. Constant strain rate or constant stress solutions are commonly used to apply external load. Control modes are implemented in NUMODIS. The desired control mode is a user-defined option.

Constant applied stress on the simulation box is the simplest loading condition in Dislocation Dynamics simulations. Although useful for some problems, it is not well suited to measure critical stress values for many phenomena. In contrast, constant strain rate control is preferred [12]. Since most performed computations in the present PhD work involve dislocation-obstacle interactions analogous to existing MD studies, the constant strain rate mode was used.

The plastic deformation is mainly related to the movement of dislocations in the materials [1,13]. It is then calculated as the summation of all segments displacements along the dislocation lines. εplastic =

segments b i n j + b j n i 2Ω * dA (2.18)
where b, n are the Burgers vector and a vector normal to the glide plane, respectively. dA is the area swept during the movement of a dislocation segment and Ω represents the simulation volume. The strain rate at each time step is defined as the derivative of equation (2.18) with respect to time. In addition, the averaged stress is evaluated at each time step. Any mismatch between the total strain and the plastic strain is by definition εelastic = εtotalεplastic . This difference is applied at the subsequent time step following Hook's law (cf. Bulatov et Cai [20]).

Recent developments of NUMODIS

During this PhD, original developments have been made in the NUMODIS code. The aim of such development is to cope with different physical phenomena and to reproduce MD simulation conditions. In this section we present such developments.

Nodal mobility laws

In section (2.1.3), we showed earlier that Arsenlis et al. introduced a term of tangential viscosity for topological nodes to prevent spurious motion along the segments to which there are connected. Alternatively, no specific viscosity is attributed to the physical nodes. Unlike topological nodes, physical nodes are not virtual and cannot be replaced. We should remind here that physical nodes are nodes at the intersection of two or more segments of different Burgers vector or glide planes. In some reactions we observed in DD simulations, physical nodes tend to propagate rapidly. Such behavior is thought to have two origins; i) a real specific viscosity exists on physical nodes, which is not implemented in NUMODIS (this possibility is discussed in chapter 4), ii) it is the result of an ill-conditioned matrix when two dislocation segments connected by a physical node are quasi collinear. Hence, to prevent the rapid propagation of physical nodes, we proposed to assign a specific viscosity term to the physical nodes. This is analogous to the tangential viscosity added in Equation (2.15) for the topological nodes.

B(l) = B g (l)(m ⊗ m) + B c (l)(n ⊗ n) + B I t δ I (l)(ξ ⊗ ξ) + B I p δ I (l) (2.19)
Four different models of physical viscosity are now proposed in NUMODIS. These models are based on the nature of the physical node; either connected to only two segments or more, or if the viscosity term is added in an isotropic or anisotropic matter. Different methods in the code can be used to apply the viscosity to the physical nodes.

• Model 1: Isotropic viscosity on physical nodes of two or more connections.

• Model 2: Isotropic viscosity on physical nodes of two connections.

• Model 3: Anisotropic viscosity on physical nodes of two connections.

• Model 4: Anisotropic viscosity on physical nodes of two or more connections.

Although these four models are implemented in NUMODIS, they are not used by default. The desired model to apply should be specified by the user in the input files. The choice of certain model is made based on the simulated mechanism. For example, Model 2 is applied only to reproduce a reaction between screw dislocation and a radiation-induced loop. More details on this additional mobility laws are given and illustrated in chapter (4).

Core reactions in NUMODIS

This section is part of the general algorithm of NUMODIS, but since I contributed to a major development in this part of the code during my PhD, it is presented here in a separate section.

General description of the "split node" algorithm

In a typical DD simulation, a dislocation undergoes certain reactions either with obstacles or with other dislocations. Hence, a DD code should be able to handle the occurrence of an interaction before it takes place e.g., collisions, displacement or junction formation. The split algorithm intervenes at the step of nodal velocity calculation. This algorithm determines the evolution of a physical node connected to multiple segments, which includes the possibility to split this node in multiple new nodes connected by new dislocation segments. Following the pioneer algorithm used in the Paradis code, this algorithm first constructs a comprehensive list of possible outcomes, then selects the most likely one based on the concept of maximum power dissipation (see [10,11]).

At each time step of the simulation, all the physical nodes are examined to determine if they are splittable or not. In NUMODIS, a splittable node must have at least three connected segments(referred to as connections in the code parlance) and must not be arbitrarily pinned (i.e., F-R source). A typical splittable node of multiple-connections is shown in Figure (2.8). The general algorithm of split node is described as follows:

• Determination of possible outcomes (thereafter called split options) for each splittable node.

• Calculation of force and velocity for each possible split option.

• Calculation of the power dissipation for each split option as well as consistency verification. The split option with the highest dissipation is chosen.

This algorithm has geometrical, crystallographic and energetic aspects. They are discussed in detail in the following sections.

Geometrical and crystallographic aspects

This part of NUMODIS corresponds to the first step of the split node algorithm. The geometrical aspects are determined as follows:

1. The connections to each splittable node are identified and sorted by glide plane.

2. If a glide plane contains more than three connections, they are sorted by the angle in that plane with respect to the first connection of the list.

3. All the possible connections dispatched into several (non-empty) groups are considered. This process explores all dispatches from one unique group for all the connections, to a number equal to the number of connections (one connection per group). Such a dispatch is referred to as an assignation. In practice, three groups at the most are explored by NUMODIS for assignation and assessment of the properties of the corresponding nodes and junctions.

The purpose of the third step is to dispatch each group of different connections of splittable nodes into assignations using a special routine. This routine is based on permutation operators (combinatorial part). Its objective is first to ensure that an assignation is not treated more than once. The second objective is to prevent any geometrical conflicts that could occur when there are four or more connections in a given plane. This possible issue is illustrated in Figure (2.9).

Crystallographic aspects

For each assignation, the NUMODIS code then builds the corresponding topological configuration in terms of new nodes and new junctions, which are referred to as virtual nodes and virtual junctions all along this algorithm. This configuration is itself call a split option. As indicated before, the number of groups (and therefore virtual nodes) goes from one to three. These split options are shown in Figures (2. Splitting a node into two virtual nodes is the canonical approach in nodal DD simulations [20] and is implemented as such in the Paradis [18] and NUMODIS codes. Interestingly, a careful examination of the interaction mechanisms observed in the MD simulations, on which this study is based on, leads to the need for a more general algorithm, which could lead to three (or even more) virtual nodes. One example of such splitting is illustration in Figure (2.12) and will be discussed in more detail in chapter (4).

Once constructed, the split options are tested from geometrical and crystallographic aspects. The following algorithm is therefore valid for split into one, two and three virtual nodes.

1. Compute the degrees of freedom of each virtual node considering the glide planes of its connections and computing their intersection.

2. Compute the properties of each virtual junction (Burgers vector, glide plane) and update the degrees of freedom of each virtual node accordingly.

3. Check the consistency of the split option: a virtual junction can obviously only be formed if at least of the virtual nodes at its ends is mobile. If not, the split option is discarded.

4. Compute the generalize coordinates of virtual nodes. They contain the dof and position.

5. Store the split option for later analysis from a dynamical point of view.

In point (1), the dof of virtual nodes is calculated without considering any virtual junction in between. The possible dof are 0,1,2.

• a zero dof is when the node is sessile and cannot move in the space, for example a pinned node.

• a dof equals to one means that the node can move only in one direction, for example if the node lies on the intersection of two different planes, it can move in the common direction between the two planes. An example of such a situation is when one segment has a constraining plane and it is sessile as in the case of sessile lock.

• a dof equals 2 when the node is restricted to move in all the plane in the space.

In point (2), the properties of a possible virtual junction are determined. Such properties are the direction, glide plane and Burgers vector. The latter is calculated by the summation of Burgers vector of all other real connections on a virtual node. A split option of two virtual nodes is shown in Figure ( In case of split into two virtual nodes, the summation of Burgers vector b 1 +b 2 = b 3 +b 4 is logical. However, in case of split into three virtual nodes or more, the summation of Burgers vector of the created junctions could be different. Moreover, the created junction between two virtual nodes can have many possible glide planes. Such planes are filtered based on a criterion that only those planes that contain the dof of nodes A and B are chosen.

In point (3), the possibility to create a virtual junction between two virtual nodes is verified. In the case of split into two virtual nodes, if the dof of both virtual nodes equal zero, then they cannot move apart. No junction is therefore formed. In the case of split into three virtual nodes, the criterion is modified. Two out of three nodes must have a non-zero dof in order to form two virtual junctions. In points (4,5), the generalize coordinates of virtual nodes are computed. They contain the dof and position of the nodes. The split option is then stored and the virtual nodes are added to the virtual graph to be later analyzed from a dynamical point of view. So far, the direction of the virtual junctions created in a split option of two or three virtual nodes is known. In some case, the direction of such junctions is unknown. In order to find a specific direction to build a split option, the non-linear solver intervenes. This point is discussed more in the following section.

In summary, at this point in the code, only geometrical and crystallographic treatment of the different split options is performed. In the next section, the physical or energetic verification is performed to choose the most favorable configuration. The notion of virtual junctions and the technique to calculate forces and dissipation is explained in the next section.

Dynamical aspects of the split node algorithm

When performing DD simulations, the system should always evolve in a way that lowers the energy. Following the original work by Cai et Bulatov, the outcome of our splitting algorithm is based on the maximum dissipation criterion. This method is therefore general and does not rely on a set of hardwired rules. We first recall how to compute the power dissipation of any given dislocation configuration based on the forces and velocities at nodes, which includes here virtual nodes using Equation (2.20).

P = i F i • v i (2.20)
where P is the dissipated power due to the formation of this configuration, F is the force acting on a segment i attached to a node n and v is the velocity of that segment. The summation is done on all existing split options.

In order to calculate the power dissipation of a split option, nodal forces and velocities should be determined. In the following, we present the importance of the core energy contribution of virtual junction to the force calculation. We also present how NUMODIS handles the situation where the direction of the virtual junction is unknown. Here for simplicity reason, we first assume the case of two virtual nodes connected by one virtual junction. Then the algorithm goes as follows:

1. Calculate the total force acting on each virtual node by adding the forces due to all its connections, as well as the core energy due to its virtual junctions.

a) Determine if the virtual junction is a predefined junction direction (i.e. at the intersection of two planes).

b) Use the non-linear solver if the direction or the sense of a virtual junction is unknown (i.e. it lies in one plane).

2. Calculate the velocities at all the nodes, including virtual nodes and compute the total power dissipation.

3. Verify that no collisions will take place in the first (∆t) between the created virtual nodes, or one node with the connections of the other node.

Here it must be noted that the core energy contribution to the force acting on virtual nodes is a key feature to avoid possible deadlocks in DD simulations.

To illustrate this point, we consider a node with four connections at time = t (cf. Figure (2.10) and assume that our algorithm decides to split this node into two virtual nodes separated by a virtual junction considering no core energy for the virtual junction. At time = t + ∆t, the force acting on the new nodes now includes a core energy that tends to reduce the junction length, possibly leading to the collapse of these two nodes. This would lead to an endless split and merge deadlock that could prevent other possible split options to be selected. This is explained in Figure (2.14). The direction of the virtual junction should be known at the moment when the core energy contribution is included, which is the more likely case. For example, assume that connections (1,4) and (2,3) in Figure (2.14a) are in glide planes (P1) and (P2), respectively. If the split option contains two virtual nodes (cf. Figure 2.14b), the virtual nodes then have one common degree of freedom at the intersection of the two glide planes. This case is or instance observed in the Lomer-Cottrell configurations as illustrated in Figure (2.15). [13]. Connections (1,3) and (2,4) are found in planes 1 and 2 respectively. The virtual nodes (in red color) have a dof equal one, they can move only in one direction. The virtual junction (in yellow color) is therefore a 1D junction.

In some reactions, the direction of the virtual junction is however unknown. For example, the interaction of two coplanar dislocations may lead to the formation of a junction whose direction is not known beforehand. It is therefore impossible to account for the core energy of the virtual junction in a simple manner, as its contribution depends directly on the velocity of its two virtual nodes which is not known a priori. In such a case, the core energy and thus the nodal forces become function of the velocity. Therefore, we have to solve a non-linear system of the form presented in Equation (2.21).

KF (V ) = V

(2.21)

The non-linear solver we use is based on the Newton-Raphson method. In order to enhance the convergence, a back-tracking algorithm is used. More details on this algorithm can be found in [27]. One of the pitfalls of the non-linear solver is its inability to converge for all proposed slip-nodes possibilities. Regardless of numerical difficulties (possible bugs, bad initial guess, etc.), this solver may actually not converge if there is no possible velocity field satisfying Equation (2.21). The convergence problem in the non-linear solver can be explained by three possible categories: programming, numerical or physical.

• programming: despite about two hundred benchmarks, implementation bugs and unforeseen situations are still present. Several bugs were identified and fixed during my PhD, but some problems in the solver probably remain.

• numerical: ill-conceived (initial guess, non-robust iterative process, etc.). since the non-linear solver is based on Newton-Raphson method, the divergence can be related to the initial guess used by the solver which is far from the solution. In such case, a possible divergence occurs despite the number of performed iterations.

• physical: a split-option can simply not happen because it is energetically not physical, at least considering the current energetic parameterization. This situation typically occurs if the core energy of the possible junctions is too high, regardless of their orientations as explained in Figure (2.14). A priori identification of such a situation is however hard to perform, and leads to spurious (and sometimes harsh) oscillations, which must be mitigated.

In the case where the solver does not converge for one of the previously mentioned reasons, the direction of the virtual junction is hard to determine. The non-linear solver is allowed to perform 50 iterations to converge toward a solution. When the solver did not converge in 50 iterations, we exit the corresponding subroutine and this split option is excluded from the list of split options. So we skip this particular split option and search for another one that converges.

In point (2), we mentioned that at this stage, the nodal velocities are computed within a finite element approach, which involves the resolution of a possibly huge system in the form KF = V . In order to expedite this calculation and additionally to the subdivision of K, we further divide the physical nodes into splittable nodes and non-splittable nodes as described in Figure 2.16. This simplifies tremendously the calculation of their velocities for all the split options using the Schur complement method [28]. We also note that the non-linear system of Equation (2.21) can be limited to this splittable part of the system and makes it computationally tractable. In point (3), we verify that the dissipated power is positive in agreement with thermodynamics laws. No energy should be created during split option reactions. In conclusion, the power dissipation should always be positive in this approach. This criterion is tested when the core energy of the virtual junction is added to the split process. In summary, the sense of the virtual junction should not be reversed when the core energy contribution is added.

In the scope of preventing deadlocks or any anomalies in the system, certain reactions should be anticipated before they occur. In point (4), a test is used to verify that no collisions will take place between the nodes formed, or one node with the segments of the other node in the first (∆t). An example of such a situation is shown in Figure (2.17).

According to nodal velocity, the position of a node can be calculated at time = t+∆t. If node B is going to interact with one of the connections of the other node, the split option that leads to this situation must be rejected. Figure 2.17 -Schematic representation of a configuration where a created virtual node has a velocity vector that leads to the interaction with one of the connections of the other virtual node in the next time step. The dashed black arrow is the velocity vector of node B.

If the split option is accepted, the nodal velocities are stored. Finally, the split option that leads to the maximum power dissipation is chosen. The graph is then updated and all virtual nodes become real nodes in the system.

Special treatment of a specific configuration

Lastly, a special case is encountered when a splittable node has three connections or more (cf. Figure 2.18a). A possible outcome may involve the splitting of this configuration into two virtual nodes, separated by a virtual junction sharing the same Burgers vector as the connection attached to a node, but gliding in another glide plane. This evolution bears similarities with a cross-slip mechanism that would not be a thermally activated mechanism in this case. As we will see later in chapter (4), this splitting mechanism is important and has been observed in MD simulations. We, nevertheless, restrict this option when the angle θ between the connection and its Burgers vector is above a critical angle. If not, it would first lead to numerical instabilities because the connection and the virtual junction are almost aligned. Here we use a local rule introduced as cosminXslip to prevent any cross-slip if the angle between the two segments is less than 15 • . It is a hard-coded value and has no particular scientific justification. This prevents any numerical instabilities in velocity calculations. Such instability was explained in section (2.1.3). The viscosity matrix becomes ill-conditioned and could result in an overestimation of the nodal velocities. (a) the initial configuration of splittable node of three connections, (b) split option into two virtual nodes, where segments 1 and 2 are connected to a virtual node and segment 3 is connected to the other virtual node. θ is the critical angle between the connection and its Burgers vector. A minimum value is authorized in order to prevent numerical instability.

The virtual junction and connection (3) in Figure (2.18b) should have the same Burgers vector simply by the summation of in and out vectors. However, the virtual junction and connection (3) must have different glide planes. This criterion ensures that the virtual junction is not aligned with the side arm of the virtual node. A cross-slip of the virtual junction is therefore prohibited. The artificial cross-slip reaction presented here is different from the more usual thermally activated cross-slip processes.

Conclusion

Dislocation dynamics is a powerful tool to study the interactions of dislocations with radiation-induced defects. It provides the link between atomistic and macroscopic scale. Our methodology and the originality of this work is the attempt to reproduce the individual screw dislocation-defects interactions previously reported in molecular dynamics studies, but within the framework of the elastic theory of dislocations.

Our DD simulations are made with the NUMODIS code that I contributed to improve during my PhD. It is a 3D nodal based code, where each dislocation is discretized into a set of topological and physical nodes. The uniqueness of this code is its ability to manage and control collisions and core reactions between dislocation segments and many different irradiation defects. This is done through a set of generic algorithms with the minimum number of local rules.

During this PhD, several models and algorithms were implemented. The three major and key developments made in NUMODIS are:

• Implementation of a new viscosity model. This model increases the stability of nodal velocity calculation. An application of this technique is illustrated in chapter (4).

• Separation into three virtual nodes and not just only two. This configuration is required in some of the reactions we investigated.

• The core energy contribution of virtual node is considered in the force summation of split options. This allows avoiding deadlock configurations between successive time steps.

One should pay attention in future simulations when the artificial viscosity on physical node is used. Since the nodal energy balance is affected and certain reactions are altered.

Several geometrical and crystallographic aspects are considered in the process of proposing different split options. Later in this manuscript these split options are verified from a dynamical point of view. Regarding the core energy contribution of the virtual nodes, the modifications made in the code allow estimating more correctly the direction of the virtual junction. In the case where the direction of the virtual node is unknown, a non-linear solver is activated to identify possible reaction path. This non-linear solver is based on Newton-Raphson method. One of the pitfalls of the non-linear solver is its inability to converge for all proposed slip-nodes possibilities. This non-linear solver needs some additional work in the future.

Introduction

One major objective of this work is to make DD as accurate as reasonably possible, accurate meaning here "as close as possible to MD". This approach was coined by J. Marian as "atomistically-informed dislocation dynamics" [1]. Before tackling the issue of radiationinduced defects, it is first important to ensure that basic features of DD simulations such as mobility laws and line tensions are well implemented and calibrated with respect to MD. The calibration method is discussed in this chapter, as well as the developments to which I contributed in the NUMODIS code to reach this objective. More specifically, in section 3.1, we present the implementation of periodic boundary conditions. Then, in section (3.2) we present the implementation and a validation of thermally activated glide of screw dislocations. Finally, in section 3.3, we reveal the importance of carefully choosing the parameters related to the dislocation energy by performing a parametric study on Orowan mechanism.

Periodic boundary conditions

Periodic boundary conditions are used in most of the MD simulations to which we want to compare our DD simulations. PBC in DD are described by the existence of image replicas of the original simulation cell in the periodic direction. It also implies the ability of nodes along the dislocation line to cross the boundary and reappear from the other side. Such boundary conditions were implemented in NUMODIS code during this PhD to enable better comparison with MD simulations. Additionally, we added the possibility to go beyond the minimum-image convention for the stress field calculation and choose the number of replicas in all three directions in the input files. For example, a node on the dislocation line can interact with one or two layers of replicas of the segments in the reference volume (see Figure 3.1). Using multidimensional replicas substantially increases the simulation time, as well as the precision of the calculated stresses. In the following chapter, the minimum image convention is used by default for all simulation results.

Thermally activated glide

In this section, we present the implementation of a thermally activated model for the glide of screw dislocations. The model is based on the double-kink mechanism. Lastly, the implementation is tested and verified.

Double-kink model in DD simulations

It is known that at low temperatures, the flow stress in α-Fe, like other BCC materials, exhibits strong temperature dependence. This dependence is due to an increasing difference of mobility between screw and edge dislocations when decreasing the temperature. Indeed, screw dislocations must face a strong lattice resistance when edge dislocations mobility is essentially temperature independent. It is due to differences in the core structure between edge and screw dislocations [2].

Screw dislocations move between local energetic minima called Peierls valleys. This motion is a process done through the nucleation and migration of kink pairs along dislocations. Kink mobility is generally known to be very high in metals. Hence, at low temperature, plasticity is thought to be controlled by the nucleation of the kink pairs along the screw dislocation lines [3,4]. The formation of a double-kink on a screw dislocation using MD simulation is shown in Figure There are several models in literature that describe the double-kink formation, i.e. in [6][7][8]. One of the double-kink models is the one presented in Equation (3.1) as in [7]. The velocity of screw parts is calculated based on this thermally dependent kink pair nucleation model. The specificity of this model comes from the fact that it considers the stress variation along screw dislocation segments. The line integral of the double-kink mechanism is therefore described by the following equation:

υ screw = hν D b l * 2 L 0 exp - ∆G(σ(x)) k B T dx (3.1)
where h is the average distance between Peierls valley, ν D is Debye frequency in the order of 10 13 s -1 , b is Burgers vector and l * is the critical length of a double-kink. The integration is done on the range [0 -L], where L is the critical length for nucleation of a kink pair. ∆G is the activation energy, σ is the average stress calculated on the screw segment, K T is Boltzmann's constant and T is the temperature.

The exponential part in Equation (3.1) represents the nucleation probability of a double-kink on the dislocation line. The activation energy is calculated by Equation (3.2) [9].

∆G(τ ) = ∆G 0 1 - τ τ 0 p q (3.2)
where ∆G 0 is the energy barrier at T=0 K. p and q are dimensionless activation parameters. This equation was later modified by Garcia 2011 [10]. They added an entropic correction factor A and the equation becomes:

∆G(τ, T ) = ∆G 0 1 - τ τ 0 p q -Aκ B T (3.3)
For pure iron, the constant A is chosen to result in an athermal transition temperature around 300 K for ∆G 0 1 eV . The activation parameters p and q used in the later simulation are equal to 0.5 and 2, respectively, as in [11].

Implementation in NUMODIS

The general algorithm for the implementation of the thermally activated glide of screw dislocations on NUMODIS is illustrated in Figure (3.3). The different steps are then discussed in details. Based on Equation (3.1), the length of the screw dislocation is an important parameter that must be evaluated to calculate the velocity of screw dislocation segments. This is because the probability to nucleate a kink pair is directly related to the effective length of screw segments. In case of curved dislocations, it is essential to determine the character of the dislocation section along the line. Furthermore, the length of the segment with a screw character must be defined.

To determine the nature of segments we apply a criterion based on the direction of tangent line of the segment. If the dislocation tangent is parallel or nearly parallel to the Burgers vector it is considered as a screw segment. This criterion is mathematically expressed as | ξ(ι) • b |= 1 + δ. In this equation ξ(ι) is the direction of the segment, b is Burgers vector and δ is a user defined value. In the present work, we consider the dislocation segment to be of screw character if the critical angle θ between Burgers vector and the tangent is less than a specific value. Otherwise, it is considered as a non-screw segment. The value of the critical angle θ is taken to be 10 • in NUMODIS. This choice of (θ) affects the length of the detected screw segment. In the case of a straight dislocation of length L 0 , if θ = 10 • then the detected length of a screw segment is equal to 1.015L 0 . According to Equation (3.1), the difference in the predicted velocity is 0.1%. In NUMODIS each line has its own tag with the nodes contained in that line. The algorithm starts scanning the dislocation line from the first node until the last node in the line. Figure (3.4a) represents a curved dislocation line where three screw segments are identified in such configuration, nodes in the range [1][2][3][4][5], and [START_REF] Terentyev | Simulation of the interaction between an edge dislocation and a < 100> interstitial dislocation loop in α-iron[END_REF][START_REF] Shi | Interaction of <100> dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF][START_REF] Terentyev | Reactions between a 1/2<111> screw dislocation and <100> interstitial dislocation loops in α-iron modelled at atomic scale[END_REF][START_REF] Liu | Molecular dynamics simulations of the interactions between screw dislocations and self-interstitial clusters in body-centered cubic Fe[END_REF][START_REF] Hirsch | Point defect cluster hardening[END_REF][START_REF] Hashimoto | Deformation mechanisms in ferritic/martensitic steels irradiated in HFIR[END_REF]. The mentioned segments satisfy the tangent criterion. They are all aligned with the Burgers vector and therefore considered as screw segments.

Figure (3.4b) represents a delicate situation in the case of a dislocation loop. One must pay attention that node [1] is a common node between the two lines [1][2][3][4][5][6][7][8][9][10] and , though they should be considered as a single line for the following reason. As stated previously in Equation (3.1), the velocity of the screw segment depends on its length. Therefore, the velocity of two adjacent screw segments of length (a) and (b) respectively, if considered as separate lines, is lower than the velocity of one screw segment of length (a+b). As a result, the two segments [1][2][3][4][5][6][7][8][9][10] and are considered as one line in the range .

In the following we discuss a special configuration that is shown in Figure (3.5). This is a definition of a zero length segment as described in Lebon 2011. A zero length segment consists of one node between two edge segments under the following condition. The two segments should be at the same side with respect to the plane perpendicular to the glide direction, as in Figure (3.5a). In order to detect zero length screw segments the nodes are consequently analyzed. The analysis is done on the ensemble of nodes n -1, n and n + 1 on each dislocation line. Figure 3.5 -Schematic representation of the configuration of a zero length screw dislocation segment. (a) the two dislocation segments are on the same side with respect to the glide direction, the red node is then a zero length segment, (b) the two dislocation segments are not on the same side with respect to the glide direction, the red node is not a zero length segment. b is the Burgers vector and d is the glide direction. The dashed line is the plane perpendicular to glide direction (after Lebon 2011) [12].

When the length of each screw segment is determined, the stress profile is then calculated. The line integral in Equation (3.1) is then calculated using the trapezoid method. The last step then is to calculate the velocity of screw dislocation parts based on the double-kink model. In NUMODIS, nodal velocities are calculated by solving a set of linear equations in a matrix form. It is done by the inversion of a viscosity matrix of the form KV = F . First, the velocity of screw segments is calculated using Equation (3.1) as proposed by Guyot et al. in [7]. Second, the calculated velocity is imposed (υ imp ) on the nodes found on screw parts similar to Dirichlet boundary conditions in solver of the KV = F matrix. The velocity is imposed by multiplying the terms k ii in matrix K by a huge number (∼ 10 10 ) and replace F i in the vector F by k ii υ imp . The system is then solved to find the values of the V vector. Such mathematical technique is used in [12] and is called matrix penalization. An advantage of this implementation is that it requires the minimum changes in the general algorithm of NUMODIS.

    k 11 . . . x 1n . . . k ii . . . k m1 . . . k mn        υ 1 υ imp υ n    =    F 1 k ii υ imp F n    (3.4)
Once the velocity of screw parts is calculated, the code advances one time step. The process is then repeated at each time step of the simulation.

Validation of the model

In order to validate the implementation of this model, a simple configuration of Franck-Read source of screw nature was tested at T = 250 K. When loading is applied, screw parts in the Frank-Read source propagates in the crystal at a lower rate that the non-screw parts. The reaction is shown in Figure (3.6). It should be noted that in the F-R multiplication process, the difference in velocity of screw and non-screw parts is observed before the critical curvature at normal and high strain rates. At very low strain rates, dislocations are at a quasi-static condition. In such condition, screw and non-screw parts are then in equilibrium and have the same velocity. The lattice friction on screw dislocation is higher at low temperatures. It tends to stabilize at higher temperature. The average difference between the velocity calculated using NUMODIS and that predicted by the double-kink model is 0.35%. As stated in Equation (3.1), screw velocity increases with temperature. This is because the probability of double-kink nucleation increases with temperature.

In conclusion, thermally activated glide mechanism was successfully implemented in NUMODIS, but further investigation and analysis are required to use this model in future DD simulations.

Line tension calibration

This work is submitted as an article to the journal of "Computational Materials Science" Shukeir et al. 2019 (in press).

Simulation of Orowan

Introduction

Dislocation Dynamics (DD) simulations inherently rely on a continuum description of dislocations based on the linear elasticity (e.g. [13]), whose accuracy in capturing the strain and stress fields around dislocations has been proven experimentally over the years (e.g. [6]). As such, this simulation method nevertheless fails at describing the discrete atomic positions in the immediate vicinity of a dislocation and consequently its energy. As pointed out by Bulatov and Cai, this remark holds even for non-singular elasticity models [14][15][16] where the stress field remains finite in the core of the dislocation. Most DD simulations therefore partition the total energy of a dislocation into an elastic and a core energy contribution [17], the latter being added to account for whatever is left by the former in the dislocation core. There are noticeable efforts in feeding atomistic information to these models [1,[18][19][20]. For example, Dang et al. fitted a 0 parameter on MD simulations for dislocation loops in Al. At 300 K, a 0 is calibrated to be equal to 0.6 of the lattice parameter [20]. Nevertheless, limited efforts have been made so far to identify the relevant parameters on atomistic data, nor to assess their influence on classical dislocation mechanisms such as the Orowan process. The objective of this study is therefore to address this issue and provide guidelines to choose the corresponding parameters.

Among the various approaches used in current DD codes ( [17, 21-25]) to account for dislocations total strain energy, the choice was made here to use the isotropic nonsingular elastic model of Cai et al. [16] in combination with a basic core energy model. In the absence of a comprehensive atomistic data base, the latter is simply based on the orientation dependent line tension approach [26]. The former offers many mathematical and computational advantages and is currently widely used. As it will be later discussed, these choices limit the relevant energy parameters to two values. We therefore believe that this study should be of use and applicable to almost all the DD codes.

The Orowan mechanism [27] was selected in this study as our reference case. This mechanism, which is the bypassing of strong precipitates by a dislocation line moving in a given glide plane, was preferred to other elementary phenomena for the following reasons. First, it has an important impact on plastic strain hardening [28,29] and has been widely studied in the literature at the mesoscopic scale using dislocation dynamics [30][31][32][33][34]. Secondly, it does not rely on nonphysical arbitrarily pinned dislocation ends by contrast with the Frank-Read source [35] and can therefore be simultaneously investigated using molecular dynamics (MD) [29,36,37]. Finally, the Orowan mechanism questions simultaneously the underlying dislocation strain energy through its effective line tension and the stress required to bow out the side arms of the dislocation, as well as its stress field and the corresponding elastic interaction between the arms as the dislocation bypasses the precipitates [30]. Previous studies investigated the effect of precipitate size and distribution on the mechanical properties, but to our knowledge, no study investigated in a systematic manner the influence of the dislocation strain energy parameters used in DD simulations. This paper is organized as follows. In the next section we present a brief overview of dislocation elastic theory and parameters appearing in the dislocation strain energy definition used in most DD simulations, as well as the standard model used to predict Orowan stress in the case of impenetrable obstacles. In section (3.3.3), we describe the details of the parametric study performed in this article. The corresponding results are presented in section (3.3.4). The last section is dedicated to a discussion and concluding remarks.

Theoretical background

In linear theory of elasticity, it is convenient to split the total energy E total of a dislocation into two separate terms, one for the long range elastic field E elastic and the other for the core energy E core [6].

E total = E elastic + E core (3.5)
Most of the total energy of a dislocation comes from the elastic strain energy contribution, while the core energy is reported to constitute a few percent of the total elastic strain energy [13,38]. Assuming isotropic elasticity, the self-energy per unit length of a dislocation stored in a cylindrical ring of inner radius r 0 (core radius) and outer radius R 0 in the case of infinite straight dislocations is

E inf elastic = µb 2 4π(1 -ν) (1 -ν cos 2 θ) ln( R 0 r 0 ) (3.6)
where µ is the isotropic shear modulus, b is the Burgers vector and ν is the Poisson's ratio. The dislocation character term θ is the angle between the Burgers vector and the tangent vector along the dislocation line. From Equation (3.6), we simply see that the elastic energy of a dislocation depends on i) the dislocation core definition via r 0 , ii) the initial and boundary condition via the definition of R 0 and iii) the dislocation character (edge or screw) via the pre-logarithmic energy constant [6,13].

Because of the inability of the elastic theory to represent the core energy of a dislocation and in the absence of comprehensive atomistic information, the choice was made here to define the core energy as being directly proportional to the elastic energy per unit length of an infinite straight dislocation through a single parameter, α core , such as:

E core = α core E inf elastic (3.7)
This approach is consistent with several available atomistic studies (e.g. [13,39]) in which the core energy amounted to a few percent of the elastic energy. For convenience, we introduce at this point a term called the core energy parameter ζ core . The latter constant quantity take the form:

ζ core = α core ln( R 0 r 0 ) (3.8)
The core energy can therefore be written as:

E core = ζ core µb 2 4π(1 -ν) (1 -ν cos 2 θ) (3.9)
The driving force controlling dislocation dynamics is mainly a function of the loading stress and internal stresses generated by each dislocation line. The calculation of the internal stress field associated with an ensemble of curved dislocations requires the integration of the following equation [6,16,40]:

σ αβ = µ 8π ∂ i ∂ p ∂ p R[b m ε imα dx β + b m ε imβ dx α ] + µ 4π(1 -ν) b m ε imk (∂ i ∂ α ∂ β R -δ αβ ∂ i ∂ p ∂ p R)dx k (3.10)
where R is the distance between points α and β on the dislocation line and ijk is Levi-Civita notation.

A major difficulty in calculating the stress field of dislocations is related to the fact that Equation (3.10) is divergent when R approaches zero. Different solutions have been proposed to eliminate such singularity as in [15,16,41]. Among those solutions, the one proposed by Cai et al. [16] features some specifications which makes it particularly convenient for DD simulations. In brief, it modifies the singular stress field solution through a mathematical transformation where R in Equation (3.10) is replaced by R a = R 2 + a 2 0 . This transformation implies the definition of a new regularization parameter a 0 coined as the dislocation core width parameter.

Hence, the dislocation strain energy definition used in this study reduces to two major parameters: the dislocation core width parameter a 0 and the core energy parameter ζ core . In order to assess the impact of these parameters, we perform DD simulations on the Orowan mechanism.

In a seminal work, Bacon et al. [30] investigated the bypass of a periodic row of impenetrable spherical obstacles by an infinite dislocation with the Orowan mechanism [27]. From an analysis of the influence of the obstacle size (D) and the inter-obstacle distance (L) on the flow stress, they defined the following equation (hereinafter referred to the BKS model):

τ = A µb L [ln D + B] (3.11)
In Equation (3.11), A is a pre-logarithmic factor that depends on the character of the dislocation. A = 1 and A = 1/(1 -ν) for edge and screw dislocations, respectively. B is a fitting parameter evaluated to 0.7 in the BKS paper. This model equation is widely used in the literature as a reference for Orowan-like interactions [33,34].

Simulation technique

Simulations were carried out using NUMODIS [25,42], a 3D nodal dislocation dynamics (DD) code based on the isotropic elastic theory of dislocations. Dislocation lines are represented by a set of nodes, interconnected by straight segments nodes, whose properties are their Burgers vector and their glide plane.

This study is performed on a single crystal of BCC Iron at 300 K. At this temperature, the lattice parameter equals 0.2855 nm, b the Burgers vector of slip systems 1 / 2 [111](1 10) equal 0.2475 nm and the shear modulus µ equals 63 GPa. The x, y and z axes of the simulated volume are oriented in the [111], [ 11 2] and [1 10] crystallographic directions, respectively. Periodic boundary conditions are applied in the [111] and [ 11 2] directions, while no particular conditions are applied to the boundary surfaces in the [1 10] direction since they are normal to the dislocation glide direction.

Following a robust methodology previously used to simulate the Orowan process [33,34,43], one 1 / 2 [111](1 10) edge or screw dislocation is introduced in the periodic volume in front of an impenetrable spherical obstacle with diameter D and cut by the dislocation glide plane at its center (Figure 3.8). This obstacle is considered as an incoherent inclusion, whose elastic properties are equal to the surrounding crystal. No specific stress field is therefore associated with the obstacle. The spacing L between periodic images of the obstacle is varied by changing the size of the simulation box along the direction parallel to the dislocation line. The dimension of the simulated volume in the glide direction is systematically adjusted to allow for Orowan loop formation before the dislocation line reach the boundary of the periodic volume. The z-axis is systematically adjusted to be three times the obstacle diameter. Pure shear stress is applied in the b direction to impose a constant strain rate which corresponds to a constant dislocation velocity in the glide plane close to 3 m/s. For reasons of simplicity and to allow for comparison with previous computations made in FCC materials, we consider a simple linear over-damped mobility law similar to [44], in the form:

v s = τ b B (3.12)
where τ is the effective resolved shear stress and B is a viscosity coefficient set to 8 × 10 -5 Pa.s. The latter quantity accounts for dissipating processes like dislocation-phonon interactions. All simulations are done in a quasi-static condition and tests have been made to verify that the viscosity coefficient value has no influence on the computed critical stress. Here, it must be noted that although thermally activated bypassing mechanisms are reported in the literature in case of small obstacles, for reasons of simplicity no thermally activated dislocation property like, dislocation cross-slip or climb, are considered in the present simulations. All the material parameters used in the DD simulations are consistent with prior molecular dynamics simulations found in [45][46][47].

In order to investigate the influence of the two core parameters, a parametric study of the Orowan mechanism was conducted. The range of the many combinations we tested are summarized in Table (3.1).

Simulation results

The results of all the calculations of the Orowan critical stress we conducted are plotted in Figure (3.9). For comparison with the BKS model, the Orowan stress is plotted as a function of the harmonic mean of the inter-obstacle spacing and obstacle diameter, (L -1 + D -1 ) -1 . A large dispersion of the values is found when changing the parameters controlling the dislocation elastic energy. In the edge dislocation case and for different values of the harmonic mean, an increase of 100% in the critical stress is observed. Such dispersion is found to be even more important in the case of screw dislocations, where the critical stress increases of approximately 280% for the different sets of parameters. The critical stress of an edge dislocation is observed to be quasi-independent of the Poisson ratio ν, while a non-linear dependence of the form 1/(1 -ν) is observed in the case of a screw dislocation. This observation is in agreement with the line tension model proposed by deWit [26] that accounts for the effect of the dislocation character and the simulation results first reported in [30]. Figure (3.10) as well shows that when ν equals zero, the Orowan stress is independent of the dislocation character and the Orowan stress is independent of the dislocation line character as expected.

The second analyzed simulation parameter is the core energy parameter ζ core . As shown in Equation (3.7), ζ core has a linear influence on the total strain energy of the dislocation and therefore on dislocation line tension. The same tendency is observed on screw and edge dislocations (see Figure 3.11). More precisely, increasing ζ core by a factor of 10 increases the critical stress by about 33%. This reveals a fairly high dependence of DD simulations on the definition of the dislocation core energy. This point has been probably been overlooked in several existing studies. ) can be motivated by the dislocation energy expression computed by Cai et al. [16] in the context of the non-singular theory. Indeed, the core width parameter a 0 appears in the denominator of the logarithmic term, as if r 0 was replaced by a 0 in Equation (3.6). Consequently, the simulated Orowan stress is decreased by approximately 18% and 16%, for screw and edge dislocations respectively, when increasing the core width parameter a 0 by a factor of 3. This result raises doubts about the relevance of DD simulations, which use regulation parameters 10 to 100 times larger than the core radius of the dislocations defined by atomic simulations for computational reasons.

Lastly, we show that in agreement with the BKS model for a given set of simulation parameters, the Orowan stress is logarithmically dependent on the harmonic mean of inter-obstacle distance and obstacle size. This result is presented in Figure (3.13) with a set of simulation parameters (a 0 = 0.45 nm, ζ core = 0.86 and ν = 0.33) and is in good agreement with the solution found in previous simulations. As discussed in [30], the evolution of the Orowan stress is here well described with the help of a harmonic mean between the inter-obstacle distance L and the obstacle diameter D. When L is much larger than D the average tends to L and the required stress to overcome the precipitate decreases, since the line tension of a dislocation in inversely proportional to its length and vice-versa.

Discussion and concluding remarks

In this work, we studied in α-iron the interaction of infinite 1 / 2 [111](1 10) screw and edge dislocations with a periodic array of impenetrable obstacles of different size and spacing. These configurations are standard simulation problems previously used to study the Orowan mechanism with MD and DD simulations. The reported results are in good agreement with previous studies, but reveal a large dispersion of results when changing the simulation parameters used to define the dislocation strain energy. More specifically, for the different solutions of simulation parameters, the calculated Orowan stress can vary by 100% and 280% for edge and screw dislocations, respectively.

In order to reveal the collective effect of the studied parameters to the Orowan stress, we propose to combine them in a new equation. This equation represents the sum of two contributions. It reflects the strain energy decomposition into elastic and core energy contributions (see Equation 3.5). The first contribution is similar to the one defined in the BKS model, but includes the contribution of the core width parameter a 0 , as it modifies the elastic strain energy in the framework of the non-singular dislocation theory. The second term accounts more specifically for the contribution of the dislocation core energy via the parameter ζ core . Indeed, as revealed in Figure (3.11), this contribution contributes significantly to the dislocation line tension and cannot be neglected when modeling phenomena involving dislocations curvature like the Orowan process. The following modified Orowan equation is the outcome of a fit made with more than 1000 simulations run with different parameters.

τ Orowan = µb L A 2π [1.23( ln D a 0 -0.18 )+ζ core ] (3.13)
The harmonic mean term of L and D appears in the numerator in conformity with the BKS model, while the core width parameter (a 0 ) appears in the denominator since τ Orowan has an inverse dependence on this term. The additional right hand side term of the equation account for the linear contribution of the core energy to the Orowan process. As illustrated in Figure (3.14), prediction made with Equation (3.13) gives excellent results with a correlation factor of 0.9965 when considering all our simulation data. It should be noted that Equation (3.13) is generic and can be used to interpolate any type of results on the Orowan mechanism. Thus, it becomes possible to perform a reverse analysis to identify the parameters of a DD simulation that can reproduce experimental data or other simulation results. In the following, an example of such adjustment is given. Equation (3.13) is used to define the parameters that are needed in our DD simulation to reproduce the results on the Orowan stress obtained by Lehtinen et al. with MD simulations [48]. These results are discussed and compared to the BKS model. The results of DD simulations when parametrized with the help of Equation (3.13) show an excellent match with the Orowan stresses computed with MD simulations. More precisely, the difference between the stress found in DD and MD simulations at different obstacles spacing is on average less than 2.4%. It is worth noting also that this difference is 53% between the BKS model and the MD simulation results. Such a difference could be interpreted as coming from the existence of thermally activated phenomena in the MD simulations that cannot be reproduced with any approach based only on the elastic theory. Our study, rather suggests that the dislocation strain energy that controls the dislocation dynamics in the MD simulations is simply different from the definition used initially in the BKS model.

In conclusion, it should be noted that the parameterization of the dislocation strain energy in the framework of the non-singular dislocation theory is governed by the idea that the total dislocation energy should remain unchanged when changing the regularization amplitude. Hence, when increasing the dislocation core width parameter in a DD simulation (to avoid using tiny time steps), the dislocation core parameter is usually also increased to keep the total energy constant. As illustrated by our study, such solution without a systematic investigation of the parameterization effect could be problematic. Indeed, we show that the dislocation dynamics are not affected in a simple manner by both parameters. For instance, some new DD simulations to be published in a forthcoming paper have shown that when modeling the interaction between a dislocation and radiation-induced loop defects, a drastic effect is observed on the dislocation dynamics by changing a 0 a few percent. Such parameterization change may even lead to the modeling of very different contact reactions.

In summary, a large dispersion is observed in DD simulation results depending on the choice of parameters used to control the dislocation strain energy. Such uncertainty has to be eliminated so that DD simulation results can be compared with MD and DD simulations or even experiment. A parametric study is proposed to incorporate the effect of essential simulation parameters in an equation form useful to predict the Orowan stress. This model reveals the existing coherence between the many simulations we performed and provides the means to calibrate rigorously future DD simulations. The model was tested and validated by a direct comparison with MD simulations. This work opens the door to more quantitative comparison between DD simulations and other simulation techniques.

Conclusion

In conclusion, some methodological aspects were successfully implemented, tested and used in NUMODIS during this PhD work. Such progresses allow for much better comparison between DD and MD simulations.

First, periodic boundary conditions are implemented. PBC are essential to minimize the calculation time as the size of the simulation cell is reduced. A specific parameter is assigned to control the number of replicas used to calculate forces based on the concept of minimum image convention. The implementation of PBC in DD simulations allows to compare our results more rigorously and precisely with MD simulations. PBC has several advantages; mainly to reduce the size of the simulation box and consequently the computation time.

Secondly, the thermally activated glide mechanism was successfully implemented in NUMODIS. We presented a robust algorithm to detect screw segments at each time step. The velocity of screw parts is then calculated based on a double-kink nucleation model in [7]. In order to validate the implementation we tested a simple case of a Frank-Read source expansion in pure iron. In addition, the effect of temperature was investigated on the glide of an infinite straight screw dislocation. The temperature dependence found using NUMODIS is in agreement with the analytical solution of the double-kink model. Although the thermally activated glide was successfully implemented, the interactions in the next chapter are based on a linear mobility law (cf. Equation 3.12). The idea was to insure that the core reactions (presented in section 2.2.2) are correct and the algorithm is working properly for different types of interactions. The thermally activated glide can be therefore used in future simulations using the NUMODIS code.

Lastly, in section (3.3), we showed that a large dispersion in DD simulation results can be observed. This dispersion is mainly related to the influence of simulation parameters associated to the dislocation strain energy definition. These parameters are the core energy parameter ζ core and the regularization parameter a 0 from the non-singular theory of elasticity. A parametric study is proposed to understand the detail influence of each parameter. This study is based on Orowan mechanism. We propose an equation to predict Orowan stress based on the simulation parameters selection. This model provides the mean to rigorously calibrate future DD simulations. The model was tested and validated by a direct comparison with MD simulations. This work opens the door to more quantitative comparison between DD simulations and other simulation techniques. The values of the main parameters ζ core and a 0 are 0.56 and 0.25 nm respectively. These values are generally used in the next chapter, unless specified elsewhere.

It should be noted that the parameterization of the dislocation strain energy in the framework of the non-singular dislocation theory is governed by the idea that the total dislocation energy should remain unchanged when changing the regularization amplitude. Hence, when increasing the dislocation core width parameter in a DD simulation (to avoid using tiny time steps), the dislocation core parameter is usually also increased to keep the total energy constant. As illustrated by the present study results, such solution without a systematic investigation of the parameterization effect is potentially the source of problems. Indeed, we show that the dislocation dynamics are not affected in a simple manner by both parameters. Such parameterization change may even lead to the modeling of very different contact reactions.

In conclusion, DD simulations include, like every mesoscopic simulations, some constitutive and partly phenomenological rules that must be educated to reproduce the results of simulations made at atomic scale. In Tables (3.2, 3.3), a list of the adjusted parameters in DD simulations to reproduce the MD results for iron at 300 K. Some of the these parameters are fitted, others are directly extracted from atomistic simulations. Burgers vector Plane B (10 -5 Pa.s) τ 0 (MPa)

1/2 1 1 1 {1 1 0} 8 10 1/2 1 1 1 {1 1 2} 8 10 1/2 1 1 1 {1 2 3} 8 10
1 0 0 {1 1 0} 80 300 1 0 0 {1 0 0} 9000 300

Introduction

In the previous two chapters, we presented several models and developments made in the NUMODIS code. In this chapter, such models and developments are used to investigate, with DD simulations, the individual interactions between screw dislocations and radiation-induced loops. Here it must be noted that, some reactions could not be reproduced correctly. We therefore present an analysis of the problems.

In this work, we focus on the interactions of the 1 0 0 square SIA loops (cf. section 4.1) and the 1/2[1 1 1] hexagonal SIA loops (cf. section 4.2). Different reaction mechanisms are observed when modeling SIA loops of different sizes and orientations. The details of the reaction process are discussed. In the following, the DD simulation results are divided in different sections corresponding to the different orientations of the SIA loop.

Several important points are discussed in this chapter. First, we highlight the importance of the split node algorithm on the reproduced reaction mechanisms in the DD simulations. In particular, we show the effect of split into three virtual nodes to investigate the interaction between a screw dislocation with a [0 0 1] square loop at 300 K. This split algorithm, presented in section (2.2.2), plays a key role in the case of interactions with large 1/2[1 1 1] SIA loops.

Secondly, a detailed analysis is presented for the formation of a helical turn in the case of interaction with small [1 1 1] SIA loops. The possible effect of an additional viscosity term on the physical nodes on the helical turn formation mechanism is here discussed. This section shows the need for additional input from molecular dynamics simulations.

In the last section, we present the results of a new MD investigation. The results of such simulations show the importance of including the thermally activated cross-slip mechanism and the twinning/anti-twinning mechanism in DD simulations.

Interaction of screw dislocation with 1 0 0 loops

In this section, the interaction of a screw dislocation with 1 0 0 loops is studied using DD simulations. The reference study for this type of interaction was discussed in section (1.5.4) (cf. [1]). We only present the calculations made at 300 K which corresponds to the room temperature. This is because of the need to understand the increase in the yield stress at room temperature.

The simulation conditions used in the DD and MD simulations are similar: The size of the simulation box is (24 * 28 * 24) nm 3 . Axes x, y and z are oriented along the crystallographic directions [1 1 2], [1 1 0] and [1 1 1], respectively. Periodic boundary conditions are applied in the [1 1 1] and [1 1 2] directions, while no particular conditions are applied to the boundary surfaces in the [1 1 0] direction since they are parallel to the simulated dislocation glide plane. A screw dislocation is introduced parallel to the x-axis and a square loop of 128 SIAs is introduced at the center of the simulation box. A constant shear strain rate ε = -10 7 s -1 is applied in the x direction. Such strain rate corresponds to a steady state velocity of 27 m/s. The simulated configuration is shown in Figure (4.14). simulations [1]. The green surface represents the glide plane of the mobile screw dislocation, while arrows at the top and bottom surfaces indicate the direction of the applied shear stress. The SIA loop is introduced at the center of the simulation box and intersects the glide plane of the mobile screw dislocation. The arrows on the loop indicate the Burgers vector.

Additional DD simulation parameters are taken as follows. The discretization method ensures that the minimum length of a segment is 2.5 Å and the maximum is 6.0 Å. A small time step equal to 0.5 ps is taken to prevent numerical instabilities due to the explicit time integration algorithm used in our simulations. The stress calculation is done within the framework of non-singular elasticity theory (cf. section 2.1.2). The regularization parameter a 0 is equal to 2.5 Å (cf. Equation 3.10), while the core parameter is set to ζ core = 6.85 in agreement with the line tension identification (cf. Equation 3.8).

In this study, two main reactions were observed according to the orientation of the burgers vector of type 1 0 0 . The first reaction (A3) results in the restoration of the original loop (see section 4.1.1), while the formation of a helical turn is observed for the reaction B3 and C3 (see section 4.1.2).

Reactions ending with the restoration of the original loop

This case is observed when a screw dislocation is interacting with a loop of b = [0 0 1] with sides lying in {1 0 0} planes (A3). The loop is subjected to strong drag coefficient of 0.8 MPa.ns and a high threshold stress of 300 MPa. The presented values are based on trial and error calculations made by Shi et al. when the interaction of an edge dislocation with [0 0 1] loops was studied using NUMODIS [2]. Hence, one must keep in mind that the mobility of such loops is reduced. The initial configuration is shown in Figure At the moment of collision, the split node algorithm is involved. Node (A) (visualized in Figure 4.2b) is a node to be split according to the criteria presented in section (2.2.2). The list of split options to be considered after the collision is shown in Table (4.1). For each split option we show the number of virtual nodes (VN), the possible dof of each virtual node, the virtual junctions (VJ), the Burgers vector and planes of the created junctions and the total energy dissipation of the split option.

#0 1 1 0 - - 0.0021 #1 2 1,1 0 - - 0.0023 #2 2 1,1 1 [1 1 1] (0 1 0),(1 1 0) solver did not converge #3 2 1,1 1 [1 1 1] (1 1 0) solver did not converge
The non linear solver did not converge in case of split options (#2, #3). We already know from MD simulations that such reactions are not expected for such configuration as they will eventually lead to a helical turn formation due to the creation of a 1 1 1 {1 1 0} junction.

In this table, split option #0 corresponds to the unchanged configuration (see 2.10). Split option #1 implies the creation of two virtual nodes but no junction in between. It basically correspond to the dislocation crossing and shearing the loop without any particular junction formation. Despite the fact that this solution is not the one observed in MD simulations, it is the solution selected by our DD simulation. Indeed, the virtual power dissipation we calculate in option #1 (0.0023 eV /ps) is found to be higher than that for option #0 (0.0021 eV /ps). The final configuration we obtain after split option #1 is shown in Figure Figure (4.3b) is useful to see that the screw dislocation mobility is affected by the stress field of the loop. The critical shear stress to unpin the dislocation is 208 MPa (∼20% lower than the CSS found in MD). This type of reaction is illustrated in Terentyev et al. [1].

The only way to reproduce the reaction observed in MD simulations (see Figure 1.37) in the DD simulation is to allow node (A) in Figure (4.2b) to dissociate into three virtual nodes connected by two virtual junctions. This development in the code was previously illustrated in section (2.2.2). At this stage, it should be pointed out that before the implementation of the split into three nodes, I "manually" forced the code to take such option by creating the configuration as in Figure (4.4a). As a result of the creation of three nodes, the dislocation propagates along the side of the loop (see Figure 4.4) in good agreement with MD results. This leads to the conclusion that the interaction was physically possible, yet our algorithm could not treat the split into three nodes and the problem was not related to a numerical problem in the non-linear solver.

The split into three nodes is successfully implemented in our DD code. An updated list of split options including a split into three nodes is illustrated in Table (4.2). With this new degree of freedom in the reaction path, there is a total of 16 possibilities to consider in the split node algorithm. Here, we highlight only a few examples to illustrate the analysis made in the split node algorithm. Further details about the outcome of the proposed split options are presented in Appendix (A).

#0 1 1 0 - - 0.0042 #1 2 1,1 0 - - 0.0043 #2 3 0,1,1 2 [1 1 1],[1 1 1] (0 1 1),(2 1 1) 0.0043 #3 3 0,1,1 2 [1 1 1],[1 1 1] (1 0 1),(1 0 1) 0.1233 #5 3 0,1,1 2 [1 1 1],[1 1 1] (1 2 1),(3 2 1) solver did not converge #6 2 0,1 1 [1 1 1] (0 1 1) solver did not converge #15 3 0,1,1 2 [1 1 1],[1 1 1] (2 3 1),( 2 1 1

) solver did not converge

As can be seen from table (4.2), some split options imply splitting into three virtual nodes connected by two virtual segments. Such solutions imply different combinations of slip systems and connections between the nodes. Split option (#3) is finally selected because this reaction solution leads to a larger energy minimization. As a result of this solution involving 3 nodes, one of the created segments is aligned with the segment 2 in the loop (see The first half of the reaction observed in MD simulation is then correctly reproduced in DD with the help of split into three nodes. The second part of the reaction implies the cross-slip of the [1 1 1](1 0 1) segment to convert the [0 0 1](0 1 0) on the right side of the loop to a segment of type [1 1 1](1 0 1). Once this step is done, the dislocation continues the glide in the (1 1 0) plane and overcomes the loop with a process similar to the Orowan mechanism. This second part of the reaction observed in MD simulations could not be reproduced in the DD simulations in the absence of cross-slip model. To end this process, the node in the lower part of the loop is to be splitted into two nodes. A [1 1 1](1 0 1) junction has to be created in between the node.

In order to provide an insight of the expected reaction, the cross-slip reaction missing in the DD simulations was manually reproduced by forcing the code to chose a specific configuration. In the following, the [0 0 1](0 1 0) segment is manually converted to a [1 1 1](1 0 1) segment. Then, the load is applied to this configuration without any relaxation of the handmade configuration. The outcome of this simulation is shown in Figure (4.6). Thanks to this handmade reaction, the original loop is restored with a bi-loop in the [1 1 1](1 0 1). This reaction is similar to MD predictions. The critical stress to overcome the new configuration is 530 MPa, which is ∼%100 higher than the CSS reported in MD. Though the original loop is restored at the end of this reaction, the gap between the critical stress found in MD and in DD is huge. No explanation could be found to justify such difference. The list of proposed split options when the first collision occurs between the dislocation and the loop is listed in Table (4.3). At the end of this reaction, the original loop is transformed into a [1 1 1] loop, with a small surviving triangle of the original character [0 1 0]. The helical turn could not be reproduced in the DD simulations. This result comes from the absence of the creation of [1 1 1] segment by the split node algorithm as shown in Table (4.3). Neither the final configuration in Figure (4.8) nor the intermediate state are similar to MD simulations.

Reactions ending in a helical turn formation

#0 1 1 0 - - 0.1057 #1 2 1,1 0 - - 0.1078 #2 2 1,1 1 [1 1 1] (1 1
The critical shear stress to overcome the helical turn in MD was reported to be equal to 430 MPa, while in DD it is 476 MPa. Despite the difference in (B3) reaction between MD and DD, the critical stress to unpin the dislocation is only ∼ 10% higher in DD than MD calculation. This result suggests that this weak reaction is dominated by long range elastic interaction and the details of the contact reaction may not be important. For this reason the solution we reproduce in DD simulation was considered acceptable for massive simulations and no further investigations was carried out for such reaction.

The second reaction, found with MD simulations, resulting in the formation of a helical turn is the interaction with a [1 0 0] loop (case C3). The initial configuration is shown in Figure (4.9a). When the stress is applied, the screw dislocation approaches the loop but no interaction occurs as illustrated in Figure The CSS to overcome the loop in MD is found to be equals 430 MPa. The CSS in DD is ∼ 73% higher than in MD simulations. This large difference in CSS is related to the creation of the [1 1 1](1 3 2) segment. In Table (4.4), we examine the list of possible tested split reactions that lead to choosing such option by the split node algorithm.

#0 1 0 0 - - 0.4049 #1 1 1 1 [1 1 1] (1 1 0) 0.4117 #2 1 1 1 [1 0 0] (0 0 1) 0.4055 #3 1 1 1 [1 1 1] (1 3 2) 0.8946
The creation of such segment in the cross-slip glide plane should not be favored. Ultimately, to prevent the previous cross-slip, the node to split (#1 in Figure 4.11) is manually changed to a pinned node. Hence, the node could not be split. As a result the original dislocation continues the glide in the (1 1 0) plane. The result of the reaction is now shown in Figure (4.12).

Conclusion

In summary, the interaction mechanism between a [1 1 1] screw dislocation with 1 0 0 square loops at 300 K has been studied. Three different orientations were investigated in DD; [0 0 1], [0 1 0] and [1 0 0] named A3, B3 and C3 respectively. The case A3 in MD is reported to be a restoration of the original loop. The first attempt in DD shows no interaction between the dislocation and the segments of the loop. This is because no proper split option was tested in the DD simulation for such reaction. However, the critical stress equals 208 MPa which is ∼ 21.7% lower than CSS in MD simulations. This reaction could not be reproduced in DD without the split into three nodes. The critical shear stress after using split into three nodes equals 530 MPa. In conclusion, the use of split in 3 nodes is not always of a great advantage from a quantitative point of view. This is because, hard configurations could be created during the process that lead to high critical stress to unpin the dislocation from the obstacle.

Regarding the case B3 and C3, the formation of a helical turn is observed in MD. This helical turn is not observed in both cases in the DD simulations. The critical stress in the case B3 is 476 MPa, which is very close to the 430 MPa found in MD. Though the reaction mechanism is not the same, this result is considered good enough for massive simulations. Indeed, when running massive simulations, the critical stress for individual interactions is more important than the details of the reaction path.

In case C3, the creation of a segment in the (1 3 2) plane favored the cross-slip of such segment. This segment expanded along the original dislocation leading to its elimination. The critical stress to unpin the new segment from the loop is high (750 MPa) compared to MD (430 MPa). When such cross-slip is manually prevented, the new critical stress is 500 MPa in much better agreement with MD simulations. A summary of the critical stress we calculated with the DD simulations in the three different cases tested in this section is shown in Figure (4.13). In conclusion, the reaction mechanism is affected by the size of the obstacle, orientation, angle of attack of the dislocation and other factors. Our goal to reproduce the same reaction with 1 0 0 loops in MD and DD simulations could not be achieved. However, for the sake of calibration of DD simulations based on the individual interactions, obtaining a critical stress in agreement with MD is sufficient. The reaction mechanism becomes more important when performing massive simulations, where each obstacle can be crossed several times by different dislocation lines. For, this reason, some phenomena such as the clear band formation can hardly be investigated with the present simulation code.

Interaction of screw dislocation with 1/2[1 11] loop

In this section, the interaction of a screw dislocation with 1/2[1 1 1] loop is studied by DD. The reference MD simulations for these types of interactions were discussed in section (1.5.4) (cf. [3]). We only present the results of the interactions at a temperature of 300 K. Two reactions are mainly observed; the formation of a helical turn and a loop planarization mechanism. The DD simulation parameters used in this section are the same as those mentioned in section (4.1). The topological discretization ensures that the minimum length of a segment is 2.5 Å and the maximum is 6.0 Å. A small time step equal to 0.5 ps is used to prevent oscillation in the nodal velocity calculation. The stress calculation is done within the framework of non-singular elasticity theory (cf. section 2.1.2). The core width parameter a 0 is equal to 2.5 Å (cf. Equation 3.10), while the core regularization parameter is set to ζ core = 6.85 (cf. Equation 3.8).

The SIA loop in this case is glissile and is found to be highly mobile. It can migrate along its crowdion axis 1 1 1 . At the beginning of the simulation, the loop migrates below the glide plane then in the direction of the dislocation. Such behavior results from the elastic attraction between the loop and the dislocation. The loop at some point comes in contact with the dislocation and a first junction is formed. This junction belongs to the [0 1 0](1 0 1) slip system. This feature is illustrated in Figure (4.15). Though the first formed junction is common for all loop sizes, the complete reaction was found to depend on the size of the SIA loop in agreement with MD simulations. A planarization mechanism is observed for the larger loops of 127 and 271 SIAs, while a helical turn is formed for the smaller loops of 37 SIAs. The details of both reaction mechanisms are discussed in the following sections.

Reactions ending with the planarization of the original loop

In section (1.5.4), we showed the MD results of the interaction between a 1/2[1 1 1] screw dislocation and large [1 1 1] loops of 127 and 271 SIAs. A planarization mechanism is observed, where the loop is transformed into a sessile [0 1 0](1 1 0) loop. In this section, we present the results of such interaction using DD simulations.

As shown in Figure (4.15), the formation of the first junction for the interaction of 1/2[1 1 1] loop with 1/2[1 1 1] screw dislocation is identical in the case of loops of 37,127 and 271 SIAs. This junction is sessile and it is a [0 1 0](1 0 1) segment type. In the case of a loop of 127 SIAs, the complete reaction mechanism observed in dislocation dynamics is shown in Figure (4.16). In the first snapshot, the formation of a helical turn is observed. This is similar to the results discussed previously in the case of a small loop of 37 SIAs. The critical shear stress to unpin the dislocation of the loop is 1006 MPa. This value is 287% higher than the CSS reported in MD for the same reaction.

The helical turn formation in the case of the interactions with large loops in DD simulations is contradictory to MD simulations. The [0 1 0] sessile segment should not be transformed into a [1 1 1] segment. Indeed, the latter facilitates the successive transformation of the segments of the loop to a helical turn of b =[1 1 1] segments. The correct mechanism found in MD is reported in Figure (4.15b). The segment (1) of the loop should cross-slip under increasing stress and react with the dislocation arm on the right side. As a result of such interaction, the interacting segment is transformed into a [0 1 0](1 1 0) segment. The other segments of the loop undergo the same transformation until the plane of the loop is completely changed.

In order to investigate the origin of the problem found in DD simulations, we precisely analyzed the solutions tested in the split node algorithm described in section (2.2.2). In this reaction, the split node into three virtual nodes is made possible. Nevertheless, no split solution with three nodes is finally investigated in the simulation. Indeed, as explained in section (2.2.2), the investigation of a particular split node solution requires the determination of the directions (tangent vector) of all virtual junctions created for this configuration. If one of the created virtual nodes has a dof=2, the node can move in two directions. In such case, the virtual junction connected to this node can have any direction in the plane. Therefore, the direction of the virtual junction cannot be determined from a crystallographic point of view, since the virtual junction is not bounded by specific crystallographic directions. This direction is then determined with the non-linear solver. Since no possible direction was found by the non-linear solver, this split option is discarded.

In order to reproduce the planarization mechanism correctly, the creation of a [0 1 0](1 0 1) segment is necessary. The desired split option is shown in Figure (4.17). The list of split options tested at the beginning of the formation of the helical turn is reproduced in Table (4.5).

#0 1 0 0 - - 8.1623 #1 2 1,1 1 [1 1 1] (1 0 1) 15.1956 #2 2 1,1 1 [1 1 1] (1 0 1) 8.5429 #3 2 1,1 1 [0 1 0] (1 0 0) solver did not converge
As can be seen from the Table (4.5), the split option of type [0 1 0](1 0 1) is not proposed in the list of configurations to be tested. Two virtual nodes connected by one virtual junction are needed for such solution. According to the analysis of split node algorithm, the first virtual node has one degree of freedom, while the other has two degrees of freedom. If one of the created virtual nodes has two dof, the direction of the virtual junction connected to such node cannot be determined from a crystallographic point of view. This direction is then determined with the non-linear solver. In this case, no possible direction was found by the non-linear solver. A possible justification is that the initial guess of the junction direction is far from the right configuration. Such divergence occurs despite the number of performed iterations. As a result, this split option is eliminated.

Results of the interaction of a screw dislocation with a loop of 127 SIAs are similar to that with 271 SIAs. A planarization mechanism should be observed in conformity with MD simulations. However, the formation of a helical turn occurs in the DD simulations. The interaction mechanism we reproduce with DD simulations is shown in Figure (4.18). In case of the interaction with 127 SIAs loop in DD simulations (cf. Figure 4.16), the [1 1 1](1 0 1) segment crosses the simulation boundaries and bows out until the critical stress. At such stress, the [1 1 1](0 1 1) segment is unpinned from the helical turn.

The difference we found between the reaction in 127 and 271 is attributed to the size of the loop. The stress field of the loops evolves rapidly with its size (≈ 1/r 3 ). The force needed to unpin the screw dislocation is therefore strongly affected by the size of the loops.

Finally, it must be noted that the planarization reaction could not be achieved in our DD simulations due to the absence of cross-slip mechanism rule in the code. Another important observation from Figure (4.18) is that the liberated dislocation glides in the (0 1 1) plane which is an anti-twinning direction as discussed earlier in section (1.4.3). The dislocation is supposed to glide in the twinning plane (1 0 1) as in Figure (4.16). Such mechanism is not implemented yet in NUMODIS, it is further discussed in section (4.3.2).

In the following analysis, the planarization mechanism is therefore manually investigated. In other words, the formation of new junction of the corresponding slip systems is preset in the entry files. The resulting configuration is shown in Figure (4.19).

Reactions ending with a helical turn formation

In chapter (1), it was reported that a helical turn is expected during the interaction of screw dislocations with small SIA loops. Reaction mechanisms that lead to the formation of helical turns are considered to have a strong contribution to hardening. This is because the dislocation motion is blocked until the helical turn is closed following Hirsch's mechanism.

The first interaction reproduced with DD simulation between a SIA loop and a dislocation is described in Figure (4.15). In this case, involving a small loop of 37 SIAs, a helical turn is formed. The formation of the helical turn starts when the junction [0 1 0] is transformed into a 1/2[1 1 1] segment at stress level equals to 134 MPa. The formation of the first [1 1 1] segment is shown in Figure (4.20). The hexagonal loop is transformed into a helical turn of five segments through progressive creation of [1 1 1] segments as a result of the split node algorithm. This process is similar to an artificial cross-slip without activation energy, for this reason it is called "easy cross-slip" in the following. At 200 MPa, the helical turn formed after the complete transformation of the hexagonal loop is shown in Figure (4.21). A hexagonal loop is finally left behind with the same Burgers vector of the original dislocation. The resolved CSS in the primary glide plane equal 750 MPa, while it is 360 MPa in the cross-slip (secondary) plane.

The critical resolved shear stress in the (0 1 1) plane is 360 MPa. Here, it must be noted that this cross-slip (secondary) plane where the screw dislocation glides after being unpinned is different from the one observed in MD simulations. This is because in the DD simulations, Hirsch's mechanism (closure of the helical turn) is observed. This mechanism simply explains why the critical shear stress in the primary plane calculated in the DD is 747 MPa, which is 126% higher than the reported value in MD (320 MPa). A much better agreement could be observed if the unpinned dislocation was in the primary plane.

In Figure (4.22), we show that the helical turn extended along the simulation box and crossed the boundary from one side to the other. As a result, the original dislocation disappeared. This disappearance of the original dislocation is not in agreement with the MD results. To elucidate such difference several calculations have been made. Here, the maximum length of the helical turn we observed in the simulation is 400 Å. The physical node at the intersection of the intersection between the first segment of the helical turn and the dislocation has a high velocity. It is equal to 2 × 10 4 m/s, which is considered as a supersonic velocity (speed of sound in iron equal 0.5 × 10 4 m/s [4]). The reason behind this huge velocity is explained in section (2.1.3). In brief, the segment Based on the simulation results of the interaction of a screw dislocation and a SIA loop, it was proposed that a specific viscosity of the physical nodes could exist and should be included in the DD simulations. When two segments connected to a physical node become nearly collinear, the viscosity matrix tends to be ill conditioned. This problem leads to a rapid propagation of the physical node when a helical turn is formed by the Hirsch's mechanism [5] and to the complete elimination of the original screw dislocation in the helical turn. In order to decrease the velocity of these fast moving nodes and to give more time to the dislocation to bend and close the loop again, we tested the possibility to add a specific viscosity on the physical nodes in the direction parallel to the dislocation line.

Such proposition was made because the physical nodes at the intersection of two quasi-collinear planes are the source of an ill-conditioned viscosity matrix KV = F . Such nodes tend to have a high velocity and propagate only in the common direction between the two planes. This configuration implies the need for an additional anisotropic viscosity term for the nodes connected to two segments in different glide plane. The corresponding viscosity model for such configuration is the model 3 presented in section ( 2 This time, the maximum length of the helical turn calculated in the DD simulations was 62 Å (see Figure 4.25c). Due to the bowing of the dislocation arms, the summation of forces exerted on the helical turn led to its closure, thus completing Hirsch's mechanism. The critical shear stress to overcome this obstacle and to complete the process is decreased to 462 MPa (including the flow stress of 22.5 MPa for the screw dislocation without any obstacle). Such value is still larger than the value obtained with MD (320 MPa (44%)). Despite the large difference in the critical shear stress between MD and DD, this result is considered satisfactory. First, the reaction mechanism we reproduce is in relatively good agreement with the observation made in MD. This is mainly, the formation and closure of the helical turn without the disappearance of the original dislocation. Second, our DD simulations for simplicity reason was made by considering a constant linear drag coefficient on the screw dislocations, while in MD screw dislocation mobility is thermally activated and controlled by the double kink mechanism.

In conclusion, the implementation of a viscosity model on physical nodes in the screw direction allows reproducing the interaction mechanism of screw dislocation with small sized SIA loops in the DD simulations in good agreement with MD. Further optimization of the model and the used values should be made in order to reach a more quantitative agreement between MD and DD. The value of the viscosity needed on the physical nodes can be obtained from MD simulations. In section (4.3), MD simulations on the flow stress of a single screw dislocation are performed. Those simulations are similar to those made by Liu et Biner.

Conclusion

In summary, the interaction of a screw dislocation with 1/2[1 1 1] hexagonal loops of different sizes was studied with DD. In the case of small loops of 37 SIAs, a helical turn formation is observed in DD in good agreement with MD simulations. In the case of large loops of 127 and 271 SIAs, a planarization mechanism is expected from MD simulations. However, in these cases we also see the formation of a helical turn. A summary of the critical shear stress found in the MD and DD simulations is given in Figure In the three investigated cases, the helical turn expansion leads to a complete transformation of the original dislocation. The segment bowing out of the helical turn has a different glide plane than the original dislocation. In the case of 127 and 271 SIA, the planarization mechanism is not observed due to the absence of the proper split option or the absence of cross-slip rule in the simulation. The crystallographic directions are correctly proposed, but split option is rejected due to the divergence of the non-linear solver. The selected split option during the DD simulation is the creation of segments of b =[1 1 1] which is energetically most favorable.

Molecular dynamics investigation

In sections (4.2, 4.1), we presented the reaction mechanisms between a screw dislocation and radiation-induced loops using DD simulations. Some reactions observed in MD could not be reproduced using DD simulations because further development of the split node algorithm and the implementation of the cross-slip mechanism are required. However, in the case SIA37 (cf. section 4.2.2) could not be achieved without an extra viscosity term on the physical nodes. Such results open questions regarding the cross-slip and the viscosity of physical nodes which is not addressed in the reference studies [1,3]. Hence, an insight investigation of the MD studies is mandatory.

In this section, MD simulations are made for comparison with the DD simulations presented in the previous sections. For simplicity, the simulation parameters used regarding the interatomic potential, temperature, simulation box, loading, etc., are the same as in Liu et Biner [3].

About the existence of viscosity on physical nodes

In section (4.2), we pointed out the need for an extra viscosity term on physical nodes in the DD simulations. This quantity is needed in the DD simulations to prevent the rapid propagation of the helical ending segments along the dislocation line. Hence, the use of a strong viscosity reducing the mobility of physical nodes appears to be a key parameter to reproduce the same mechanism in MD and DD simulations.

In an attempt to verify the existence of a specific viscosity on physical nodes we compared the flow stress in the case of single edge dislocation and a jogged dislocation. The latter results from the interaction of an edge dislocation with a [1 0 0] loop (see Terentyev et al. [6]). The configuration of the jogged dislocation is shown in Figure (4.27). As can be seen in Figure (4.28), the difference in the flow stress when comparing the mobility of a jogged and a straight dislocation lines is low. This result does not support the possible existence of large viscosity on the physical nodes at high strain rates. The difference in the flow stress we found ranges between (1.9 -13.7 MPa) for velocities between (20 -160 m/s). Such difference is small and can be explained as the outcome of the extra length of dislocation line coming from the presence of a super-jog.

We are now convinced that an extra viscosity term on physical nodes has no significant meaning. Nevertheless, we are still unable to understand the disappearance of the original dislocation in DD simulations after the formation of a helical turn in Liu et Biner. In order to have an insight of interaction mechanism we performed the same simulations using molecular dynamics. This investigation was made using the Lammps code [7]. Results are visualized using the Paraview code [8]. The size of the simulation box is (30 * 16 * 30 nm 3 ), while a constant strain rate of ε = 8 × 10 6 s -1 is applied. The input file for Lammps contains the positions of atoms of the dislocation line and the SIA loop. In order to create a dislocation lines and loop, a specific module in NUMODIS is used to benefit from the manipulation of crystallographic directions already implemented in NUMODIS. The simulated initial configuration is relaxed before the load is applied. Ackland et al. EAM inter-atomic potential is used [9]. The simulation results are shown in Figure (4.29), while the shear stress as function of shear strain is plotted in Figure (4.30).

The critical stress to overcome the loop is 410 MPa. Once the dislocation is liberated, it glides in the (3 2 1) plane as shown in Figure (4.29d). The corresponding glide velocity of the screw dislocation is 43 m/s. Such velocity is higher than the initial glide velocity of a single dislocation in the simulation box at the same applied strain rate. A possible explanation is discussed later in this section. The observed formation of the helical turn is in good agreement with Liu et Biner. However, the helical turn is found to extend rapidly along the dislocation line leading to the elimination of the original dislocation. The [1 1 1](0 1 1) segment on the helical turn bows out under increasing stress. As a result, the original dislocation is found to disappear unlike what was previously reported in Liu et Biner. They only discussed the process of unpinning the dislocation from the loop. However, the disappearance of the original dislocation and the existence of a cross-slip on the unpinned dislocation in the (3 2 1) direction was not mentioned in Liu et Biner.

The cross-slip of the original dislocation was discussed with Liu and Biner through a private communication in an attempt to understand the reaction mechanism with a loop of 37 SIA. Regarding the glide of a single dislocation, without the existence of any obstacle in the simulation box, they observed the following: "For 300 K case, the dislocation travels closer to the surface after some strain level so it must be the cross-slip effect, since it is most likely that the shift is a gradual process rather than an abrupt one" [10]. This confirms our results found using MD simulations of the same interaction. More importantly, such results open questions regarding the cross-slip and the stability of screw dislocations in the glide plane even with the highest resolved shear stress which is discussed in the next section.

In conclusion, in section (4.2.2) we emphasized the need for an extra viscosity term on physical nodes to reproduce the interaction reported by Liu et Biner. Based on the results of the MD simulations made to test such possibility we are now convinced that there is no evidence of an additional viscosity on physical nodes. Rather, DD simulations appear to be in good agreement with MD simulations without including such viscosity.

Cross-slip and twinning mechanisms in α-iron

The next point we investigated is the abnormal cross-slip of a dislocation segment along the observed helical turn we observed in MD simulations. In the simulations we made, a screw dislocation is created in the middle of the simulation box. The size of the simulation volume is (35*35*16 nm). Different shear strain rates are applied that result in a dislocation glide velocity of 5, 10 and 20 m/s. The simulation temperature is set to 300 K. Results of such simulations are given in Figure (4.31). In Figure (4.32b), the flow stress is found to increase by ∼ 53% when the velocity is increased by 300%. At low velocity (5 m/s), the screw dislocation glides in the primary (1 1 0) plane, when at intermediate and high velocity (10 and 20 m/s), the screw dislocation glides in the secondary (3 2 1) plane. Domain et Monnet in [11] studied the temperature dependence of the CRSS for screw dislocations using Mendelev et al. [12] EAM interatomic potential. Although the used potential is slightly different than the one used in our simulations [9], a simple comparison can be still made. The velocity of the screw dislocation in their study is 4 m/s. Such result can be compared to our results for the velocity of 5 m/s at 300 K. At such temperature and dislocation velocity, no cross-slip of the dislocation was mentioned. This is not in agreement with our MD simulations and with the information obtained through private communication with Liu and Biner. The MD simulations we made show that the dislocation preferentially glides in a secondary plane at high velocity. This is in agreement with the study of Chaussidon, in which the same phenomenon was reported. As explained by these authors, dislocation motion can occur both in the primary plane but also in a secondary plane where the RSS is half of the applied stress in the twinning direction.

At low velocity (5 m/s), the screw dislocation glides in (1 1 0) plane but with a behavior similar to the pencil glide. To illustrate this point, the dislocation motion in the In a second step, we reversed the sign of the applied strain rate imposed to the simulation to see if the shifts observed in the MD simulation take place in both dislocation glide direction. The result of such simulation is reproduced in Figure (4.35). Glide shifts are observed again along the positive (1 1 0) direction despite inverting the direction of the applied stress. This behavior can easily be explained by the well-known existence of twinning and anti-twinning directions in iron.

Lastly, twinning and anti-twinning directions in BCC iron were tested with a set of MD simulations. A screw dislocation is created along the [1 1 1] direction. The crystal orientation is adjusted to force the glide of the dislocation in the plane with the maximum resolved shear stress. The results of this MD study are reproduced in Figure (4.36). As can be seen in Figure (4.36), twinning and anti-twinning directions are alternating each 60 • . The flow stress takes the shape of a sinusoidal function of local extrema alternating between twinning and anti-twinning directions respectively. The maximum difference in the flow stress between the twinning and anti-twinning directions is ∼ 63%.

The interaction between a screw dislocation with large hexagonal loops (127 and 271 SIAs) of b =1/2[1 1 1] is investigated in section (4.2.1) using DD simulations. In Figure (4.18), a segment in the (1 3 2) glide plane is created during the reaction. Such segment lies in the anti-twinning direction. The split node algorithm should be modified to privilege split options in the twinning directions since they are more energetically favorable.

Conclusion

In this chapter, a comparison is made between the individual interactions of screw dislocations with SIA loops between DD and MD simulations. The studied SIA loops are square (b = 1 0 0 ) of different orientations and hexagonal (b =[1 1 1]) of different size. The dynamics of each reaction is studied and compared to MD simulations to calibrate our DD code. The aim of this work is to provide a generic algorithm capable of treating different configurations to reproduce the MD results.

The observed reaction mechanisms for such interactions are the restoration of the original loop, planarization and helical turn formation. An insight analysis of the split node algorithm is done to determine whether the configuration found in the DD simulations agrees with MD predictions. In some cases, the correct split configuration could not be reproduced by the split node algorithm due to the inability of the non-linear solver to determine the direction of the junctions involved in the reactions. In addition, we show that the split into three virtual nodes is a key algorithmic element to reproduce some reactions. Finally, we showed that in other cases, the correct configuration could not be obtained due to the absence of cross-slip mechanisms in DD simulations.

In the following a summary of the studied reactions and how to reproduce them in DD in agreement with MD simulations:

• Interaction with 1 0 0 square SIA loops:

-Case A3, the original loop is not restored in DD as in MD simulations. This is due to the absence of a split into three virtual nodes the configurations proposed by the split node algorithm.

-Case B3, the formation of the helical turn could be reproduced in the DD simulations. This result comes from the absence of the creation of [1 1 1] segment by the split node algorithm.

-Case C3, this reaction is not in agreement with MD simulations. This is because of the artificial cross-slip due to the creation of a (1 3 2) segment which should be eliminated.

• Interaction with 1/2[1 1 1] hexagonal SIA loops:

-Cases SIA127 and SIA271, the planarization mechanism could not be achieved because the proper split option cannot be determined from a crystallographic point of view in the split node algorithm. This direction is then determined using the non-linear solver. No possible direction is found due to convergence problem in the non-linear solver.

case SIA37, a helical turn is formed in agreement with MD simulations.

In a nutshell, further development of the split node algorithm is needed in order to reproduce most of the reaction in DD simulations in agreement with MD simulations. Moreover, the implementation of the cross-slip mechanism is inevitable.

In conclusion, the reaction mechanism is affected by the size of the obstacle, orientation, angle of attack of the dislocation and other factors. When performing massive simulations, where each obstacle can be crossed several times by different dislocation lines, the details of the reaction mechanism become important. For this reason, some phenomena such as the clear band formation cannot be studied. numerical or dynamical. The numerical origin can be directly related to the fact that a Newton-Raphson solver is used. In the case where the initial guess is far from the solution, a possible divergence occurs despite the number of iterations. On the other hand, if the problem is dynamical, as explained in Figure (2.14), then no possible physical solution is found due to inherent oscillation of creation and removal of a virtual junction.

• the implementation of the cross-slip and twinning/anti-twinning mechanisms is inevitable. The individual reactions in certain DD simulations cannot be regarded and compared to MD simulations without modelling such dislocation properties.

The reaction mechanism is affected by the size of the obstacle, orientation, angle of attack of the dislocation and other factors. When performing massive simulations, where each obstacle can be crossed several times by different dislocation lines, the details of the reaction mechanism become important. For this reason, some phenomena such as the clear band formation cannot be studied.

Perspectives

A materials multi-scale modeling provides the link between different physical phenomena and processes at a different time and space scales. At the mesoscopic scale, the ultimate objective of this work is to enrich existing crystal plasticity laws using dislocation dynamics. This is done by taking into account the nature of the interaction with the different types of radiation-induced defects. So far in this work, we have treated the loops and the homogeneous precipitates. In the future work, it is essential to provide an assessment of the contribution of cavities and solute atoms to crystal plasticity.

In order to provide a robust crystal plasticity law, massive simulations on representative samples are required at the scale of DD. Therefore, the validation of the individual interactions of different radiation-induced defects with dislocations is important before performing massive simulations.

In this work, we identified the problem of some interactions between dislocations and radiation loops. The next step, is to implement the necessary models, such as cross-slip and twinning-anti twinning mechanisms, in NUMODIS to validate such individual interactions in DD. On the long term, all simulations should be performed using the thermally-activated model of dislocations instead of a linear viscosity law.

Regarding the MD simulations, it is a direct method to provide an insight of the individual interactions between dislocations and radiation-induced defects. It is of a great importance to identify the origin of the abnormal cross-slip observed of screw dislocations at different velocities. One way to identify such problem is simulate the interactions using other interatomic potentials in MD, as Ackland 2003or Mendelev 1997. An important part of this work is centered on the comparison between MD and DD simulations, in an attempt to pass, from the former to the latter, all the elementary features to capture the dislocation-loop interactions. Beyond the technical and missing physical mechanisms that were identified so far, the intrinsic limits of atomistically-informed dislocation dynamics should be questioned. One may in particular wonder whether this (almost) deterministic method can capture nano-scale phenomena, in which the role of thermal fluctuations may be critical. An example of such short-coming is the cross-slip mechanism in faced-centered cubic materials, which involves the constriction of the two Shockley partial dislocations of a screw dislocation through thermal fluctuations. These thermal fluctuations, and more generally stochasticity, are undoubtedly missing in our approach and should be considered in the future.

Adopting a stochastic approach as pioneered by Hiratani et al. [1] is one interesting direction. A general review of this topic is found in [2,3]. Dudarev et al. in [4] and Li et al. in [5] used stochastic DD to study the interaction and diffusion of prismatic dislocation loops in irradiated iron. One promising option is the use of Langevin dynamics (e.g. [6]), where random force terms, related to the dislocation mobility through the fluctuation-dissipation theorem (e.g. [7]), are added to dislocation segments to account for the effective temperature effect.

Finally, when the previous discussed points are verified and validated, they should be implemented in the parallel version of NUMODIS. This is essential to increase the computational power and facilitate performing massive simulations.

1

  Schematic diagram of the major components of a nuclear power plant of type 2 Irradiation effect on the stress-strain curve and the DBT temperature of iron 3 Time-length scale in the materials multi-scale modeling from atomistic to macroscopic scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Schematic diagram of the reactor pressure vessel and the internal components, including fuel and control rods . . . . . . . . . . . . . . . . . . . . . 1.2 Schematic diagram of the nuclear island composed of a RPV connected to four steam generators, four coolant pumps and one pressurizer . . . . . . . 1.3 Continuous cooling transformation diagram of the RPV . . . . . . . . . . . 1.4 Bainitic microstructure of 16MND5 RPV steel observed with optical microscope and SEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 SEM observation of carbides in the RPV formed by a chemical segregation during the solidification process showing austenitic intra-granular sites and at inter-granular site in the ferrite . . . . . . . . . . . . . . . . . . . . . . . 1.6 Stress-strain curve for 16MND5 steel at different temperatures under a strain rate of 5 × 10 -4 s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Ductile and brittle fracture morphology of RPV steel . . . . . . . . . . . . 1.8 Evolution of the ductile-brittle Charpy transition curve in irradiated and unirradiated states of RPV steel under a fluence of Φ = 4.65 × 10 19 n/cm 2 1.9 Irradiation effect on the yield stress as function of strain of pure iron at T = 320 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Schematic diagram of the formation process of edge and screw dislocations in crystalline materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 A screw dislocation part in a curved dislocation undergoes a cross-slip process to a secondary plane . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Arrangement of atoms in twinning / anti-twinning in BCC metals . . . . . 1.13 Principle crystallographic orientations in crystals. The twinning and antitwinning directions of each plane are indicated as (T) or (AT) . . . . . . . 1.14 Microstructural evolution of the number density as function of neutron dose of different types of radiation-induced defects for several alloy compositions using different characterization techniques . . . . . . . . . . . . . . . . . . 1.15 Microstructural evolution of radiation-induced defects for different alloy compositions using different characterization techniques at dose = 0.1 dpa . 1.16 The obstacle strength of the different types of defects in the investigated alloys and the hardening contribution of each defect in function of the hardening measured by tensile tests . . . . . . . . . . . . . . . . . . . . . . 1.17 Concentration profile of Cu precipitates with distance using atomic probe tomography of binary iron-copper alloy (Fe-1.4%Cu) irradiated with neutrons at 288 • C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.18 MD results of the strengthening due to Cu precipitates as function of their diameter and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.19 MD result of the temperature influence on the critical shear stress of an v List of figures edge dislocation to overcome a periodic array of voids of 2 nm size . . . . . 1.20 Interaction of an edge dislocation with a periodic array of voids in α-Fe using MD simulations based on Ackland 2004 interatomic potential . . . . 1.21 CRSS as function of temperatures of the interaction between edge dislocation with voids of 1 or 2 nm diameter for different interatomic potentials using MD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.22 Dependence of the critical shear stress, τ c using static MD simulations at T = 0 K on void's diameter . . . . . . . . . . . . . . . . . . . . . . . . . . 1.23 Solute atoms distribution around an edge dislocation dipole. The chemical composition of the alloy is Fe-0.7%Ni-1.4%Mn. The segregation temperature is 300 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.24 TEM images of the microstructure of RPV steel at different irradiation doses. SIA loops are only observed for high irradiation dos > 0.2 dpa . . . 1.25 Evolution of density and diameter of dislocation loops as function of the dose when irradiated at T = 340 K . . . . . . . . . . . . . . . . . . . . . . 1.26 Schematic representation of the initial configuration of the interaction between an edge dislocation and a hexagonal loop . . . . . . . . . . . . . . 1.27 Different stages of static MD simulation at T = 0 K between edge dislocation and a hexagonal loop of 37 SIAs. A glissile superjog is formed as a result of this interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.28 Different stages of dynamic simulation at T = 300 K between edge dislocation and a hexagonal loop of 331 SIAs . . . . . . . . . . . . . . . . . . . . . . . 1.29 Temperature effect on the critical shear stress for static and dynamic simulations for loops of different sizes . . . . . . . . . . . . . . . . . . . . . 1.30 Configuration C2 of MD simulation result of the interaction between edge dislocation and a square SIA loop introduced in the middle of the glide plane. A double superjog is formed in this mechanism . . . . . . . . . . . . 1.31 Dislocation dynamics results of the C2 configuration . . . . . . . . . . . . . 1.32 Configuration C4U of MD simulation result of the interaction between edge dislocation and a square SIA loop . . . . . . . . . . . . . . . . . . . . . . . 1.33 Dislocation dynamics results of the C4U configuration . . . . . . . . . . . . 1.34 Comparison of the stress-strain curve for MD and DD results for C2, C4U and C5 configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.35 Geometrical configuration of the simulation box that contains a screw dislocation and SIA loops of different orientations . . . . . . . . . . . . . . 1.36 Schematic representation of the first step of the interactions between a screw dislocation and square radiation-induced loops of 128 SIAs in α-iron at different temperatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.37 Snapshots of reaction A3 between a screw dislocation and a square SIA loop of at 300 K using MD simulations. The original loop is restored at the end of this reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.38 Interactions steps between a screw dislocation and a SIA loop at 300 K, the original loop is absorbed as helical turn and restored when the dislocation breaks away . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.39 Normalized critical stress as function of temperature for loops of different orientations and number of SIAs . . . . . . . . . . . . . . . . . . . . . . . . 1.40 Interaction mechanism between screw dislocation and loop of 37 SIAs using vi MD. The original loop is transformed into a helical turn, the dislocation is liberated as in Hirsch's mechanism . . . . . . . . . . . . . . . . . . . . . . 1.41 Interaction mechanism of screw dislocation with a loop of 127 or 271 using MD. A planarization mechanism is observed where the loop is transformed into a faulted loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.42 Schematic representation of the planarization interaction mechanism . . . . 2.1 Discretization schemes of dislocation lines in 3D DD codes . . . . . . . . . 2.2 Scheme of the main loop in the NUMODIS DD code . . . . . . . . . . . . 2.3 Description of the nodal discretization of dislocation lines in NUMODIS . . 2.4 Forces acting on a dislocation according to the line tension model . . . . . 2.5 Schematic representation of the discretization of topological operations on dislocation segments . . . .

3. 1

 1 Schematic representation of the implementation of the minimum image convention concept in NUMODIS where different number of layers can be chosen for nodal force calculations . . . . . . . . . . . . . . . . . . . . . . . 3.2 Formation of a double-kink on a screw dislocation using MD simulation . . 3.3 Thermally activated glide algorithm implemented in NUMODIS . . . . . . 3.4 Schematic representation of two different cases treated by the screw segment detection algorithm (a) curved dislocation line (b) continuous dislocation loop 3.5 Schematic representation of the configuration of a zero length screw dislocation segment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Simulation result of the propagation of Frank-Read source of using NU-MODIS based on a thermally activated model . . . . . . . . . . . . . . . . 3.7 Velocity as function of temperature for a screw dislocation that glides in Fe crystal calculated by NUMODIS compared to the theoretical model . . . . 3.8 Illustration of the DD Simulation volume containing a dislocation and a spherical obstacle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Result of the Orowan stress calculation with screw and edge dislocations in reduced units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Variation of the Orowan stress as function of the Poisson ratio for screw and edge dislocations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Variation of the Orowan stress as function of the core energy parameter ζ core 3.12 Variation of the normalized Orowan stress as function of core width parameter a 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Normalized Orowan stress vs the harmonic mean of inter-obstacle length and obstacle diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Comparison between simulation results and Equation (3.13) prediction for the 1000 combinations of parameters taken from Table (3.1) and tested in the present study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.15 Comparison between critical Orowan stress as function of the inter-obstacle distance L for obstacle of D = 2 nm found using MD and BKS model and DD model in Equation (3.13) . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Illustration of the DD simulation volume used to reproduce Terentyev et al. simulations between a screw dislocation and a square SIA loop . . . . . . . 4.2 Visualization of the first interaction between a screw dislocation and 128 SIAs loop (case A3 in Terentyev et al.) as reproduced in DD simulations . 4.3 Snapshots of the interaction between a screw dislocation and 128 SIAs loop of (case A3 in Terentyev et al.) using DD simulations . . . . . . . . . . . . 4.4 Visualization of the split node reaction selected in the case (A3) when 3 nodes reactions are considered in the DD simulations . . . . . . . . . . . . 4.5 Visualization of the successful creation of the first half of the reaction observed in the MD simulation with the DD simulations thanks to the split option into three nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Snapshots of the interaction of a screw dislocation in case (A3) when the cross-slip is manually produced . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Snapshots of the interaction between a screw dislocation and a 128 SIAs loop (case B3 in Terentyev et al.) as regarded in DD simulations . . . . . . 4.8 Snapshots of the interaction between a screw dislocation and a 128 SIAs loop (case B3 in Terentyev et al.) as reproduced with DD simulations . . . 4.9 Snapshots of the initial configuration of the interaction between a screw dislocation and a 128 SIAs loop (case C3 in Terentyev et al.) as observed in DD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Snapshots of the interaction between a screw dislocation and a 128 SIAs loop (case C3 in Terentyev et al.) as regarded in DD simulations . . . . . . 4.11 Visualization of the node to be pinned in order to prevent the creation of a split option in the (1 3 2) plane in case C3 . . . . . . . . . . . . . . . . . . 4.12 Snapshots of the critical and final configurations of the interaction between a screw dislocation and a 128 SIAs loop case C3 . . . . . . . . . . . . . . . viii 4.13 Comparison of the critical shear stress between MD and DD of the interaction of a screw dislocation with 1 0 0 square loops of 128 SIAs . . . . . . 4.14 Illustration of the DD Simulation volume for the interaction between a screw dislocation and hexagonal loop as in Liu et al. . . . . . . . . . . . . . 4.15 Snapshots of the formation of the first junction as a result of the interaction between a screw dislocation and a hexagonal SIA loop . . . . . . . . . . . . 4.16 Snapshots of the interaction between a screw dislocation and 127 SIAs loop as observed in DD simulations in the current version of the code . . . . . . 4.17 Schematic representation of the node to split at the time step where split node algorithm proposes different split options . . . . . . . . . . . . . . . . 4.18 Snapshots of the interaction between a screw dislocation and 271 SIAs loop 4.19 DD simulation of the planarization mechanism manually produced to show the outcome of MD simulation for this reaction (a) side view, (b) top view 4.20 Visualization of the transformation of the hexagonal loop into a helical turn after the appearance of the first 1 1 1 segment . . . . . . . . . . . . . . . 4.21 The 37 SIA loop is completely transformed into a helical turn of segments with burgers vectors 1 1 1 segments . . . . . . . . . . . . . . . . . . . . . 4.22 Visualization of the configuration after the complete disappearance of the original dislocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Visualization of the final stage of the reaction between a screw dislocation and a 37 SIA loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 Snapshots of the interaction between a screw dislocation and 37 SIAs loop in a large simulation box . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Interaction between a screw dislocation and 37 SIAs loop when the viscosity on physical nodes is increased to 30 MPa.ns in the screw direction . . . . . 4.26 Comparison of the critical shear stress between MD and DD simulations for the interaction of a screw dislocation with a hexagonal loops of 37, 127 and 271 SIAs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Configuration of the super-jog on an edge dislocation as a result of the interaction with a square loop . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Flow stress as function of initial dislocation velocity in the case of a straight and jogged edge dislocations . . . . . . . . . . . . . . . . . . . . . . . . . . 4.29 Visualization of the interaction between a screw dislocation with a hexagonal loop with the same simulation conditions as in Liu et Biner . . . . . . . . . 4.30 Stress-strain curve of the interaction between a screw dislocation with a loop of 37 SIAs using our MD simulations with Ackland 2004 potential . . 4.31 Glide path of a screw dislocation in iron at different velocities at 300 K using MD simulations with Ackland et al. potential . . . . . . . . . . . . . 4.32 Glide of a screw dislocation in iron at 300 K at different velocities in MD simulations using Ackland 2004 potential . . . . . . . . . . . . . . . . . . . 4.33 Glide of a (z-axis) screw dislocation as function of time . . . . . . . . . . . 4.34 Schematic representation of a BCC unit cell of iron with lattice parameter a 0 equal 0.2856 nm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.35 Glide of a (z-axis) screw dislocation as function of time . . . . . . . . . . . 4.36 The flow stress of a screw dislocation for different glide planes in BCC iron.

1. 1

 1 Comparison between the different types of steels used in nuclear reactors regarding the crystallographic structure, mechanical properties and corrosion resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Chemical composition of RPV in French nuclear reactors in weight percent as specified in the French code of design and construction (RCC-M) . . . . 1.3 Reaction summary between a screw dislocation and square radiation-induced loops of 128 SIAs in α-iron at different temperatures . . . . . . . . . . . . 1.4 Critical hardening stress and strength coefficient α for different sizes of SIA loops at T = 300 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Energy and geometrical simulation parameters we explored to calculate the Orowan stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 List of the adjusted parameters in DD simulations to reproduce the MD results for iron at 300 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Identification of mobility law parameters in iron at 300 K. The drag coefficients and the friction stress of different glide systems are indicated . . . . 4.1 Split options of node (A) in Figure (4.2b) at the time step after the collision between a screw dislocation and a square loop . . . . . . . . . . . . . . . . 4.2 Split options of node (A) shown in Figure (4.2b) after the collision between a screw dislocation and a loop . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Split options at the time step after the collision between a screw dislocation with a square loop (case B3) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Split options tested at the collision between a screw dislocation with a square loop (case C3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Split options tested at node (1) in Figure (4.17) at the time step before the beginning of the formation of the helical turn for the reaction of a screw dislocation with a hexagonal loop of 127 SIAs . . . . . . . . . . . . . . . .

  Figure (2a) illustrates the effect of irradiation on the tensile stress of polycrystalline Iron. Figure (2b) shows the irradiation effect on the ductile to brittle transition (DBT) temperature.

Figure 2 -

 2 Figure 2 -(a) Irradiation effect on stress-strain curve of polycrystalline Fe [7].(b) Irradiation effect on the ductile to brittle transition temperature of iron[8]. One of the major design basis accidents in a nuclear reactor is the loss of coolant accident (LOCA) from the primary cooling circuit. When the DBT temperature decreases due to irradiation, in case of a LOCA for example, the temperature of the RPV might drop to values below the DBT temperature. It leads therefore to the initiation of a brittle fracture in the RPV. Hence, it is important to provide the end-users (operators) and decision makers with accurate models that describes the microstructural evolution of the RPV. Such models are used to predict the residual lifetime of the existing nuclear reactors to ensure safe operation in conformity with national and international regulations.

Figure 3 -

 3 Figure 3 -Time-length scale in the materials multi-scale modeling from atomistic to macroscopic scale. The upper half shows the experimental techniques, while the lower half represents different simulation techniques (after SOTERIA) [9]

1 . 2 1 . 3 1 . 4

 121314 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Steel in the nuclear industry . . . . . . . . . . . . . . . . . . . . . Generalities on the RPV steel . . . . . . . . . . . . . . . . . . . . 1.2.1 RPV as a primary component in nuclear reactors . . . . . . . 1.2.2 Chemical composition and heat treatment . . . . . . . . . . . 1.2.3 Microstructure of the RPV . . . . . . . . . . . . . . . . . . . . Mechanical behavior of RPV . . . . . . . . . . . . . . . . . . . . . Dislocations and plastic deformation . . . . . . . . . . . . . . . . 1.4.1 Dislocations in general . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Cross-slip mechanism . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Twinning and anti-twinning directions . . . . . . . . . . . . . 1.5 Radiation-induced defects . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Copper-rich precipitates . . . . . . . . . . . . . . . . . . . . . 1.5.2 Cavities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.3 Solute Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.4 Radiation-induced loops . . . . . . . . . . . . . . . . . . . . . 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -Schematic diagram of the reactor pressure vessel and the internal components including fuel and control rods (after Kok 2009) [5].

Figure 1 . 2 -

 12 Figure 1.2 -Schematic diagram of the nuclear island. The reactor pressure vessel is connected to four steam generators. A reactor coolant pump is connected to each steam generator, while only one pressurizer is connected to the hot leg of one steam generator (after Kok 2009) [5].

Figure 1 . 3 -

 13 Figure 1.3 -Continuous cooling transformation diagram of the RPV. Based on the cooling rate and the initial temperature of the ingot, three microstructures are identified in this diagram, martensite, bainite and ferrite-bainite. The blue arrows indicate the different cooling rates of the inner and outer sides of the RPV due to its considerable thickness (after Raoul 1999) [12].

Figure 1 . 4 -

 14 Figure 1.4 -Bainitic microstructure of 16MND5 RPV steel. Observed with (a) optical microscope, (b) SEM at lower scale. The dark constituent is a bundle of bainitic ferrite laths, while the lighter phase is made of carbides (after Hausild 2002) [14].

Figure 1 . 5 -

 15 Figure 1.5 -SEM observation of carbides in the RPV formed by a chemical segregation during the solidification process. (1) at austenitic intra-granular site (2) at inter-granular site in the ferrite (after Carassou 2000) [16].

Figure 1 . 6 -

 16 Figure 1.6 -Stress-strain curve for 16MND5 steel at different temperatures under a strain rate of 5 × 10 -4 s -1 (after Libert 2007) [18].

Figure 1 . 7 -

 17 Figure 1.7 -Fracture morphology of RPV steel (a) brittle fracture (after Bouchet et al. 2005) [21], (b) ductile fracture (after Tanguy et al. 2008) [22].

Figure 1 . 8 -

 18 Figure 1.8 -Evolution of the ductile-brittle Charpy transition curve in irradiated and unirradiated states of RPV steel. Results for irradiated state are done under a fluence Φ = 4.65 × 10 19 n/cm 2 . ∆T 56J is the reference temperature for Charpy energy as specified by the French nuclear authority (after Bouchet et al. 2005) [21].

Figure 1 . 9 -

 19 Figure 1.9 -Irradiation effect on the yield stress as function of strain of pure iron at T = 320 K (after Victoria et al. 2000) [23].

Figure 1 .

 1 Figure 1.10 -(a) Model of a simple cubic lattice, the atoms are represented by filled circles, (b) positive edge dislocation DC formed by inserting an extra half-plane of atoms in ABCD, (c) left-handed screw dislocation DC formed by displacing the faces ABCD relative to each other in direction AB, (d) spiral of atoms adjacent to the line DC in (c) (after Hull et Bacon 2001) [28].

Figure 1 .

 1 Figure 1.11 -A screw dislocation part in a curved dislocation undergoes a cross-slip process to a secondary plane (after Hull et Bacon 2001) [28].

Figure 1 . 13 -

 113 Figure 1.13 -Principle crystallographic orientations in crystals. The twinning and anti-twinning directions of each plane are indicated as (T) or (AT) respectively (after Chaussidon 2007) [32].

  (see Figure 1.14) for several chemical compositions.

Figure 1 . 14 -

 114 Figure 1.14 -Microstructural evolution of the number density as function of neutron dose of different types of radiation-induced defects for several alloy compositions using different characterization techniques (after Meslin et al. 2010) [36].

Figure 1 . 15 -

 115 Figure 1.15 -Microstructural evolution of radiation-induced defects for different alloy compositions using different characterization techniques at dose = 0.1 dpa. (a) obstacle size, (b) Defect concentration (after Lambrecht et al. 2010) [37].

Figure 1 .

 1 Figure 1.16 -(a) The obstacle strength from Equation (1.2) of the different types of defects in the investigated alloys. (b) the hardening contribution of each defect in function of the hardening measured by tensile tests (after Lambrecht et al. 2010) [37].

Figure 1 . 17 -

 117 Figure 1.17 -Atomic probe tomography of binary iron-copper alloy (Fe-1.4%Cu) irradiated with neutrons at 288 • C (a) concentration profile of Cu precipitates with distance through the precipitate (b) cluster of Cu atoms (after Pareige et Auger 1994) [44].

Figure 1 .

 1 Figure 1.18 -MD results of the strengthening due to Cu precipitates (a) CSS as function of Cu precipitate diameter for different lengths of the dislocation, (b) CSS as function of temperature for different sizes of the Copper precipitates (after Bacon et al. 2004) [47].

Figure 1 .

 1 Figure 1.19 -MD result of the temperature influence on the normalized critical shear stress of an edge dislocation to overcome a periodic array of voids of 2 nm size (after Osetsky et al. 2005) [59].

Figure 1 . 20 -

 120 Figure 1.20 -Interaction of an edge dislocation with a periodic array of voids in α-Fe using quasi-static MD simulations based on Ackland 2004 interatomic potential at T = 1 K. Black and green lines represent the critical configuration of dislocation line shape in the glide plane for voids with D = 2 nm and (L + D) = 21 or 41 nm respectively (after Terentyev et al. 2008) [55].

Figure 1 .

 1 Figure 1.21 -CRSS as function of temperatures of the interaction between edge dislocation with voids of 1 or 2 nm diameter for different interatomic potentials as in [60, 61] using static and dynamic MD simulations (after Terentyev et al. 2008) [55].

Figure 1 . 22 -

 122 Figure 1.22 -Dependence of the critical shear stress, τ c using static MD simulations at T = 0 K on void's diameter and dislocation length (after Osetsky et al. 2003) [40].

Figure 1 . 23 -

 123 Figure 1.23 -Solute atoms distribution around an edge dislocation dipole. The chemical composition of the alloy is Fe-0.7%Ni-1.4%Mn. The segregation temperature is 300 K. The dislocation dipole is presented in black lines, while colors of solute atoms are: Cu in green, Mn in violet and Ni in red (after Pascuet et al. 2017) [66].

Figure 1 .

 1 Figure 1.24 -TEM images of the microstructure of RPV steel at different irradiation doses. SIA loops are only observed at high irradiation dose > 0.2 dpa (after Zinkle et Singh 2006) [50].

Figure 1 . 25 -

 125 Figure 1.25 -Experimental observation of the evolution of density and diameter of dislocation loops as function of displacement damage (dpa) when irradiated at T = 340 K (after Zinkle et Singh 2006) [50].

Figure 1 . 26 -

 126 Figure 1.26 -Schematic representation of the initial configuration of the interaction between a b =1/2[1 1 1] edge dislocation and a b =1/2[1 1 1] hexagonal loop used in MD simulations. The applied shear direction is indicated with arrows on the top and bottom of the simulation box (after Bacon et al. 2006) [79].

Figure 1 . 27 -

 127 Figure 1.27 -Different stages of static MD simulation at T = 0 K between an edge dislocation and 1/2[1 1 1] hexagonal loop of 37 SIAs. A glissile superjog is formed as a result of this interaction (after Bacon et al. 2006) [79].

Figure 1 . 28 -

 128 Figure 1.28 -Different stages of dynamic simulation at T = 300 K between edge dislocation and a 1/2[1 1 1] hexagonal loop of 331 SIAs (after Bacon et al. 2006) [79].

Figure 1 . 29 -

 129 Figure 1.29 -Calculation with molecular static and molecular dynamics of the temperature effect on the critical stress defining the interaction strength between a screw dislocation and loops of 37 and 331 SIAs under different strain rates (indicated by sr) (after Bacon et al. 2006) [79].

  30 -1.33) are snapshots of two different cases studied using MD[START_REF] Terentyev | Simulation of the interaction between an edge dislocation and a < 100> interstitial dislocation loop in α-iron[END_REF] and compared directly to DD[START_REF] Shi | Interaction of <100> dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF]. A good agreement was found between the two simulation methods.

Figure 1 .

 1 Figure 1.30 -Configuration C2 of MD simulation result of the interaction between edge dislocation and a [1 0 0] SIA loop introduced in the middle of the glide plane. A double superjog is formed in this mechanism (after Terentyev et al. 2008) [80].

Figure 1 .

 1 Figure 1.31 -C2 configuration in [80] reproduced using dislocation dynamics simulations (after Shi et al. 2015) [81].

Figure 1 . 32 -

 132 Figure 1.32 -Configuration C4U of MD simulation result of the interaction between edge dislocation and a [100] SIA loop. The lower corner of the loop is placed in the glide plane. Two loops are left behind of b =1/2[1 1 1] and 1/2[1 1 1] connected by a shared segment of b =[1 0 0] (after Terentyev et al. 2008) [80].

Figure 1 .

 1 Figure 1.33 -C4U configuration in [80] reproduced using dislocation dynamics simulations (after Shi et al. 2015) [81].

Figure 1 . 34 -

 134 Figure 1.34 -Comparison of the stress-strain curve for MD and DD results for C2, C4U and C5 configurations (after Shi et al. 2015) [81].

Figure 1 .

 1 Figure 1.35 -Geometrical configuration of the simulation box that contains a screw dislocation and SIA loops of different orientations. Arrows in top and bottom of simulation box indicate the direction of the applied strain rate (after Terentyev et al. 2010) [82].

Figure 1 . 36 -

 136 Figure 1.36 -Schematic representation of the first step of the interactions between a 1/2[1 1 1] screw dislocation and square 1 0 0 radiation-induced loops of 128 SIAs in α-iron at different temperatures. (a) reaction C1, (b) reaction A3, (c) reaction B3 [82].

Figure 1 .

 1 Figure 1.37 -Snapshots of reaction A3 between a 1/2[111] screw dislocation and a square SIA loop of b = [001] at 300 K using MD simulations. The original loop is restored at the end of this reaction (after Terentyev et al. 2010) [82].

Figure 1 .

 1 Figure 1.38 -Interactions steps between a screw dislocation and a SIA loop of b = [0 1 0] at 300 K using MD simulations. The original loop is absorbed as helical turn and restored as 1/2[1 1 1] loop when the dislocation breaks away (after Terentyev et al. 2010) [82].

Figure 1 .

 1 Figure 1.39 -A summary of the normalized critical stress as function of temperature for loops of different orientations and number of SIAs (after Terentyev et al. 2010) [82].

Figure 1 . 40 -

 140 Figure 1.40 -Interaction mechanism between screw dislocation and loop of 37 SIAs using MD simulations. The original loop is transformed into a helical turn, the dislocation is liberated as in Hirsch's mechanism (after Liu et Biner 2008) [83].

Figure 1 .

 1 Figure 1.41 -Interaction mechanism of screw dislocation with a loop of 127 or 271 using MD simulations. A planarization mechanism is observed where the loop is transformed into a [1 0 0] faulted loop (after Liu et Biner 2008) [83].

Figure 1 .

 1 Figure 1.42 -Schematic representation of the planarization interaction mechanism (after Liu et Biner 2008) [83].

  In general, DD codes can be classified into two main categories: lattice-based or nodal codes. This classification is based on the discretization scheme of the dislocation lines. Different approaches of discretization in DD codes are shown in Figure2.1).

Figure 2 . 1 -

 21 Figure 2.1 -Discretization schemes of dislocation line in 3D DD codes. (a) Lattice-based discretization of type Edge-Screw, (b) Lattice-based discretization of type Edge-Mix-Screw, (c) Nodal discretization.

  (2.2). The steps listed in this figure are repeated at each time step until the end of the simulation.

Figure 2 . 2 -

 22 Figure 2.2 -Scheme of the main loop in the NUMODIS DD code at each time step during the simulation.

Figure 2 . 3 -

 23 Figure 2.3 -Nodal discretization of dislocation line in NUMODIS. The Burgers vector b, the tangent unit vector ξ and the normal to the glide plane n are indicated. Physical nodes are red and topological nodes are blue (after Shi 2015) [12].

Figure 2 . 4 -

 24 Figure 2.4 -Forces acting on a dislocation according to the line tension model. ds is an infinitesimal part of the dislocation (a) the applied stress field on the dislocation by Peach-Koehler force, where such force is balanced by the line tension, (b) the line tension is the gradient of the forces exerted on a dislocation segment of length ds (after Dupuy et Fivel 2002) [15].

Figure 2 . 5 -

 25 Figure 2.5 -Schematic representation of the discretization of topological operations on dislocation segments. A long segment is split when a new node is added in the middle. Short segments are merged when the node in between is removed (after Bulatov et Cai 2006) [20].

Figure 2 . 6 -

 26 Figure 2.6 -Schematic representation of the concept of periodic boundary condition. A simulation unit cell is indicated by the dashed square and it has periodic replicas in a two-dimensional space. The interaction of any element in the unit cell is calculated with the nearest neighbor using the minimum image convention (after Frenkel et Smith 2001) [23].

Figure 2 . 7 -

 27 Figure 2.7 -Illustration of the problem of dislocation auto-annihilation induced by PBC in DD simulation (after Madec 2001)[26]. When a cubic (here a square in 2D) volume is used, a dislocation loop expanding in a plane of type [1 1 0] self-annihilates after passing two times the boundaries of the simulated volume. This is because the dislocation line is re-introduced in its initial glide plane.

Figure 2 . 8 -

 28 Figure 2.8 -A susceptible splittable node i (in black color) with n connections.

Figure 2 . 9 -

 29 Figure 2.9 -Illustrative example of the verification of any geometrical incompatibilities for two cases: (a, b) connections of the splittable node are in the same glide plane, (c) connections of the splittable node have different glide planes.

  10, 2.11, 2.12).

Figure 2 . 10 -Figure 2 . 11 -

 210211 Figure 2.10 -Schematic representation of a no split case. In this case, the splittable node in Figure (2.8) is transformed to a virtual node. The four connections are connected to the same virtual node.

Figure 2 . 12 -

 212 Figure 2.12 -Schematic representation of split options for an initial configuration of four connections into three virtual nodes (B,C and D). Segments 5 and 6 (in yellow color) are virtual junctions connected to the virtual nodes.

  2.13). The summation of Burgers vector for nodes B and C is b 1 + b 2 = b B , where b B is the Burgers vector of the created junction B and b 3 + b 4 = b C , where b C is the Burgers vector of the created junction C.

Figure 2 . 13 -

 213 Figure 2.13 -Schematic representation of a split option for a node of five connections, where Burgers vector of the connections and the virtual junctions are indicated: (a) the first set of connections to the first virtual node (b) the second set of connections to the second virtual node.

Figure 2 . 14 -

 214 Figure 2.14 -Schematic representation of an oscillation process between consecutive time steps for an example case of split option. (a) A splittable node at time = t has a tendency to split due to the forces of its connections (in black arrows), (b) the split option into two virtual nodes is chosen with a junction in between, in this case back-stress forces points to the direction opposite to the forces exerted by the connections. The back-stress forces are elevated and this configuration has a tendency to merge the two nodes.

Figure 2 . 15 -

 215 Figure 2.15 -Schematic diagram of a Lomer-Cottrell junction[13]. Connections (1,3) and (2,4) are found in planes 1 and 2 respectively. The virtual nodes (in red color) have a dof equal one, they can move only in one direction. The virtual junction (in yellow color) is therefore a 1D junction.

Figure 2 .

 2 Figure 2.16 -The KV = F matrices divided into sub-matrices of similar characters (splittable, non-splittable nodes) to facilitate the calculation of the Schur complement.

Figure 2 . 18 -

 218 Figure 2.18 -Schematic representation of a split option for a node of three connections.(a) the initial configuration of splittable node of three connections, (b) split option into two virtual nodes, where segments 1 and 2 are connected to a virtual node and segment 3 is connected to the other virtual node. θ is the critical angle between the connection and its Burgers vector. A minimum value is authorized in order to prevent numerical instability.

Figure 3 . 1 -

 31 Figure 3.1 -Schematic representation of the implementation of the minimum image convention concept in NUMODIS. Different number of layers can be chosen for nodal force calculations.

  (3.2).

Figure 3 . 2 -

 32 Figure 3.2 -Formation of a double-kink on a [111] screw dislocation using MD simulation. Screw dislocation glides in the [ 12 1] direction in the ( 101) (after Narayanan et al. 2014) [5].

Figure 3 . 3 -

 33 Figure 3.3 -Thermally activated glide algorithm implemented in NUMODIS.

  At each time step in NUMODIS, screw segments are detected and their length is stored. Lebon showed three different situations that should be treated to determine the length of each screw segment. Two possible configurations are shown in Figure(3.4).

Figure 3 . 4 -

 34 Figure 3.4 -Schematic representation of two different cases treated by the screw segment detection algorithm (a) a curved dislocation line, it contains three screw segments between nodes [1-5], [20-50] and [80-85] (b) continuous dislocation loop, it contains two screw segments between nodes [20-50] and [80-10]. The blue arrow is the Burgers vector direction. (after Lebon 2011) [12].

Figure 3 . 6 -

 36 Figure 3.6 -Simulation result of the propagation of Frank-Read source of using NUMODIS. The Frank-Read source initially has a screw character. The thermally activated glide is based on the double-kink model. Screw dislocation parts (a-b), (c-d)and (e-f) has a low glide velocity compared to other non-screw parts.

  Another test to validate the implementation of the thermally activated glide is to study the effect of temperature. It is shown in Figure (3.7).

Figure 3 . 7 -

 37 Figure 3.7 -Velocity as function of temperature for a screw dislocation that glides in Fe crystal calculated by NUMODIS compared to the theoretical model. The red line represents the velocity as predicted by the double-kink model in [7].

  's mechanism using dislocation dynamics based on the non-singular elastic theory M. Shukeir 1 , L. Dupuy 1 * , B. Devincre 2

Figure 3 . 8 -

 38 Figure 3.8 -Illustration of the DD Simulation volume. The green surface represents the glide plane of the mobile dislocation, while arrows at the top and bottom surfaces indicate the direction of the applied shear stress. The straight dislocation represents the initial dislocation configuration while the curved one is the critical configuration observed before the formation of the Orowan loop.

Figure 3 . 9 -

 39 Figure 3.9 -Result of the Orowan stress calculation with screw and edge dislocations in reduced units. Each mark is a unique combination of the parameters reported in Table(3.1). Continuous lines are the prediction of the BKS model (Equation 3.11).

  Figure 3.9 -Result of the Orowan stress calculation with screw and edge dislocations in reduced units. Each mark is a unique combination of the parameters reported in Table(3.1). Continuous lines are the prediction of the BKS model (Equation 3.11).

  Figure 3.9 -Result of the Orowan stress calculation with screw and edge dislocations in reduced units. Each mark is a unique combination of the parameters reported in Table(3.1). Continuous lines are the prediction of the BKS model (Equation 3.11).

Figure 3 .

 3 Figure 3.10 -Variation of the Orowan stress as function of the Poisson ratio for screw and edge dislocations. Other simulation parameters are L = 1000b (247 nm), D = 100b (24.7 nm), a 0 = 0.4 nm and ζ core = 0.86.

Figure 3 .

 3 Figure 3.11 -Variation of the Orowan stress as function of the core energy parameter ζ core . Other simulation parameters are set to L = 100b (24.7 nm), D = 10b (2.47 nm), a 0 = 0.15 nm and ν = 0.435.

Figure 3 . 12 -

 312 Figure 3.12 -Variation of the normalized Orowan stress as function of core width parameter a 0 . Other simulation parameters are set to L = 1000b (247 nm), D = 100b (24.7 nm), ζ core = 0.86 and ν = 0.435.

Figure 3 .

 3 Figure 3.13 -Normalized Orowan stress vs the harmonic mean of inter-obstacle length and obstacle diameter. The simulation parameters controlling the elastic energy are set to a 0 = 0.45 nm, ζ core = 0.86 and ν = 0.33. The continuous lines are the BKS model prediction with the fit parameter calculated in [30].

Figure 3 . 14 -

 314 Figure 3.14 -Comparison between simulation results and Equation (3.13) prediction for the 1000 combinations of parameters taken from Table (3.1) and tested in the present study.

Lehtinen

  et al. studied the interaction between 1 / 2 [111](1 10) edge dislocation with non-coherent cementite precipitate (Fe 3 C) of different spacing using molecular dynamics. These calculations were produced using the interatomic potential H13 proposed by Henriksson et al. describing the FeCrC system [49]. Parameters that can be directly shared between MD and DD simulations are µ = 75 GPa, ν = 0.379 at T = 750 K, b = 0.2502 nm and obstacle diameter D = 2 nm. Use is then made of Equation (3.13) to identify the missing parameters needed in the DD simulations to correctly adjust the dislocation strain energy and to be in agreement with the prediction of the interatomic potential used in the MD simulations. Such adjustment gives for the dislocation core width parameter a 0 = 0.45 nm and for the dislocation core energy parameter ζ core = 0.51. Both parameters value are physically meaningful and allow running DD simulations to either reproduce quantitatively the MD simulation results in a few seconds or to upscale the MD simulation results. To illustrate this point, a comparison is made in Figure (3.15) between the initial results of Lehtinen et al., the results we obtained with the adjusted DD simulation and the BKS model predictions.

Figure 3 . 15 -

 315 Figure 3.15 -Comparison between critical Orowan stress as function of the inter-obstacle distance L for obstacle of D = 2 nm found using MD and BKS model and DD model in Equation (3.13). DD simulation parameters are a 0 = 0.45 nm, ζ core = 0.51.

Figure 4 . 1 -

 41 Figure 4.1 -Illustration of the DD simulation volume used to reproduce Terentyev et al. simulations [1]. The green surface represents the glide plane of the mobile screw dislocation, while arrows at the top and bottom surfaces indicate the direction of the applied shear stress. The SIA loop is introduced at the center of the simulation box and intersects the glide plane of the mobile screw dislocation. The arrows on the loop indicate the Burgers vector.

  (4.2a). As a result of the applied stress, the dislocation moves in the glide plane to meet the loop at the point (A) as shown in Figure(4.2b).

Figure 4 . 2 -

 42 Figure 4.2 -Visualization of the first interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop of b=[0 0 1] (case A3 in Terentyev et al. [1]) as reproduced in DD simulations. (a) the initial configuration, (b) the dislocation encounters the loop and node A is considered a node to be splitted.

Table 4 . 1 -

 41 Split options of node (A) in Figure (4.2b) at the time step after the collision between the 1/2[1 1 1] screw dislocation with a [0 0 1] loop. * VN: number of virtual nodes, ** VJ: number of virtual junctions, dof are 0,1 or 2 Split option VN* dof VJ** Burgers Plane Dissipation rate (eV/ps)

  (

4 . 3 )

 43

Figure 4 . 3 -

 43 Figure 4.3 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop (b=[0 0 1] case A3 in Terentyev et al. [1]) as reproduced in DD (a) no junction is formed and the dislocation continues gliding, (b) the dislocation approaches the second node to be splitted but no interaction takes place similar to point (a), (c) the original loop is restored and the CSS to pass the loop equals 208 MPa.

Figure 4 .

 4 2b). The outcome of this reaction is shown in Figure (4.4a).

Figure 4 . 4 -

 44 Figure 4.4 -(a) Visualization of the split node reaction selected in the case (A3) when 3 nodes reactions are considered in the DD simulations. Two junctions are created on the [1 1 1](1 0 1) and the [1 1 1](1 0 1) slip systems in agreement with MD simulations. The virtual power dissipation due to this configuration equals 0.1233 eV /ps. (b) the same interaction as observed in MD simulations.

Figure 4 .

 4 Figure 4.4a) glides along the loop segment until it reaches the lower corner of the loop. Then at this step, a segment in the loop is completely changed. As illustrated in Figure (4.5), this segment has a new Burgers vector and a new glide plane.

Figure 4 . 5 -

 45 Figure 4.5 -Visualization of the successful creation of the first half of the reaction observed in the MD simulation with the DD simulations thanks to the split option into three nodes.

Figure 4 . 6 -

 46 Figure 4.6 -Snapshots of the interaction of 1/2[1 1 1](1 1 0) screw dislocation in (A3) when the cross-slip is manually produced. The original loop is restored with a bi-loop in the [1 1 1](1 0 1). The CSS to unpin the mobile dislocation is 530 MPa.

  According to MD, the interaction of a screw dislocation with a loop of b = [0 1 0] (B3) and [1 0 0] (C3) should result in a helical turn formation. In this section, we present the results we obtained for both cases in DD simulations. The initial configuration of case (B3) is shown in Figure (4.7a). A repulsive force exists initially between the dislocation and the loop. As a result of the applied stress, the dislocation moves in its glide plane to meet the loop and forms a [1 1 1](1 1 2) junction as shown in Figure (4.7b).

Figure 4 . 7 -

 47 Figure 4.7 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop of b=[0 1 0] (case B3 in Terentyev et al. [1]) as regarded in DD. A formation of a [1 1 1](1 1 2) junction is observed.

  of a first [1 1 1](1 1 2) segment, several contact reactions are tested with the split node algorithm. Such reactions are equivalent to test easy cross-slip reaction at a specific location along the dislocation loop. The most favorable sequence of the reaction is shown in Figure (4.8).

Figure 4 . 8 -

 48 Figure 4.8 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop (b=[0 1 0] case B3 in Terentyev et al. [1]) as reproduced with DD simulations. No formation of the helical turn is observed in this reaction in DD simulations. The CSS to bypass the loop equals 476 MPa.

  (4.9b).

Figure 4 . 9 -

 49 Figure 4.9 -Snapshots of the initial configuration of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop (b=[1 0 0] case C3 in Terentyev et al. [1]) as observed in DD simulations. No interaction collision takes place between the screw dislocation and the first interacting segment in the loop.

Figure 4 . 10 -

 410 Figure 4.10 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop (b=[1 0 0] case C3 in Terentyev et al. [1]) as regarded in DD simulations.The CSS to overcome the loop equals 750 MPa.

Table 4 . 4 -

 44 Split options tested at the collision between the 1/2[1 1 1] screw dislocation with a [1 0 0] loop (case C3) where the creation of (1 3 2) segment is observed. * VN: number of virtual nodes, ** VJ: number of virtual junctions, dof are 0,1 or 2 Split option VN* dof VJ** Burgers Plane Dissipation rate (eV /ps)

Figure 4 . 11 -

 411 Figure 4.11 -Visualization of the node to be pinned in order to prevent the creation of a split option in the (1 3 2) plane in case C3.

Figure 4 . 12 -

 412 Figure 4.12 -Snapshots of (a) the critical configuration, (b) the final configuration of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 128 SIAs loop b=[1 0 0] case C3. No creation of a (1 3 2) segment is observed due to the pinning of the node to split in Figure (4.11).

Figure 4 . 13 -

 413 Figure 4.13 -Comparison of the critical shear stress between MD and DD of the interaction of a 1/2[1 1 1](1 1 0) screw dislocation with 1 0 0 square loops of 128 SIAs. The loops have different Burgers vector [0 0 1], [0 1 0] and [1 0 0] named A3, B3 and C3 respectively.

  The simulation conditions in DD are taken similar to MD simulations. The size of the simulation box is (30 * 16 * 30) nm3 and the crystallographic orientation of x, y and z axes are parallel to the [1 1 2], [1 1 0] and [1 1 1] directions, respectively. Periodic boundary conditions are applied in the [111] and[ 11 2] directions, while no particular conditions are applied to the boundary surfaces in the [11 0] direction since those surfaces are normal to the mobile dislocation glide direction. A screw dislocation is introduced parallel to the z-axis and a SIA loop is introduced at the center of the simulation box. A constant shear strain rate ε = 8 × 10 6 s -1 is applied. This strain rate corresponds to a steady state dislocation velocity of 30 m/s. The DD simulation box is shown in Figure(4.14).

Figure 4 . 14 -

 414 Figure 4.14 -Illustration of the DD Simulation volume. The green surface represents the glide plane of the mobile screw dislocation, while arrows at the top and bottom surfaces indicate the direction of the applied shear stress. The SIA loop is introduced at the center of the simulation box and intersects the glide plane of the dislocation. The arrows on the loop indicate its Burgers vector.

Figure 4 . 15 -

 415 Figure 4.15 -Snapshots of the formation of the first junction as a result of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and a 1/2[1 1 1] hexagonal SIA loop. The junction belongs to the [0 1 0](1 0 1) slip system and it is a sessile junction.

Figure 4 . 16 -

 416 Figure 4.16 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 127 SIAs loop (b=[1 1 1]) as observed in DD in the current version of the code. A helical turn is formed in the process instead of the planarization of the loop along the glide plane (1 1 0).

Figure 4 . 17 -

 417 Figure 4.17 -Schematic representation of the node to split at the time step where split node algorithm proposes different split options. (a) the initial configuration, (b) the desired split option into 2 virtual nodes, the created junction should be of type [0 1 0](1 0 1).

Figure 4 .

 4 Figure 4.18 -(a -d) Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 271 SIAs loop (b=[1 1 1]). Here, a helical turn is formed instead of the planarization of the loop along the glide plane (1 1 0) as observed in MD simulations. The dislocation in (d) breaks away in the anti-twinning direction as shown in (e), (e) twinning and anti-twinning directions in BCC crystals.

Figure 4 .

 4 Figure 4.19 -DD simulation of the planarization mechanism manually produced to show the outcome of MD simulation for this reaction (a) side view, (b) top view.

Figure 4 . 20 -

 420 Figure 4.20 -Visualization of the transformation of the hexagonal loop into a helical turn after the appearance of the first 1 1 1 segment. The Burgers vectors and the glide plane of each segment are illustrated. The resolved shear stress on the (1 1 0) glide plane is 134 MPa.

Figure 4 . 21 -

 421 Figure 4.21 -The 37 SIA loop is completely transformed into a helical turn of segments with burgers vectors 1 1 1 segments. The resolved shear stress on the principal glide plane is 200 MPa. The resolved CSS in the [1 1 1](0 1 1) is 100 MPa.

Figure 4 . 22 -

 422 Figure 4.22 -Visualization of the configuration after the complete disappearance of the original dislocation. The resolved CSS in the [1 1 1](1 1 0) equal 522 MPa, while the resolved CSS on the extended helical turn segment equal 100 MPa.

Figure 4 . 23 -

 423 Figure 4.23 -Visualization of the final stage of the reaction. The [1 1 1](0 1 1) segment glides in the cross-slip plane, the helical turn is then closed following Hirsch's mechanism.A hexagonal loop is finally left behind with the same Burgers vector of the original dislocation. The resolved CSS in the primary glide plane equal 750 MPa, while it is 360 MPa in the cross-slip (secondary) plane.

First, we

  expanded the length of the helical turn in order to see if the original dislocation does not disappear with larger dimensions of the simulation box. In order to do such test, a 1/2[1 1 1] loop and a screw dislocation are placed in a larger simulation box with dimensions (100, 16, 100) nm 3 . Results of this test are shown in Figure (4.24).

Figure 4 . 24 -

 424 Figure 4.24 -Snapshots of the interaction between 1/2[1 1 1](1 1 0) screw dislocation and 37 SIAs loop (b=[1 1 1]) in a large simulation box. The dimensions of the simulation box are (100, 16, 100) nm 3 . The helical turn is formed and closed in agreement with MD results.

[1 1 1 ]

 1 (0 1 1) in Figure(4.21) is almost aligned with the dislocation arm. Hence, it leads to an ill-conditioned nodal viscosity matrix. When this matrix is inverted to calculate the nodal velocities, the velocity of the node is tremendous.

  .1.3). In this DD simulation, the applied viscosity on physical nodes is set equal to 30 MPa.ns. As expected, the two physical nodes at the extremities of the helical turn faced an extra resistance to their motion along the dislocation line. With this additional rule in the DD simulation, the formation and closure of the helical turn are reproduced as shown in Figure (4.25).

Figure 4 . 25 -

 425 Figure 4.25 -Interaction between 1/2[111](1 10) screw dislocation and 37 SIAs loop (b = [1 11]) when the viscosity on physical nodes is increased to 30 MPa.ns in the screw direction.

  (4.26).

Figure 4 . 26 -

 426 Figure 4.26 -Comparison of the critical shear stress between MD and DD simulations for the interaction of a 1/2[1 1 1](1 1 0) screw dislocation with [1 1 1] hexagonal loops of 37, 127 and 271 SIAs.

Figure 4 . 27 -

 427 Figure 4.27 -Configuration of the super-jog on an edge dislocation as a result of the interaction with a [1 0 0] square loop. The red dots are topological nodes, while black dots are physical nodes by definition.

Figure 4 .

 4 Figure 4.28 -Flow stress as function of initial dislocation velocity in the case of a single straight edge dislocation and a jogged dislocation made after interacting with a square loop. The total length of the dislocation segments is higher in the case of a jogged dislocation.

Figure 4 . 29 -

 429 Figure 4.29 -Visualization of our MD results of the interaction between a 1/2[1 1 1](1 1 0) screw dislocation with a 1/2[1 1 1] hexagonal loop with the same simulation conditions as in Liu et Biner. The green surface represents the glide plane of the original dislocation. The inset figure in (a) represents the twinning/anti-twinning directions in BCC crystal.

Figure 4 .

 4 Figure 4.30 -Stress-strain curve of the interaction between a [1 1 1] screw dislocation with a 1/2[1 1 1] loop of 37 SIAs using our MD simulations. The used interatomic potential is Ackland 2004. The velocity of the screw dislocation is 20 m/s. The red line indicates the averaged shear stress as function of strain.

Figure 4 . 31 -

 431 Figure 4.31 -Snapshots of our MD simulations of the glide path of a 1/2[1 1 1] screw dislocation in iron at different velocities 5, 10 and 20 m/s. The MD simulations are made at 300 K using Ackland et al. EAM interatomic potential. The colored lines indicate possible crystallographic glide planes. (a) the initial configuration of three screw dislocations in the simulation box, (d) screw dislocations at high velocity (10, 20 m/s) glide in the (3 2 1) plane, but in the (1 1 0) plane at low velocity. The inset figure in (a) represents the twinning/anti-twinning directions in BCC crystal.

Figure ( 4 .

 4 Figure(4.31) shows that the simulated screw dislocation glides in different directions (glide planes) at different applied strain rates. The simulated stress-strain curves are plotted in Figure (4.32a) and the average flow stress (plateau stress) as function of the imposed velocity is shown in Figure (4.32b).

Figure 4 . 32 -

 432 Figure 4.32 -Glide of a 1/2[1 1 1] screw dislocation in iron at 300 K at different velocities using Ackland 2004 potential. (a) stress strain diagram, (b) flow stress as function of velocity.

[1 1 0

 1 ] and [1 1 2] directions as function of simulation time is traced in Figure(4.33).

Figure 4 . 33 -

 433 Figure 4.33 -Glide of 1/2[1 1 1] (z-axis) screw dislocation along Y=[1 1 2] (in violet) and Z=(1 1 0) (in green) as function of time. A positive shear strain rate ε = 10 7 s -1 is applied in the xz direction where the MRSS is in the (1 1 0) plane.

Figure 4 . 34 -

 434 Figure 4.34 -Schematic representation of a BCC unit cell. The lattice parameter a 0 equal 0.2856 nm in iron. The arrow indicates the [1 1 0] direction.

Figure 4 . 35 -

 435 Figure 4.35 -Glide of 1/2[1 1 1] (x-axis) screw dislocation along Y=[1 1 2] (in violet) and Z=(1 1 0) (in green) as function of time. A negative shear strain rate ˙ = -10 7 s -1 is applied in the xz direction where the MRSS in in the (1 1 0) plane.

Figure 4 . 36 -

 436 Figure 4.36 -The flow stress of a screw dislocation for different glide planes in BCC iron. Twinning and anti-twinning directions are indicated for each plane. The [1 1 1] direction is perpendicular to the surface of the circle. Each inner circle represents 50 MPa flow stress.

  . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Schematic representation of the concept of periodic boundary condition . . 2.7 Illustration of the problem of dislocation auto-annihilation induced by PBC in DD simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 A susceptible splittable node i (in black color) with n connections . . . . . 2.9 Illustrative example of the verification of any geometrical incompatibilities when the connections of the splittable node have different or similar glide planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Schematic representation of a split into one virtual node where the splittable node is transformed into a virtual node . . . . . . . . . . . . . . . . . . . . 2.11 Schematic representation of split options of initial configuration of four connections into two virtual nodes . . . . . . . . . . . . . . . . . . . . . . . 2.12 Schematic representation of split options for an initial configuration of four connections into three virtual nodes . . . . . . . . . . . . . . . . . . . . . . 2.13 Schematic representation of a split option for a node of five connections, where Burgers vector of the connections and the virtual junctions are indicated 2.14 Schematic representation of an oscillation process between consecutive time steps for an example case of split option . . . . . . . . . . . . . . . . . . . 2.15 Schematic diagram of a Lomer-Cottrell junction . . . . . . . . . . . . . . . 2.16 The KV = F matrices divided into sub-matrices of similar characters (splittable, non-splittable nodes) to facilitate the calculation of the Schur complement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Schematic representation of a configuration where created node interacts with one of the connections of the other node . . . . . . . . . . . . . . . . 2.18 Schematic representation of a split option for a node of three connections .

Table 1 .

 1 1 -Different types of steels used in nuclear reactors. Details of the crystallographic structure, mechanical properties and corrosion resistance are indicated for each type. Examples of the components for each type of steel is also listed. (after[2][3][4])

		Carbon Steel Low-Alloy Steel	Austenitic S.S.
	Crystal structure	BCC	BCC	FCC
	Microstructure	Ferrite, Bainitic Ferrite, Bainitic	Austenitic
	Main alloying elements Mn, Si	Mn, Mo, Ni, Cr	∼10%Ni, ∼18%Cr
	Yield stress (MPa)	250-450	250-450	∼200
	Toughness/Ductility	High	High	Very high
	DBT	Yes	Yes	No
	Irradiation resistance	N/A	Moderate	High
	Corrosion resistance	Risk for FAC	Moderate	High

Table 1 .

 1 2 -Chemical composition of RPV in French nuclear reactors in weight percent as specified in the French code of design and construction (RCC-M)[8].

	1.2. Generalities on the RPV steel

Table 1 .

 1 3 -Reaction summary between 1/2[1 1 1] screw dislocation and square 1 0 0 radiation-induced loops of 128 SIAs in α-iron at different temperatures. The critical stress represents the stress needed to overcome the obstacle without the flow stress of the dislocation alone (after Terentyev et al. 2010)[START_REF] Terentyev | Reactions between a 1/2<111> screw dislocation and <100> interstitial dislocation loops in α-iron modelled at atomic scale[END_REF].

	Reaction Temperature (K) Burgers	Outcome	CSS (µb/L)
	A1	100	[001]	Restored as [001] loop	0.40
	B1	100	[010]	No reaction	0.09
	C1	100	[100]	Bi-loop	0.56
	A3	300	[001]	Restored as [001] loop	0.42
	B3	300	[010]	Absorbed as helix	0.68
	C3	300	[100]	Absorbed as helix	0.68
	A6	600	[001]	Bi-loop	0.59
	B6	600	[010]	Absorbed as helix	0.62
	C6	600	[100]	Absorbed as helix	0.62

Table 1 .

 1 4 -Critical hardening stress and strength coefficient α for different sizes of SIA loops at T = 300 K (after Liu et Biner 2008)[START_REF] Liu | Molecular dynamics simulations of the interactions between screw dislocations and self-interstitial clusters in body-centered cubic Fe[END_REF].

		∆σ (MPa) α
		300 K 300 K
	SIA 37	320	0.51
	SIA 127	260	0.38
	SIA 271	250	0.38

Table 3 .

 3 1 -Energy and geometrical simulation parameters we explored to calculate the Orowan stress.

	Parameter	Symbol	Range	Unit
	Regularization parameter	a 0	1.5, 2.5, 4.0, 5.0	Å
	Core energy parameter	ζ core	0.09, 0.43, 0.86, 1.28	-
	Poisson ratio	ν	0.0, 0.2, 0.33, 0.435, 0.495	-
	Inter-obstacle distance	L	100.0, 316.2, 1000.0	b
	Obstacle size	D	10.0, 31.62, 100, 316.2	b

Table 3 .

 3 2 -List of the adjusted parameters in DD simulations to reproduce the MD results for iron at 300 K.

	Parameter	Symbol Unit Value	Ref.
	lattice parameter	a	nm	0.286 Terentyev et al. 2008
	shear modulus	µ	GPa	65	Terentyev et al. 2008
	Poisson ratio	ν	-	0.48	Terentyev et al. 2008
	core energy parameter	ζ core	-	0.56	Shukeir et al. 2019
				0.26	Shi 2015
	core width parameter	a 0	nm	0.25	Shukeir et al. 2019
				0.14	Shi 2015

Table 3 .

 3 3 -Identification of mobility law parameters in iron at 300 K. The drag coefficients and the friction stress of different glide systems are indicated (after Shi 2015[50]).

Table 4 .

 4

	2 -Split options of node (A) shown in Figure (4.2b) after the collision between
	the 1/2[1 1 1] screw dislocation with a [0 0 1] loop.
	* VN: number of virtual nodes, ** VJ: number of virtual junctions, dof are 0,1 or 2
	Split option VN* dof VJ**	Burgers	Plane	Dissipation rate (eV /ps)

Table 4 .

 4 3 -Split options at the time step after the collision between the 1/2[1 1 1] screw dislocation with a [0 1 0] loop (case B3). * VN: number of virtual nodes, ** VJ: number of virtual junctions, dof are 0,1 or 2

Split option VN* dof VJ** Burgers Plane Dissipation rate

  (eV /ps)

Table 4 .

 4 5 -Split options tested at node (1) in Figure (4.17) at the time step before the beginning of the formation of the helical turn for the reaction of 1/2[1 1 1] screw dislocation with a 1/2[1 1 1] loop of 127 SIAs. * VN: number of virtual nodes, ** VJ: number of virtual junctions, dof are 0,1 or 2

Split option VN* dof VJ** Burgers Plane Dissipation rate

(eV /ps)

Latin word which means "from the beginning"[17].

MD simulations can also be performed using ab-initio evaluations of effective interatomic forces

Topological operations are not needed in the case of MD simulation or in case of

2D DD[20].

Acknowledgments

First and foremost, I would like to thank my advisors Laurent Dupuy and Benoit Devincre for giving me the opportunity to do research with them. Laurent is not just a fine researcher with great dedication to the field of dislocation dynamics, but he is also passionate and always present when in need. I learned many things from you, especially persistence and dedication. It was a great honor for me to be your PhD student. My supervisor Benoit Devincre, a great mind with a lot of knowledge to share. I would like to thank you deeply for the critical thinking that I learned from you, and how to push each discussion to its limits and far beyond, because this is the essence of research. I am grateful for the knowledge I gained because of you during this PhD.

Acknowledgments

This project has received funding from the Euratom research and training program 2014-2018 under Grant Agreement N • 661913. This work reflects only the author's view and that the Commission is not responsible for any use that may be made of the information it contains. The authors wish to thank Dr. Ghiath Monnet (EDF Lab Renardières) for the fruitful discussions.

General conclusion and perspective

General conclusion

Irradiation has a major influence on the mechanical properties degradation of the RPV due to the existence of radiation-induced defects. There are four types of defects: precipitates, SIA loops, cavities and solute clusters. In this work we concentrate on the interaction between screw dislocations and the SIA loops. The global objective of this PhD work, is to reproduce the individual interactions of screw dislocations in dislocation dynamics in a good agreement with molecular dynamics simulations. Such agreement is characterized by two aspects, the dynamics of the reaction and the critical stress to overcome the obstacles.

A dislocation dynamics code is used, which is a powerful tool to study the interactions of dislocations with radiation-induced defects. The uniqueness of this code is its ability to manage and control collisions and core reactions between dislocation segments. This is done through a set of generic algorithms with the minimum amount of local rules.

In a primary step we identified the parameters in the DD model that allows for better comparison with MD simulations. First, the periodic boundary conditions are implemented. PBC are essential to minimize the calculation time as the size of the simulation cell is reduced. The minimum image convention is used by default. Second, the thermally activated glide mechanism was successfully implemented in NUMODIS. We presented a robust algorithm to detect screw segments at each time step. The implementation was tested with the help of a simulation of a Frank-Read source expansion in pure iron. The temperature dependence is also tested. It is in a good agreement with the results of the analytical solution of the used double-kink model. Finally, in section (3.3), we showed that a large dispersion in DD simulation results can be observed depending on the values of core energy parameter ζ core and the regularization parameter a 0 from the non-singular theory of elasticity.

In chapter (4), a comparison is made of the individual interactions of screw dislocations with SIA loops between DD and MD simulations. The studied SIA loops are square (b = 1 0 0 ) of different orientations and hexagonal (b =[1 1 1]) of different size. The dynamics of each reaction is studied and compared to MD simulations to calibrate our DD code.

Unlike in MD simulations, the studied phenomena should be implemented as predefined models in DD codes. In order to reproduce the studied interactions in this thesis in agreement with MD simulations the following improvements are required:

• an insight analysis of the split node algorithm, whether the problem comes from missing possible split options or simply related to the non-linear solver.

• the convergence problem in the non-linear solver can be divided into two categories;

Appendix: Outcome of Split node algorithm

This is the split options for the first splittable node when the dislocation meets the loop at its corner in case (A3) as reported in section (4.