Numerical Calculations of the Bias Factor of an Edge Dislocation

(1)

Numerical Calculations of the Bias Factor of an Edge

Dislocation

Karl Samuelsson

Reactor Physics Department Royal Institute of Technology

Stockholm, Sweden

(2)

The ability to predict the behavior of the physical and mechanical properties of reactor components during operation is of great impor- tance. It is believed that the edge dislocation’s tendency to preferably attract and annihilate interstitial atoms over vacancies is one of the reasons for the swelling of crystalline materials under irradiation. The bias factor of an edge dislocation can be calculated by solving the diffusion equation with a drift term if the interaction energy between the dislocation and a mobile defect is known; however, an analytical solution only exists for some idealized interaction energies. Due to the development of faster and more powerful computers, it is possible to map the potential energy landscape surrounding an edge dislocation.

Using the results from such calculations, the bias factor would have to be calculated numerically. In this thesis, the possibility of numerically calculating the bias factor of an edge dislocation using the Partial Dif- ferential Equation Tool-box^TMin MATLAB has been investigated. In order to have a reference value, in this method, an interaction energy for which an analytical solution to the diffusion equation exists has been used. This numerical method tends to overestimate the bias factor as the interaction energy is increased.

(3)

Sammanfattning

Förm˚agan att kunna förutsäga hur reaktorkomponenters fysikaliska och mekaniska egenskaper förändras under bestr˚alning är värdefull.

En kantdislokations benägenhet att attrahera och annihilera inter- stitialer framför vakanser tros vara en av anledningarna till varför kristallina material sväller när de utsätts för bestr˚alning. Kantdis- lokationens biasfaktor (p˚a engelska bias factor ) kan beräknas genom att lösa diffusionsekvationen med en driftterm givet en interaktionsenergi mellan kantdislokationen och en punktdefekt. En analytisk lösning existerar emellertid endast för vissa idealiserade interaktionsenergier. Tack vare utvecklingen av kraftfullare datorer är det möjligt att kartlägga interaktionsenergin kring en kantdislokation. Ifall resul- tat fr˚an s˚adana beräkningar skulle användas för att beräkna biasfaktorn för en dislokation, skulle diffusionsekvationen behöva lösas numeriskt. I denna avhandling har möjligheten att beräkna biasfaktorn numeriskt för en kantdislokation genom att använda Partial Differ- ential Equation Tool-box^TM i MATLAB undersökts. För att ha ett referensvärde har en speciell interaktionsenergi betraktats, för vilken det existerar en analytisk lösning till diffusionsekvationen. Metoden som användes visar en tendens att överskatta biasfaktorn för högre interaktionsenergier.

(4)

In this short section I would like to thank Nils Sandberg, who has been my advisor during this work. Nils is also the one who introduced me to the field of radiation material science. I would also like to thank the people at the reactor physics department at KTH for providing me with a helpful and rewarding work environment.

(5)

List of Figures

1.1 A vacancy and an interstitial atom. . . 2

1.2 An edge dislocation is created by inserting a half-plane into an ideal lattice. . . 4

1.3 The Burgers vector. . . 5

1.4 The position of a point defect in polar coordinates with the dislocation in origo. . . 11

1.5 Size interaction energy. . . 12

1.6 Size interaction energy. . . 13

2.1 Flowchart of the procedure of numerically calculating the bias factor. 18 2.2 The transition from 3-dimensional to 2-dimensional geometry of an edge dislocation. . . 21

2.3 Geometry of the solution space. . . 22

2.4 Discretization of a geometry. . . 23

2.5 Geometry of the discretized solution space. . . 24

2.6 The procedure of approximating a solution to a differential equation using the finite element method. . . 28

2.7 Geometry in a square mesh. . . 30

3.1 Results from geometry 1. . . 34

3.6 Comparison between different triangular mesh sizes. . . 39

(8)

3.7 Comparison between different triangular mesh sizes. . . 40 3.8 Comparison between different square mesh sizes. . . 41

(9)

List of Tables

2.1 A list of the different geometries used. . . 20 2.2 Different triangle mesh resolutions. . . 31 2.3 Different square mesh resolutions. . . 32

(10)

(11)

1 Introduction

1.1 General Introduction

Ever since the birth of the first nuclear reactors in the 1940s, physicists have been concerned about the impact on the properties of the materials in the reactor from the high radiation levels that are inherently present during operation [1]. In the case of nuclear reactors, the majority of this radiation is neutrons hailing from atoms that undergo fission in the fuel. When a non-fissile crystalline material is bombarded by neutrons, a collision between an atom and a neutron may, if the neutron has enough kinetic energy, result in a displacement of the target atom, i.e. the atom will be displaced from its original position in the lattice.

The displaced atom may, if enough energy was transferred to it, collide with and displace another lattice atom, and this may continue in several steps, resulting in a displacement cascade [2]. A displacement is equivalent to the creation of an interstitial atom, an atom that is not positioned at a lattice site; and a vacancy, the empty space where the displaced atom was previously positioned. This is illustrated in figure 1.1 (it should be noted that in this work, all figures depicting an atomic lattice shows the case of the simple cubic structure, which is chosen to make the figures easier to understand). Interstitial atoms and vacancies are known as point defects. Any crystalline material at a temperature above 0 K will contain defects due to local fluctuations in the energy [2], however, the point defect concentration is drastically increased during irradiation due to interaction between radiation and lattice atoms. The consequences of the increased number

(12)

vacancy

interstitial

Figure 1.1: A vacancy and an interstitial atom.

of point defects in the reactor material have been the subject of study from both experimental and theoretical scientists for the last 70 years, and continue to puzzle the scientific community in its quest for materials with greater radiation resistance.

1.2 Behavior of Materials Under Irradiation

All crystalline material in a reactor, e.g. the fuel, the cladding, structural material etc. undergoes physical changes during operation. These changes include swelling, segregation, growing and phase change [2]. It is not not difficult to realize that if a structural material in a reactor changes dimensions during operation, it may pose a threat to the integrity and stability of that reactor. In addition to the thermal creep that, due to the high temperatures and pressures present in a reactor, in itself is a problematic behavior, irradiation creep may also occur [3, 4]. This is something that, if not taken into account may further jeopardize the safety of nuclear reactors. These phenomena all change the physical proper-

(13)

1.3 Sinks and Point Defects

ties of the material that is being irradiated, and changes in physical properties might lead to changes in mechanical properties, such as decreased ductility [2].

If the structural material loses ductility it will become brittle, which is most un- wanted during operation and especially in an accident scenario. It is thus is easy to realize that the ability to predict the physical behavior of reactor components during operation is of great value.

If materials for which radiation has little or no impact on the physical and mechanical properties were to be developed, nuclear power could potentially become safer and more economically attractive [5]. It could also help the development of the next generation of nuclear reactors, which are supposed to create less radioac- tive waste. Radiation resistant materials could thus help improve nuclear energy in the aspects where it is most commonly criticized: the safety aspect, the economical aspect, and the environmental aspect. Perhaps it should be noted that the development of the next generation of reactors does not have its motivation from an economical standpoint since it is not expected to decrease the cost of electricity [6], but rather from a safety and environmental related standpoint.

With this in mind, the potential positive impact on nuclear power combined with scientific curiosity is the reason why efforts are made to expand and improve the field of radiation material science.

1.3 Sinks and Point Defects

Point defects in a material have a tendency to group together and form so-called dislocations [7], and these dislocations have a tendency to attract free point defects. Dislocations are also produced when the material is plastically deformed [8]. When a point defect is included into a dislocation, it disappears and this is called annihilation. Anything in a material that, when comes in contact with mobile point defects annihilates them, is called a sink. Sinks can be dislocations, but not necessarily; for example, voids, grain boundaries, and precipitates can also act as sinks.

(14)

1.3.1 The Edge Dislocation

In this work, the sink properties of the so-called edge dislocation are being studied.

An edge dislocation can be seen as the addition of one half-plane of atoms into an atomic lattice, something that will lead to deformation of the atomic planes around it. This could be visualized by taking a book, opening it, cutting out half of a page, and closing the book. This causes the surrounding intact pages to curve around the edge of the cut page. The distortion of the surrounding atoms causes the mobile point defects to drift towards or away from the dislocation. In figure 1.2 an edge dislocation is created by inserting a half-plane into an ideal lattice.

Figure 1.2: An edge dislocation is created by inserting a half-plane into an ideal lattice.

(15)

1.3 Sinks and Point Defects

1.3.2 The Burgers Vector

A concept that is important when defining dislocations is the Burgers vector. It is a measurement of how much and in what direction a lattice is deformed by a dislocation. It is defined as follows [7]:

Imagine that a dislocation line is encircled by taking atom-to-atom steps (see figure 1.3 to the left). Then imagine that the dislocation is removed, but that the same “steps” are made (see figure 1.3 to the right). What previously was an encirclement is no longer a closed circuit (known as a closure failure), and the step that must be taken in order to re-close the circuit is the Burgers vector, b.

b

Figure 1.3: The Burgers vector.

(16)

1.4 The Rate Equations

A simple way to describe the global concentration of interstitials and vacancies in a material is the so-called chemical rate equations:

dC_i

dt = K₀− K_ivC_iC_v − K_isC_iC_s, (1.1) dC_v

dt = K₀− K_ivC_iC_v− K_vsC_vC_s (1.2) where K₀ is the production rate of point defects, K_iv is the interstitial-vacancy recombination coefficient, K_is is the interstitial-sink reaction rate coefficient, and K_vs the vacancy-sink reaction rate coefficient.

Point defects are produced either by interaction between atoms and radiation (if the material is being irradiated), or by emission of point defects from sinks.

This rate of production will depend on the radiation levels and temperature, but also on the concentration and size distribution of sinks.

A point defect is lost when it recombines with another opposing point defect, or when it is absorbed by a sink. A point defect’s susceptibility to annihilation at a sink depends on the type of sink. A sink that has a large attracting force will absorb more point defects than one with a weak force. The force between a sink and a point defect will not only depend on the kind of sink, but also on the kind of point defect. In the case of the edge dislocation, which is the case being studied in this work, the force between the dislocation and an interstitial atom is larger than that between the dislocation and a vacancy. This means that an edge dislocation will attract and absorb more interstitial atoms than vacancies, leaving excess vacancies in the lattice, some of which will group together and form voids.

Another way to describe this would be to say that the interstitial-sink reaction rate coefficient is larger than the vacancy-sink reaction rate coefficient.

This tendency for a dislocation to prefer to absorb interstitial atoms rather than vacancies is believed to be one of the reasons for the swelling of irradiated crystalline materials [9, 10].

(17)

1.5 The Migration of Point Defects

Here, it should be noted that if one tries to predict the swelling of nuclear fuel under irradiation, the rate equations as written in equations 1.1 and 1.2 are not sufficient. The swelling will depend on, among other things, the formation and size distribution of voids in the material. The predictions of these values requires the rate equations to be stated in a more sophisticated manner (see e.g. [11]).

In the rate equations, the loss of point defects due to sinks depends on the sink concentration. One may instead want to study the properties of a certain local sink, that is, regardless of the global sink concentration in the material.

For this reason, a property called the sink strength has been introduced. The sink strength is a measurement of how well a specific sink attracts mobile point defects. The concept of sink strength will be further explained in section 1.5.1, but before that a few other concepts must be introduced.

1.5 The Migration of Point Defects

The migration of mobile point defects in a crystalline material is governed by two processes: diffusion and drift [12]. The diffusion term comes from the random movement of the point defects in the material, and can be expressed by the following equation:

j = −D∇C_α (1.3)

where j is the diffusion flux, D is the diffusion coefficient, and C_α is the concentration of either interstitials or vacancies.

This equation is called Fick’s first law and states that there will be a flow of, in this case, point defects in the direction where the point defect concentration is lower. If the concentration is constant everywhere, ∇C_α will be zero and there will be no net flow of point defects.

The second governing process, known as the drift term, arises from any difference in potential energy throughout the material. A change in the potential

(18)

energy causes the point defect to drift towards the direction where the potential energy is lower. This phenomenon can be compared with a ball rolling down a hill, i.e. the direction where the potential energy is lower. In this work, it is the stress field around a dislocation that causes the change in potential energy, but the change can also be due to applied stress on the material, temperature gradients or chemical potentials. Adding this drift term to the diffusion equation above results in the following equation [9]:

j = −D∇C_α− DC_α∇E k_BT =

β ≡ 1 k_BT

= −D∇C_α− βDC_α∇E (1.4) where k_B is Boltzmann’s constant and T is the absolute temperature. E is the interaction energy between the point defect and the dislocation.

If one assumes steady-state point defect concentrations around the dislocation, the following equation will describe the point defect flux:

dCα

dt = ∇ · j = 0 (1.5)

which, combined with equation 1.4 can be written as:

∇²(DC_α) + β∇(DC_α) · ∇E + βDC_α∇²E = 0. (1.6)

1.5.1 Sink Strength and Bias Factor

If equation 1.4 is solved, and the flux of point defects is obtained, the total point defect current, J_tot into a cylinder around the dislocation can be calculated by:

J_tot = −r₀ Z

dS_r₀br · j. (1.7)

Here S_r₀ is the dislocation surface, br is the normal vector to the dislocation surface, and r₀ is the radius of the circle. This radius can be interpreted as the distance at which a point defect is close enough to the dislocation to be absorbed.

The sink strength of a dislocation is defined by Wolfer [9] as the constant of proportionality between the currents of point defects into a sink and the difference

(19)

1.5 The Migration of Point Defects

in the diffusion potential function far from and at the sink. The diffusion potential function, Φ, is defined as:

diffusion potential = Φ = DC_αe^kBT^E = DC_αe^βE. (1.8) It should perhaps be noted that this potential does not have a meaningful physical interpretation, but is merely a re-written form of DC_α, which simplifies equation 1.4 into:

j = −e^−βE∇Φ (1.9)

and equation 1.5 becomes:

∇²Φ = β(∇E) · ∇Φ. (1.10)

Thus, the sink strength, commonly called k², can be written as:

k² = J_tot

Φ∞− Φ₀. (1.11)

Wolfer continues to define the dislocation bias factor, Z, as the ratio between sink strength with and without a stress field around the dislocation:

Z = k²

k₀² (1.12)

where k₀² is also known as the ideal sink strength. In their paper [12], Wolfer and Ashkin regard Z as “a factor which describes completely the effect of the interaction between point defects and an edge dislocation on the steady-state defect current. Thus we may appropriately call it the bias factor of an edge dislocation.”

The net bias, ∆B, is then defined as the ratio between the interstitial dislocation bias factor and the vacancy dislocation bias factor minus one:

∆B = Zinterstitial

Z_vacancy − 1. (1.13)

Unfortunately, in this field there seems to be a lack in the standardization of the nomenclature. What Wolfer calls net bias is sometimes referred to as the bias

(20)

factor, and what Wolfer calls bias factor is sometimes called the sink capture efficiency [13, 14]. Furthermore, in some cases the sink strength is denoted by a Z [15]. This can easily lead to confusion.

1.6 The Analytical Solution

In the case of an edge dislocation, Z can not be calculated analytically if the surrounding stress field is used in its entirety [12]. However, the majority [9], [16]

of the interaction between the edge dislocation and a point defect comes from the so-called size interaction, which can be approximated as:

E(r, θ) = E_size= −(1 + ν) (1 − ν)

µbV 3π

sin(θ)

r = A ·sin(θ)

r (1.14)

where ν is Poisson’s ratio, µ is the shear modulus, b is the length of one Burgers vector, and V is the relaxation volume. θ and r mark the position of the point defect if the dislocation is in the origin (see figure 1.4). The resulting energy landscape can be seen in figures 1.5 and 1.6. The relaxation volume can be seen as the difference between the point defect volume before it is inserted into the material, and the hole in the material into which the point defect is inserted [10, 17]. This concept is probably easier to imagine if the material is seen as an elastic medium instead of a lattice of atoms. In fact, the expression for the size interaction is derived by replacing the atomic lattice with an isotropic elastic medium with certain elastic moduli [16], which means that the size interaction energy, as given by equation 1.14, is not exact in the real case of an atomic lattice.

If one only uses the size interaction in equation 1.4 , the steady-state solution for j can be found and Z can be calculated.

Frank S. Ham describes in his paper [18] from 1959 a way to analytically calculate the point defect current into an edge dislocation with the assumption

(21)

1.6 The Analytical Solution

edge dislocation line

point defect r

θ y

x

Figure 1.4: The position of a point defect in polar coordinates with the dislocation in origo.

that the size interaction is the only interaction between the point defect and the dislocation. Ham defines the following function:

Ψ(r, θ) = C_αeAβ sin(θ)/2r (1.15) and puts it in equation 1.5 together with equation 1.4 , which results in:

∇²Ψ − Aβ 2r²

2²

Ψ = 0 (1.16)

and has the solution:

e^inθ

α_nI_n Aβ 2r

+ β_nK_n Aβ 2r

. (1.17)

Here, n is an integer, α_nand β_n are constants, and I_n and K_nare modified Bessel functions of the first and second kind.

(22)

−15

−10

−5 0 5 10 15

−10 −15 0 −5

10 5 15

−1.5

−1

−0.5 0 0.5 1 1.5

y [b]

x [b]

β⋅E_size

Figure 1.5: Size interaction energy.

Ham does not use this solution to calculate the so-called bias factor, instead he uses it to calculate an “effective capture radius”, which he defines as the radius a cylinder surrounding the dislocation line would have to have in order to absorb point defects at the same rate, were the potential field around it zero. This result, however, shall not be used in this work. The expression for the effective capture radius derived by Ham is based on the following boundary conditions:

limr→0C_α= 0, (1.18)

limr→0r[∂C_α/∂r + βC_α(∂E/∂r)] 6 0. (1.19)

(23)

1.6 The Analytical Solution

−15

−10

−5 0

5 10

15

−15

−10

−5

0

5

10

15

β⋅Esize

x [b]

y [b]

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

Figure 1.6: Size interaction energy.

Wolfer and Ashkin [12] use a very similar approach, but with a slightly different size interaction;

E(r, θ) = A · cos(θ)

r (1.20)

and a slightly different diffusion potential function, η :

η(r, θ) = DCe^βE/2 (1.21)

together with their own boundary conditions:

η(r₀) = D₀C₀e−Aβ cos θ/2r0, (1.22)

(24)

η(R) = D∞C∞eAβ cos θ/2R (1.23) where R is half the average distance between two dislocations and r0 is the radius of the dislocation.

Using this approach, Wolfer and Ashkin end up with an expression for the dislocation bias factor that they define as in equation 1.12:

Z_edge = ln(R/r₀)

∞

P

n=0

(2 − δ_n0)

Kn(ξr0/R)

In(ξr0/R) − ^K_Iⁿ^(ξ)

n(ξ)

(1.24)

where

ξ =

βA 2r₀

(1.25) and A is a material-dependent constant (see equation 1.14).

1.7 The Numerical Solution

In order to solve the equation

j = −D∇Cα− βDCα∇E, (1.26)

one must first know the gradient of the potential field ∇(E), around the dislocation, which is not trivial. In recent years, due to the development of more powerful computers, calculating the potential field has become possible. This is done using a method called molecular statics. The idea behind this method of calculating the interaction energy is to construct an ideal lattice consisting of a large number of atoms, called a supercell [19], in a computer, where the electronic configuration of every atom is known. Based on the interatomic interactions, all atoms are iteratively moved to their respective “resting location”, i.e. where the forces from the surrounding atoms cancel each other out. The total energy of the system is then calculated. By knowing this energy as well as the total energy of a cell consisting of an edge dislocation and a point defect located at a certain lattice

(25)

1.7 The Numerical Solution

site, it is possible to calculate the interaction energy between the dislocation and the point defect (it is also necessary to know the total energies of the cells with only a dislocation an only a point defect [20]). This procedure would have to be repeated for all possible locations of the point defect in order to map the energy landscape around the dislocation. A molecular statics model does not take into account the thermal vibrations of the atoms, meaning that formally, the model simulates the material at 0 K.

Using the molecular statics method will yield a discrete expression for the energy surface [20, 21]. If one wants to solve the diffusion equation analytically, one must first fit the potential to a continuous surface, and even then, it is not necessarily possible to get an analytical solution, since it might not exist.

This suggests that a more convenient way to obtain the diffusion flux, and ultimately the dislocation bias, is to use the discrete values of the potential energy from molecular statics calculations and solve the diffusion equation numerically.

In this work, the possibility of using this approach is explored more closely by calculating the bias factor of an edge dislocation using the Partial Differential Equation Toolbox in MATLAB. The procedure of the calculation is described in chapter 2.

(26)

(27)

2 Method

As mentioned in the previous chapter, the goal of this work is to examine the possibility of numerically solving the steady-state diffusion equation with the drift term, and using the solution to calculate the bias factor. In order to assess the calculated bias factor using this method, it is compared with the bias factor that is produced using the numerical solution of Ham’s method [18], which means that only the size interaction is taken into account. The bias factor obtained using Ham’s method is taken from [12]. The numerical solver used is the Partial Differ- ential Equation Toolbox in MATLAB, and the calculation process is represented by figure 2.1.

2.1 Stating the Governing Equations

The PDE solver used in the Partial Differential Equation Toolbox numerically finds a u that solves the following differential equation:

−∇ · (c∇u) + au = f on Ω (2.1)

where c,a,f are user-defined functions and Ω is the user-defined geometry. As stated in section 1.5.1, combining equation 1.4 with 1.5 and 1.8 results in the equation:

(28)

Create geometry with boundary conditions and discretize it

State the potential field on the discretized geometry

State the diffusion equation and solve it on the geometry using pdenonlin

State the diffusion equation for the ideal sink and solve it using pdenonlin

Calculate the gradients in the x- and y- directions using pdegrad

Calculate the point defect currents and

sink strengths Calculate the bias factor, Z E(x,y)

Φ(x,y)

Φ

⁰

(x,y)

Φ _x

⁰

(x,y), Φ _y

⁰

(x,y), Φ (x,y), Φ (x,y),

k

²₀

, k

²

x y

Figure 2.1: Flowchart of the procedure of numerically calculating the bias factor.

∇²Φ = β(∇E) · ∇Φ.

If this is to be written on the form of equation 2.1 using the size interaction energy, the parameters c, a and f become:

c = −1, (2.2)

(29)

2.2 Geometry and Boundary Conditions

f = Aβ

u_x

−2xy (x²+ y²)²

+ u_y

x²− y² (x²+ y²)²

, (2.3)

a = 0. (2.4)

In the case of the ideal sink, since there is no interaction between the dislocation and the point defects, the parameters c₀, a₀ and f₀ become:

c₀ = −1, (2.5)

f₀ = 0, (2.6)

a₀ = 0. (2.7)

This differential equation is also known as Laplace’s equation.

2.2 Geometry and Boundary Conditions

The migration of point defects in a crystalline material is essentially a 3-dimensional problem. However, in the case of an edge dislocation that absorbs point defects, the problem can be regarded as 2-dimensional if one chooses to look at the diffusion in a single atomic layer. This is only possible under the assumption that that the edge dislocation line goes straight through the lattice. The transition from 3-dimensional to 2-dimensional geometry is illustrated in figure 2.2.

The geometry on which the differential equation is solved in this work consists of an annulus and can be seen in figure 2.3. The necessity for an inner radius comes from the fact that the analytical expression for the size interaction (see equation 1.14) diverges as limr→0. This inner boundary is also regarded as where point defects are absorbed by the dislocation. According to [9], the diffusion potential function (see equation 1.8) is constant on the inner and outer boundary. In this work the boundary conditions have been chosen to be as follows:

(30)

Φ(r₀) = 0 s⁻¹, (2.8)

Φ(R) = 1 s⁻¹. (2.9)

The Partial Differential Equation Toolbox has a GUI where the user can define the geometry and boundary conditions. In order to better assess the usefulness of the numerical method, a number of geometries with different inner and outer radii are used. A typical geometry as seen in the Partial Differential Toolbox can be found in figure 2.5 A list of all geometries used can be found in table 2.1.

r0 [b] R [b] R/r0

Geometry 1 2 200 100

Geometry 3 5 50 10

Table 2.1: A list of the different geometries used.

In every PDE solver that uses the finite element method, the domain where the solution is calculated is divided, or discretized, into smaller parts. Once this is done, assuming a 2-dimensional domain, the geometry will be defined by a mesh consisting of triangles and nodes (the corners of the triangles), this is illustrated in figure 2.4. The numerical solution of the equation is given at each node, which means that the more triangles the geometry consists of, the higher the resolution.

Higher resolution generally corresponds to longer computation time.

(31)

Figure 2.2: The transition from 3-dimensional to 2-dimensional geometry of an edge dislocation.

(32)

r

₀

R

Ω

∂Ω

Figure 2.3: Geometry of the solution space.

(33)

triangle node Ω

discretization

∂Ω

Figure 2.4: Discretization of a geometry.

(34)

−10 −5 0 5 10

−10

−8

−6

−4

−2 0 2 4 6 8 10

x [b]

y [b]

Figure 2.5: Geometry of the discretized solution space. Note that not all of the geometry is shown here, but only the area closest to the inner boundary.

(35)

2.3 The Solver

The numerical solver in Partial Differential Toolbox uses the so-called finite element method (FEM) to numerically approximate the solution to differential equations [22]. The procedure of solving a non-linear differential equation with the finite element method can be summarized in figure 2.6, and is explained in more detail in this section. The explanation for the procedure can be found in [22].

In order to examine if u does indeed solve the differential equation:

−∇ · (c∇u) + au = f on Ω one defines the function r:

r(u) = −∇ · (c(u)∇u) + a(u)u − f (u) (2.10) which is equal to zero if u is the solution to the differential equation.

Multiplying r with an arbitrary test function v, and integrating it over the domain Ω will, if u is the solution to the differential equation, result in:

Z

Ω

−(∇ · (c(u)∇u))v + a(u)uv dx = Z

Ω

f (u)v dx. (2.11) Using Green’s formula, the equation above can be re-written as:

Z

Ω

(c(u)∇u) · ∇v + a(u)uv dx − Z

∂Ω

n · (c(u)∇u)v ds = Z

Ω

f (u)v dx (2.12)

where n is the outward pointing unit normal vector. The boundary conditions can be defined in the general form:

n · (c(u)∇u) + q(u)u = g(u) on ∂Ω. (2.13)

(36)

Combining equation 2.12 with the boundary conditions defined for the geometry yields the following expression:

Z

Ω

(c(u)∇u) · ∇v + a(u)uv dx − Z

∂Ω

(−q(u)u + g(u))v ds = Z

Ω

f (u)v dx (2.14) and if we want to find a solution to r(u) = 0, a u must be found that satisfies:

Z

Ω

(c(u)∇u) · ∇v + a(u)uv − f (u)v dx − Z

∂Ω

(−q(u)u + g(u))v ds

= 0.

(2.15) Up until this point, all functions are continuous, which is not the desired case if one wants to solve the differential equation numerically. Hence, u and v must be discretized. Let uh be a piecewise linear approximation of u, with as many pieces as the geometry has nodes:

u_h =

Nn

X

j=1

U_jφ_j and v = φ_i (2.16)

where N_n is the number of nodes that defines the discretized geometry, φ_j(x) is 1 at the j:th node and zero everywhere else, and U contains the values for u_h at every node.

Equation 2.15 can now be written in the discrete form:

Nn

X

j=1

Z

Ω

Ujc∇φj· ∇φi+ Uja · φj · φidx − Z

∂Ω

Ujq · φj· φids Z

∂Ω

Ujqφiφj ds

−

− Z

Ω

f φ_idx − Z

∂Ω

gφ_ids = 0 for all indices i. (2.17) If one introduces the following functions:

Ki,j = Z

Ω

c∇φj· ∇φidx, (2.18)

M_i,j = Z

Ω

aφ_j· φ_i dx, (2.19)

(37)

2.3 The Solver

Q_i,j = Z

∂Ω

qφ_j · φ_ids, (2.20)

F_i = Z

Ω

f φ_i dx, (2.21)

Gi = Z

∂Ω

gφids (2.22)

equation 2.17 can be re-written as:

(K + M + Q)U − (F + G) = 0. (2.23)

Note that this is only true if u is the solution to the differential equation. The residual of any solution vector U , regardless of whether or not it solves the differential equation, is defined as

ρ(U ) = (K + M + Q)U − (F + G). (2.24) If the norm of ρ(U ) is zero, then u is the solution to r(u) = 0.

The PDE solver needs an initial guess for U in order to start finding a solution.

If Uⁿ is the first guess, the solver calculates the residual vector for this proposed solution, and if the residual norm is not close enough to zero, a better solution is needed. An improved solution Uⁿ⁺¹ can be found using the following procedure:

∂ρ(Uⁿ)

∂U

(Uⁿ⁺¹− Uⁿ) = −αρ(Uⁿ) (2.25) where α is known as the damping coefficient. If α is chosen to be small enough, then:

ρ(Uⁿ⁺¹)

< ||ρ(Uⁿ)|| (2.26)

and Uⁿ⁺¹ can be calculated as:

Uⁿ⁺¹= Uⁿ+ αp_n (2.27)

(38)

where p_n is called the descent direction and is defined as:

p_n = ρ(Uⁿ) · 1

∂ρ(Uⁿ)

∂U

. (2.28)

The residual vector of Uⁿ⁺¹ can now be calculated, and if necessary an improved solution Uⁿ⁺² can be obtained from Uⁿ⁺¹. This iteration process can now be performed until the residual norm is close enough to zero.

Calculate residual vector ρ(U

ⁿ

)

Geometry and boundary conditions

K,M,Q,F,G

U

ⁿ

ρ(U

ⁿ

) Calculate descent direction p

_n

Choose damping coefficient α

U

^{n + 1}

= U

ⁿ

+ α·p

_n

Guess U

Figure 2.6: The procedure of approximating a solution to a differential equation using the finite element method.

2.4 The Gradients

Once an acceptable numerical solution has been obtained, the gradients of this solution must be calculated in order to calculate the point defect current. The

(39)

2.4 The Gradients

Partial Differential Toolbox has a function (pdegrad) that, given a solution vector, returns the gradients of the solution in the x- and y-directions. Before the point defect current can be calculated the gradient vectors are mapped onto a quadratic mesh (see 2.7). This is done in order to make the flux calculations more convenient: before the mapping, the gradients are given at each triangle center as a vector with the same length as the number of triangles in the mesh. Since there is no simple way to see which triangles are closest to the inner boundary (where we want to calculate the flux), a transformation to a square mesh is performed.

This transformation is done by the MATLAB function tri2grid.

Once the gradients are mapped on the square mesh, the flux in the x- and y-direction at every square can be calculated by:

j_x(x, y) = −u_x(x, y) · e^−E(x,y), (2.29)

j_y(x, y) = −u_y(x, y) · e^−E(x,y) (2.30) and the total current of point defects is obtained by summing the flux at all squares, nn (nearest neighbor), closest to the inner boundary:

J_tot = constant ·X

nn

j_x· ˆr_x+ j_y · ˆr_y. (2.31) Note that this is the point defect current multiplied by a constant that depends on the inner radius and the resolution of the square mesh. However, this constant is unimportant since the same constant shows up when calculating the ideal sink current, causing the constant to cancel itself out when calculating the bias factor.

With the point defect currents being calculated, the sink strengths (remember:

one sink strength with and one without the interaction between the dislocation and point defect) are simply calculated using the equation:

k² = J_tot Φ∞− Φ₀ and the bias factor can finally be calculated by:

(40)

Figure 2.7: Geometry in a square mesh. Notice that not all of the geometry has been mapped onto the square mesh, but only the area closest to the inner boundary.

Z = k² k₀².

2.5 The Resolution

The impact on the calculated bias factor, caused by changes in the resolution has been investigated. Here, we talk about two different kinds of resolutions:

(41)

2.6 Other Input Parameters

the triangle mesh resolution and the square mesh resolution. The triangle mesh resolution refers to the amount of triangles that defines the geometry on which the differential equation is solved. The square mesh resolution refers to the number of squares that the gradients are mapped onto.

All resolution comparisons have been made on the same geometry; the one with an inner radius of 6 b and an outer radius of 600 b. This is Geometry 5 in table 2.1.

2.5.1 The Triangular Mesh Resolution

The different triangle mesh resolutions used can be found in table 2.2. Each mesh refinement splits every triangle into four new ones.

No. of triangles No. of nodes

Resolution 1 810 425

Resolution 5 207360 104000

Resolution 6 829440 415360

Table 2.2: Different triangle mesh resolutions used for the geometry with R = 600 b and r₀ = 6 b (Geometry number 5 in table 2.1). Note that a higher resolution index indicates a higher resolution.

2.5.2 The Square Mesh Resolution

The different square mesh resolutions used can be found in table 2.3.

2.6 Other Input Parameters

In the calculations, different values for Aβ have been used, where A is a measurement of the magnitude of the interaction between the point defect and the

(42)

Mesh size Total no. of squares

Square mesh 1 20×20 400

Square mesh 2 100×100 10000 Square mesh 3 200×200 40000 Square mesh 4 400×400 160000

Table 2.3: Different square mesh resolutions used for the geometry with R = 600 b and r₀ = 6 b (Geometry number 5 in table 2.1). Note that a higher square mesh index indicates a higher resolution.

dislocation (see equation 1.14), and β = _k¹

BT. The different values used for Aβ are: 1, 5, 10, 15, and 20 b (b being the length of one Burgers vector). These values are typical values for some metals found in steels at half of their respective melting temperatures [9].

(43)

3 Results and Discussion

Using the aforementioned method, numerical calculations of the bias factor of an edge dislocation have been performed. All calculations have been performed under the assumption that the size interaction is the only interaction present between the point defects and the dislocation. The results of these calculations have been compared to the expression for the bias factor that Wolfer and Ashkin [12] derived using the method conceived by Ham [18]. The significance of the two different resolutions mentioned in section 2.5 has also been investigated.

3.1 Calculating the Bias Factor

Results from the calculations of the bias factor for different geometries (see section 2.2) and interaction magnitudes (see section 2.6) are presented below. In all cases, square mesh 2 (see table 2.3) has been used together with a triangle mesh resolution approximately the same as resolution 6 in table 2.2 (the number of triangles and nodes are not identical in the different geometries, but the meshes have all been refined the same number of times). In all figures, the analytical value refers to the expression found in equation 1.24.

As can be seen in figures 3.1 - 3.5, the numerical method used produces bias factors that resemble their corresponding analytical values. However, as the interaction magnitude increases, the method starts to overestimate the bias factor. For all geometries, the largest overestimation is when |Aβ| is equal to 20 b.

(44)

The largest of these overestimations is in the case of geometry 1 (see figure 3.1) where the calculated bias factor is 4% larger than its corresponding analytical value. This is the geometry that has the smallest inner radius, and thus the largest interaction energy at the inner boundary.

1 5 10 15 20

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7

Comparison (R = 200b, r

0 = 2b)

|A0β| [b]

Zedge

analytical value numerical calculation

Figure 3.1: Results from geometry 1. Note that the overestimation is larger for larger interaction magnitudes.

(45)

3.1 Calculating the Bias Factor

1 5 10 15 20

1 1.05 1.1 1.15 1.2 1.25

|A0β| [b]

Zedge

0 = 4b)

(46)

1 5 10 15 20 1

1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16 1.18

|A0β| [b]

Zedge

0 = 5b)

(47)

3.1 Calculating the Bias Factor

1 5 10 15 20

1 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16 1.18

|A0β| [b]

Zedge

0 = 5b)

(48)

1 5 10 15 20 1

1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.15

|A0β| [b]

Zedge

0 = 6b)

(49)

3.2 The Triangular Mesh Resolution

The results from the triangular mesh resolution investigation are found in this section. All calculations here have been performed using the square mesh 2 in table 2.3.

0 5 10 15 20

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2

|Aβ| [b]

Zedge

Comparison between different resolutions (R = 600b, r

0 = 6b)

analytical value resolution 1 resolution 2 resolution 3 resolution 4 resolution 5

Figure 3.6: Comparison between different triangular mesh sizes. A higher resolution index corresponds to a higher resolution. The specifics of all triangle resolutions can be found in table 2.2.

As seen in figures 3.6 and 3.6, the calculated bias factors seem to approach the analytical values as the triangle mesh becomes finer. It is not entirely unexpected that a better result is produced by increasing the number of nodes, since in the

(50)

MATLAB function pdegrad, the gradients are evaluated in each triangle center by checking the value of the solution at each node [22]. This means that a finer triangle mesh will produce a better approximation of the gradients.

0 5 10 15 20

1 1.05 1.1 1.15

|Aβ| [b]

Zedge

Comparison between different resolutions (R = 600b, r

0 = 6b)

analytical value resolution 5 resolution 6

Figure 3.7: Comparison between different triangular mesh sizes. A higher resolution index corresponds to a higher resolution. The specifics of all triangle resolutions can be found in table 2.2.

(51)

3.3 The Square Mesh Resolution

In this section, results from the square mesh resolution investigation are presented. Here, triangle mesh resolution 6 (see table 2.2) has been used for all calculations.

0 5 10 15 20

1 1.05 1.1 1.15

|Aβ| [b]

Zedge

Comparison between different square mesh resolutions (R = 600b, r

0 = 6b)

analytical value square mesh 1 square mesh 2 square mesh 3 square mesh 4

Figure 3.8: Comparison between different square mesh sizes. A higher square mesh index corresponds to a higher resolution. The specifics of all square mesh resolutions can be found in table 2.3.

One interesting result from the mesh size investigation is that the bias factor does not approach the analytical value when the square mesh becomes finer (given a fixed triangular mesh resolution), see figure 3.8. This tendency remains unexplained, but could probably be avoided altogether if the square grid mapping

(52)

could be avoided. It should probably be noted that the square mesh resolution that yields the best value (square mesh 2) contains roughly the same amount of squares as the triangle mesh contains nodes. As mentioned in section 2.4, the reason for the square mapping to begin with is the fact that there is no easy way to know which triangle is closest to the boundary in the way that the triangle mesh is defined in MATLAB. It is, however, still possible.

(53)

4 Conclusions

As previously stated, the purpose of this work was to investigate the possibility of calculating the bias factor numerically given an interaction energy, e.g. calculated from molecular statics (see [21] for example). This was done using an interaction energy corresponding to a known analytical expression for the bias factor. The results of this work indicate that, given a high enough refinement of the triangular mesh, a value for the bias factor can be calculated. The value of the bias factor, however, is overestimated for high interaction magnitudes. The overestimation is reduced if the triangle mesh resolution is increased.

Perhaps the largest limitation connected to this method is the fact that the Partial Differential Toolbox is not able to generate a triangular mesh on a geometry if the ratio between the outer and inner radius exceeds a certain value (∼ 100). Part of the motivation for this work was to find a way to calculate the bias factor for an arbitrary 2-dimensional geometry, in addition to interaction energy, and as long as this limitation remains, the method used in this work will be unable to do so.

If this calculation was to be performed for an arbitrary interaction energy, repeated calculations with increasing resolution would have to be carried out until the bias factor stabilizes, i.e. does not change when the resolution is increased.

Some of the code would also have to be re-written in order to use a calculated

(54)

energy landscape as input. Furthermore, it would probably be worthwhile to find a way to avoid the square mesh mapping used in this method.

Finally, it should be mentioned that there are probably other methods one could use in order to calculate the bias factor. The question whether or not these other methods are more useful than the one presented in this work (or if they even exist) should be addressed , especially if results from this code, or one similar to it, were to be used in any attempt to predict the behavior of irradiated materials.

This question, however, is not part of the scope of this work.

(55)

Bibliography

[1] E. P. Wigner, “Theoretical physics in the metallurgical laboratory of Chicago,” Journal of Applied Physics, vol. 17, no. 11, pp. 857–863, 1946.

1

[2] G. S. Was, Fundamentals of Radiation Materials Science: Metals and Alloys.

Springer, 2007. 1, 2, 3

[3] L. Mansur, “Theory and experimental background on dimensional changes in irradiated alloys,” Journal of Nuclear Materials, vol. 216, no. 0, pp. 97 – 123, 1994. 2

[4] K. Ehrlich, “Irradiation creep and interrelation with swelling in austenitic stainless steels,” Journal of Nuclear Materials, vol. 100, no. 13, pp. 149 – 166, 1981. 2

[5] L. Hallstadius, S. Johnson, and E. Lahoda, “Cladding for high performance fuel,” Progress in Nuclear Energy, vol. 57, no. 0, pp. 71 – 76, 2012. Nuclear Materials: Selected articles from the E-MRS 2011 Spring Meeting. 3

[6] E. Schneider and W. C. Sailor, “Nuclear fission,” Science & Global Security, vol. 14, no. 2-3, pp. 183–211, 2006. 3

[7] D. Hull and D. J. Bacon, Introduction to Dislocations, Fourth Edition.

Butterworth-Heinemann, 2001. 3, 5

[8] J. P. Hirth and J. Lothe, Theory of Dislocations. Krieger Pub Co, 1992. 3 [9] W. G. Wolfer, “The dislocation bias,” Journal of Computer-Aided Materials

Design, vol. 14, no. 3, pp. 403–417, 2007. 6, 8, 10, 19, 32

(56)

[10] P. T. Heald and M. V. Speight, “Point defect behaviour in irradiated materials,” Acta Metallurgica, vol. 23, no. 11, pp. 1389–1399, 1975. 6, 10

[11] B. N. Singh, S. I. Golubov, H. Trinkaus, A. Serra, Y. N. Osetsky, and A. V.

Barashev, “Aspects of microstructure evolution under cascade damage conditions,” Journal of Nuclear Materials, vol. 251, pp. 107–122, 1997. 7 [12] W. G. Wolfer and M. Ashkin, “Diffusion of vacancies and interstitials to edge

dislocations,” Journal of Applied Physics, vol. 47, no. 3, pp. 791–800, 1976.

7, 9, 10, 13, 17, 33

[13] I. W. Chen, “Anisotropic diffusion of point defects to edge dislocations,”

Journal of Nuclear Materials, vol. 125, no. 1, pp. 52–63, 1984. 10

[14] S. Golubov, A. Barashev, B. Singh, and R. Stoller, “Dislocation bias revis- ited.” Fusion Materials Semiannual Progress Report for the Period Ending June 30, 2010, DOE-ER-0313/48, Prepared by Oak Ridge National Labora- tory for the U.S. Department of Energy. 10

[15] D. Mordehai and G. Martin, “Enhanced annealing of the dislocation network under irradiation,” Physical Review B - Condensed Matter and Materials Physics, vol. 84, no. 1, 2011. 10

[16] R. Bullough and R. C. Newman, “The kinetics of migration of point defects to dislocations,” Reports on Progress in Physics, vol. 33, no. 1, pp. 101–148, 1970. 10

[17] E. J. Mittemeijer, Fundamentals of Materials Science: The Microstructure- Property Relationship Using Metals as Model Systems (Graduate Texts in Physicsna). Springer, 2011. 10

[18] F. S. Ham, “Stress-assisted precipitation on dislocations,” Journal of Applied Physics, vol. 30, no. 6, pp. 915–926, 1959. 10, 17, 33

[19] C. Domain, C. S. Becquart, and J. Foct, “Ab initio study of foreign interstitial atom (c, n) interactions with intrinsic point defects in α-Fe,” Phys. Rev.

B, vol. 69, p. 144112, Apr 2004. 14

(57)

BIBLIOGRAPHY

[20] E. Clouet, “The vacancy-edge dislocation interaction in fcc metals: A comparison between atomic simulations and elasticity theory,” Acta Materialia, vol. 54, no. 13, pp. 3543–3552, 2006. 15

[21] E. Kuramoto, K. Ohsawa, and T. Tsutsumi, “Computer simulation of defects interacting with a dislocation in Fe and Ni,” Journal of Nuclear Materials, vol. 283-287, no. PART II, pp. 778–783, 2000. 15, 43

[22] “Partial Differential Equation Toolbox^TM User’s Guide.” http://www.

mathworks.com/help/pdf_doc/pde/pde.pdf, June 2012. MathWorks, Inc.

25, 40

Numerical Calculations of the Bias Factor of an Edge Dislocation