Simulation of the Interaction Mean Free Path between Neutrons and TRISO Particles

(1)

Degree Programme in Energy Technology Bachelor’s Thesis

SIMULATION OF THE INTERACTION MEAN FREE PATH BETWEEN NEUTRONS AND TRISO PARTICLES

15^th October, 2013 Lauri Halla-aho

(2)

Lappeenrannan teknillinen yliopisto Teknillinen tiedekunta

Energiatekniikan koulutusohjelma Lauri Halla-aho

Neutronien ja TRISO-partikkelien v¨alisen vuorovaikutusvapaamatkan mallin- nus

Kandidaatinty¨o 2013

39 sivua, 9 kuvaa, 7 taulukkoa ja 5 liitett¨a Tarkastaja: Ville Rintala

Hakusanat: neutroni, vapaa matka, kuulakekoreaktori, polttoainekuula Keywords: neutron, free path, pebble bed reactor, fuel pebble

Neutroneiden ja TRISO-partikkelien välistä vuorovaikutusvapaamatkaa simuloidaan käyt- täen MATLABissa kirjoitettuja ohjelmia reaktorifysiikkakoodi Serpentissä kuulakekoreak- torien simuloinnin tuloksien pakkaustiheyden mukana kasvavan virheen ratkaisemiseksi.

Neutroneiden kulkua seurataan sek¨a rajoittamattomassa ett¨a rajoitetussa avaruudessa.

Neutronien rata lasketaan lineaarisesti ohjelmasta riippuen suoraan neutronien paikkavek- torien ja kaikkien polttoainepartikkelien pintayhtälöiden avulla; jakamalla avaruus ali- avaruuksiin, joissa jokaisessa on murto-osa kokonaispartikkelimäärästä, ja valitsemalla niiden aliavaruuksien partikkelit, joiden läpi neutroni kulkee; tai valitsemalla neutronin kulkusuuntaan muodostetun äärettömän sylinterin sisälle jäävät partikkelit. Vertailuko- hteena käytetään Serpentissä nykyisin käytettävän eksponentiaaliseen jakaumaan perus- tuvan analyyttisen mallin vapaan matkan arvon estimaattia. Serpentin implisiittisen mallin tulokset viittaavat liian suureen vapaan matkan arvoon korkeilla pakkaustiheyk- sillä. Saadut tulokset tukevat tätä havaintoa antamalla 17 % pakkaustiheydellä noin 2.46 % pienemmän vapaan matkan kuin referenssimalli. Tätä tukee neutronin kokema pakkaustiheys, jonka simulaation tuloksena saatiin 17.29 % pakkaustiheys. Tämän lisäksi havaittiin, että polttoainehiukkasten pinnoilta lähtevät neutronit eivät hidasteesta lähtevi- en tavoin noudata täysin eksponentiaalista jakaumaa. Käytetyn mallin jakauma ei siten ole sellaisenaan kelvollinen neutroneiden vapaiden matkojen määrittämiseen.

(3)

Lappeenranta University of Technology Faculty of Technology

Degree Programme in Energy Technology Lauri Halla-aho

Simulation of the Interaction Mean Free Path between Neutrons and TRISO Particles

Bachelor’s Thesis 2013

39 pages, 9 figures, 7 tables and 5 appendices Examiner: Ville Rintala

Keywords: neutron, free path, pebble bed reactor, fuel pebble

The interaction mean free path between neutrons and TRISO particles is simulated using scripts written in MATLAB to solve the increasing error present with an increase in the packing factor in the reactor physics code Serpent. Their movement is tracked both in an unbounded and in a bounded space. Their track is calculated, depending on the program, linearly directly using the position vectors of the neutrons and the surface equations of all the fuel particles; by dividing the space in multiple subspaces, each of which contain a fraction of the total number of particles, and choosing the particles from those subspaces through which the neutron passes through; or by choosing the particles that lie within an infinite cylinder formed on the movement axis of the neutron. The estimate from the current analytical model, based on an exponential distribution, for the mean free path, utilized by Serpent, is used as a reference result. The results from the implicit model in Serpent imply a too long mean free path with high packing factors. The received results support this observation by producing, with a packing factor of 17 %, approximately 2.46

% shorter mean free path compared to the reference model. This is supported by the packing factor experienced by the neutron, the simulation of which resulted in a 17.29 % packing factor. It was also observed that the neutrons leaving from the surfaces of the fuel particles, in contrast to those starting inside the moderator, do not follow the exponential distribution. The current model, as it is, is thus not valid in the determination of the free path lengths of the neutrons.

(4)

LIST OF SYMBOLS AND ABBREVIATIONS

Latin

a, b, c quadratic equation variable [cm]

D diffusion coefficient [cm]

D mean free distance [cm]

D₀ average transmitting length [cm]

i,j,k variable [−]

L average transmitting length [cm]

l free path length [cm]

l length [cm]

¯l mean free path [cm]

˜l median free path [cm]

M o mode free path [cm]

N amount [−]

ˆ

n normal unit vector [cm]

n number density [1/cm³]

f factor [−]

hR²i mean of squares of particle radii [cm²]

r position vector [cm]

r radius [cm]

r⁰ reduced radius [cm]

S neutron source [1/cm³s]

s angular neutron source [1/cm³s]

s distance [cm]

s⁰ virtual spacing [cm]

t time parameter [−]

v velocity [cm/s]

x,y, z coordinate [cm]

(7)

Greek

error [cm, %]

η reduced density [−]

ξ random variable [−]

ρ number density [1/cm³]

Σ macroscopic cross section [1/cm]

σ standard deviation [cm]

σ² variance [cm²]

φ neutron flux [1/cm²s]

Ωˆ direction unit vector [cm]

Subscripts

0 original

a absolute

C chord-length

C cross section

c cell

l leaked

max maximum

min minimum

n neutron

o observed

o optical

p packing

p plane

pb pebble

pt particle

r relative

s scattering

t total

x x-axis

y y-axis

z z-axis

(8)

Abbreviations

ACE a compact ENDF

BC boundary condition

LWR light water reactor

PBMR pebble bed modular reactor

PBR pebble bed reactor

PuC plutonium carbide

PuO₂ plutonium dioxide

PyC pyrolytic carbon

SiC silicon carbide

ThC thorium carbide

ThO₂ thorium dioxide TRISO tristructural isotropic

UC uranium carbide

UO2 uranium dioxide

(9)

1 INTRODUCTION

Approximately 13.5 % of the annual electricity consumed in the world is produced with nuclear power with 430 commercial reactors across the world [15]. The majority of the operating reactors were built in the 1970s and 1980s and are nearing the end of their operational lives. Thus building new reactors to compensate for the future shutdowns of the present reactors is important. An improved fuel efficiency, better safety and more ef- ficient waste disposal are among the main priorities in the design of new types of reactors.

Improved fuel efficiency is achieved by, for example, increasing the operating temperature of the reactor or using the fissile and fertile materials of the nuclear waste. The latter reduces the leftover nuclear waste significantly in addition to improving the efficiency.

One of the main focuses of the fourth generation nuclear reactors has been their safety.

Passive or inherent safety has been achieved theoretically with a number of them. The reactor type discussed in this Bachelor’s thesis is the pebble bed reactor (PBR). Reactors of this type are designed for inherent safety, an example of which is the ability of its fuel to withstand the highest temperatures achievable in a transient caused by a loss-of-coolant accident (LOCA) before the Doppler broadening causes the thermal power output to decrease.

Various reactor physics codes are used to simulate the behaviour of reactors in operating conditions. These simulations determine for example the neutron flux and power distributions, the burnup of the fuel and the criticality of the reactor. In Serpent the fuel is currently simulated with known particle distributions of fuel pebbles. It is an accurate approach to determine the nuclear interactions between randomly moving neutrons and the fuel material but an implicit one, where the particles are not given a specific location but are randomized as the neutron travels through the fuel, would be faster.

The goal of this thesis is to try to validate or reject the analytical model currently used in Serpent for the interaction mean free path between neutrons and TRISO particles. This will be achieved by writing MATLAB programs to generate empirical estimations for the mean free paths in identical circumstances and then compare them with the estimates given by the used model.

(10)

2 PEBBLE BED REACTORS

PBRs are one of the fourth generation reactor types currently researched. Among its strengths are inherent safety, a possibility to build units with a lower thermal power, continuous refuelling, a higher thermal efficiency and an intermediate helium cycle.

The inherent safety of a PBR arises from the ability of the fuel to withstand temperatures of 1,900 K, which is the maximum temperature the reactor could achieve without proper cooling. Building small 100–200 MW modular power plants would aid the use of nuclear power internationally because they would be significantly cheaper to build than the conventional 1,000–1,500 MW power plants. A PBR is refuelled continuously by removing used fuel pebbles from the reactor and checking them for faults and burnup before either reinserting them in the reactor or replacing them. This type of refuelling removes the need for annual refuelling shutdowns. Using thermal efficiencies of 45 % can be achieved with the use of helium as a coolant in comparison to an efficiency of 32–36

% of light water reactors (LWR). The intermediate helium cycle allows the power plant to produce electricity or heat for hydrogen conversion. [4, p. 331]

Though designed for inherent safety, the PBRs also have weaknesses. They include the release of graphite dust within the primary circuit, leak of oxygen into the reactor, the difficulty of detecting helium leaks and the increased power demand of pumping the required amount of helium through the reactor in order to ensure the necessary cooling.

One of the more probable weaknesses is the formation of graphite dust. It is formed when the fuel pebbles grind against each other during operation and their outer graphite layers are eroded, causing a release of dust particles. These dust particles become acti- vated, causing an additional risk of radiation leakage in the case of a coolant leak, and can reduce the effectiveness of the heat exchangers by covering the heat exchange surfaces.

[11, p. 1] As an example of the unlikely weaknesses, an oxygen leak would ignite the carbon used as a moderator and the burning of the fuel would release the gaseous decay products into the atmosphere with the CO₂-emissions.

2.1 Reactor core

An example of a PBR reactor was the South African Pebble Bed Modular Reactor (PBMR). The PBMR core consists of 350,000 fuel pebbles randomly distributed in an annular vessel. Its core is annular with the control rods in the graphite centre. This

(11)

is due to the uncertainty of the lasting integrity of the central column of fuel pebbles.

The shift of the peak power radially further away from the centre allows for a higher power output while remaining within the margins set, for example, by the highest allowed power density and the heat conductivity of the fuel. The whole reactor is surrounded with graphite reflectors with an approximate thickness of 600 mm. The use of graphite brings the risk of the accumulation of Wigner energy. This is avoided by operating the reactor constantly in high temperatures, for example in PBMR between 536 and 900°C for the inlet and outlet temperatures respectively. [3, p. 33–35]

The helium coolant flows through the bed downwards using the cavities between the fuel pebbles as channels. Helium is an ideal coolant because it is chemically neutral, that is it does not react chemically and thus does not corrode the graphite. In addition it is, in practice, invisible to the neutrons, that is its total macroscopic cross section is small enough in order to be disregarded when examining the neutron physical characteristics of the reactor. In the case of an accident where the active cooling of the reactor is not possible, the decay heat can be removed passively with heat conduction, radiation and free convection. This is due to from the capability of the pebbles to withstand the highest attainable temperature without failing, thus enabling the removal of heat passively.[3, p. 33]

2.2 Fuel

The PBR fuel consists of spherical fuel pebbles that include a design-dependent amount of fuel particles, typically 10,000–15,000, with tristructural isotropic (TRISO) coating.

The pebbles are 60 mm in diameter and they are made of a 50 mm graphite matrix con- taining the fuel particles and a 5 mm thick graphite coating. The particles are usually approximately 1 mm in diameter as shown in Figure 2.1. The graphite of the pebbles and particles acts as the primary moderator.

The TRISO particles consist of a central fuel sphere, filled with a varying mixture of uranium dioxide (UO2), plutonium dioxide (PuO2), thorium dioxide (ThO2), uranium carbide (UC), plutonium carbide (PuC) and thorium carbide (ThC) depending on the design. The fuel centre is covered with layers of carbonaceous buffer, pyrolytic carbon (PyC), silicon carbide (SiC) and another layer of PyC. The low-density carbonaceous buffer layer is designed to retain the gaseous fission products and is 95 µm thick. The following inner PyC layer prevents the gases from leaking out of the particle and is 40 µm thick. The purpose of the SiC layer is to contain the fission products with which the PyC layers struggle. The SiC layer is 35 µm thick. The last layer is the outer PyC layer

(12)

d=50mm d=60mm

Pebble TRISO particle

1

2 43 5

Figure 2.1. A PBR fuel pebble and a TRISO particle, with the layers: 1) fuel (250 µm), 2) carbonaceous buffer (95 µm), 3) inner pyrolytic carbon (40µm), 4) silicon

carbide (35µm), and 5) outer pyrolytic carbon (400 µm).

and is equivalent to the inner PyC layer. [3, p. 11]

To be viable for use a fuel pebble cannot contain more than on the order of 0.001 % defective TRISO particles [2, p. 10]. This applies for the fuel pebbles throughout their use and thus tests have been performed to find different ways in which fuel particle failure can occur. Identified failure mechanisms are pressure induced failure, kernel-coating interaction, SiC dissociation and increase of the porosity of SiC, SiC–fission product interaction and irradiation effects on coating integrity. [3, p. 11–13]

The pressure inside a fuel particle increases during operation as fission gases accumu- late leading to an increasing tension in the SiC layer. When the stress is greater than what the layer can endure the failure of the SiC layer leads to a failure across all the coating layers. The displacement of carbon due to a thermal gradient inside the fuel particle can lead to a contact between the SiC layer and the fuel kernel and thus coat degradation.

Damage to SiC through dissociation or an increase in its porosity under high operating temperatures leads in an increased fission product leakage. High thermal gradients can lead to chemical reactions between fission products, especially palladium, and the SiC layer eventually leading to a coating failure. Lastly bonding between the outer PyC layer and the graphite of the graphite matrix can lead to portions of the outer PyC layer being torn off due to stresses induced by irradiation. [3, p. 11–13]

(13)

3 REACTOR PHYSICS CALCULATION

Nuclear reactor analysis and reactor physics calculation are mostly centred on solving the distribution of neutrons, or neutron fluxes, in a reactor. This is because the neutron flux determines the rates of nuclear reactions present in the reactor. For the determination of the distribution of neutrons their motion through the reactor between their origin and their eventual absorption in or escape out of the reactor must be tracked. This motion is called the neutron transport. [1, p. 103]

The neutron transport can be treated as a diffusion of neutrons similar to the diffusion of gas molecules or calculated using the neutron transport equation. With diffusion the neutrons moving through the reactor are thought to constitute a ”neutron gas” in which neutrons diffuse from highly concentrated areas to the areas with a lower neutron density. This, however, loses accuracy because instead of frequent collisions, as with gases, neutrons tend to travel relatively long distances between collisions due to the small size of nuclei. [1, p. 103]. The diffusion model does have its uses in reactor physics calculations, for example with the energy-dependent diffusion equation [1, p. 141]:

1 v

∂φ

∂t − ∇ ·D(r, E)∇φ+ Σt(r, E)φ(r, E, t)

= Z ∞

0

dE⁰Σ_s(E⁰ →E)φ(r, E⁰, t) +S(r, E, t),

(3.1)

wherev is the velocity of a neutron, [cm/s];φ the neutron flux, [1/cm²s];D the diffusion coefficient, [cm]; Σ_t and Σ_s the macroscopic total and scattering cross sections respectively, [1/cm]; and S the neutron source term, [1/cm³s]. This equation is often used in the derivation of multigroup diffusion equations. More precise results are received with the use of the neutron transport equation [1, p. 113]:

∂n

∂t +vΩˆ ·n+vΣ_tn(r, E,Ω, t)ˆ

= Z

4π

dΩˆ⁰ Z ∞

0

dE⁰v⁰Σ_s(E⁰ →E,Ωˆ⁰ →Ω)n(r, Eˆ ⁰,Ωˆ⁰, t) +s(r, E,Ω, t),ˆ

(3.2)

where n is the number density, [1/cm³]; and s the direction-dependent neutron source term, [1/cm³s]. This equation is the equivalent to the Boltzmann equation associated with gases. Due to the complexity of the equation it is usually solved through various simplifications followed by a deterministic calculation or using Monte Carlo methods to simulate the neutrons and deduce the solution from the statistical results the method produces.

(14)

3.1 Monte Carlo

The Monte Carlo method was first used by Georges Louis Leclerc, Count of Buffon, in 1777 when he used it to find the value of π by dropping toothpicks on a sheet of paper filled with parallel equidistant lines. The actual Monte Carlo method originates in the Second World War during the years 1940–1945 when computers were used to help with the projects of United States Department of Defense aiming to construct atomic bombs.

[5, p. 1]

The Monte Carlo method produces an expected value for a wanted quantity. This is achieved by averaging over a number of computed sample values each dependent on a randomly distributed variable. [5, p. 3] For an example consider the mean free pathl of a neutron inside a fuel pebble. Its value is received from the mean of a number of independent free path lengths λ_i (i = 1,2, . . . , n) which depend on the packing factor f_p of the fuel particles, the radius of the fuel particlesr_pt, the position and direction vectors of the neutrons r_n and Ωˆ_n and the position vectors of the fuel particles r_pt,j (j = 1,2, . . . , k).

Each of these variables can be assumed to be a function of a uniformly distributed random variableξ (0≤ξ≤1). This yields, as the expression of the free path length of a neutron, when all k fuel particles are considered:

λ_h = min

j=1,...,k

f

f_p(ξ), r_pt(ξ),r_n(ξ),Ωˆ_n(ξ),r_pt,j(ξ)

(3.3)

In this case the variables f_p and r_pt are constants for all values of ξ. This calculation is then repeated n times resulting in a sequence of independent variables that, with n sufficiently large, resemble variables from a normal distribution in accordance with the central limit theorem. This distribution then has an expected value as follows:

l=E(λ) = 1

n(λ₁+· · ·+λ_n) (3.4) The error of the resulting expected value converges at a rate of 1/√

n. Due to the slow convergence one of the downsides of Monte Carlo method is its demand for resources and time during calculation.

Monte Carlo reactor physics codes follow the histories of individual neutrons inside a reactor from their ”birth” in a fission to their ”death” in an absorption or leak. The estimation of their reaction rates during the histories can be done analogously or implicitly. Analogous approach treats each simulation of a neutron similarly to a real neutron

(15)

and counts the interactions it undergoes and these interactions are then collected from a plethora of simulations to form statistical estimates for the reaction rates. The implicit approach gives each neutron a statistical weight that is reduced after each interaction instead of terminating its history, thus enabling the neutrons to reach less accessible regions. [6, p. 127–128]. Two methods used to track the movement of the neutrons and interactions are the ray tracing and delta-tracking methods.

3.1.1 Ray tracing

The ray tracing method calculates the interactions of a neutron given its location, and thus the region where it is located, and direction. The optical distance to the next collision is sampled using the macroscopic total cross section of the region. If the distance is greater than the optical distance to the boundary of the nearest material region, the neutron is moved to the given boundary and a new distance is sampled using the properties of the next region. [6, p. 98–99]

If the destination of the neutron is within the current material region the reaction is sampled after the neutron has been moved. In cases of leak and absorption the history is terminated and if scattering occurs the energy of the neutron is reduced by a sampled amount and its new direction is sampled.

3.1.2 Delta-tracking

The delta-tracking method uses a majorant macroscopic cross section for all material regions instead of individual cross sections with the addition of a virtual collision cross section to the macroscopic cross sections of different materials. This means that the interaction probability in all the material regions is the same, thus removing the need to detect surface crossings. When a collision is sampled the ratio of the actual macroscopic cross section of the current region and the majorant is used as the probability of a real collision. This is followed by a sampling to select the reaction that occurred. [6, p. 101]

Delta-tracking thus simplifies the calculation of the optical distance to the nearest surface to a simple check to see if the neutron lies within a surface with the sign of the surface function. Using the majorant to sample the distances between collisions allows modelling any area with a value for the macroscopic total cross section. The presence of absorber materials cause a significant increase in time wasted with computing virtual collisions.

This can be avoided by using the ray tracing method in the vicinity of such regions or using two majorants, one of which excludes any absorber materials, in the calculation [6, p. 102, p. 151].

(16)

3.2 Serpent

Serpent is a program developed for continuous-energy reactor physics burnup calculation and is under development in VTT since 2004. It utilizes the Monte Carlo method and the extended delta-tracking method which is similar to the aforementioned version using two majorants. One of its main properties is the ability to generate group constants for deterministic simulator codes. Serpent uses the ACE-format (A Compact Evaluated Nu- clear Data File) nuclear data libraries to construct its own cross section data on an evenly spaced universal grid to speed up the calculation, due to the lack of a need to iteratively calculate the energy-dependent cross sections on the default irregular energy grid, while requiring an increased amount of memory. [12]

The main purpose of Serpent is the modelling of LWR fuel assemblies but it is suit- able for almost every configuration, including PBR simulation. Currently Serpent utilizes two different models, explicit and implicit, to determine the locations of neutron–particle interactions. If the collisions are solved explicitly, the exact locations of particles within a pebble are known. This produces accurate results but is slow and requires a lot of com- putational memory. Thus a more ideal approach would be an implicit approach where the collisions are solved on the fly. This leads to a significant increase in simulation speed and causes every fuel pebble to have a different fuel particle distribution.

There is, however, an obscure source of error, the origin of which is currently unknown, with one or more of the models currently used with the impicit model in Serpent. The error increases with the increase of packing factor of the fuel particles inside the fuel pebbles. An expected source of this error is the analytical model used for the approximation of the mean free path between neutron–particle interactions. To study this, both the mean free paths of neutrons and the apparent packing factor observed by neutrons have to be simulated.

3.3 Prior studies

Prior studies regarding the mean free paths and mean chord lengths include those of Lu and Torquato [9], Torquato et al. [14] and Liang et al. [8]. Lu and Torquato [9, p. 6479–6481] derived expressions for the exclusion probabilities of particles and mean chord-lengths between the surfaces of those particles. An exact representation was available only for the one-dimensional chord-length distribution and an approximation was used for the calculation of the three-dimensional chord-length distribution and thus its mean chord-length. The resulting mean chord-length was:

(17)

l_C = 1 πρ

1−η

hR²i, (3.5)

where l_C is the mean chord-length, [cm]; ρ the particle number density, [1/cm³]; η the reduced density andhR²i the ensemble average of the squares of particle radii, [cm²]. An application of the mean chord-length they derived is currently used in Serpent.

Torquato et al. [14] studied the nearest-neighbour distribution of spheres from an arbitrary point in space. An exact solution for the three-dimensional case was again unavailable and three approximations were used. The resulting nearest-neighbour distribution cannot, however, be used to represent the mean chord-length because the nearest-neighbour approach does not take into account the direction in which the nearest neighbour lies.

Liang et al. [8] studied the effects of packing factor correction on the criticality calculations using the mean chord-length from the equation (3.5). It was found that the most accurate results were received from the use of a track length estimator to simulate the packing factor experienced by a neutron passing through the region [8, p. 144].

(18)

4 SOLUTION APPROACH

One possible major source for error in the simulations is, as previously mentioned, the interaction mean free path between neutrons and TRISO particles. Murata et al. [10, p. 107] have derived an equation for the packing factor as a function of average transmitting lengths through a sphere and between spheres:

fp = L

L+D₀, (4.1)

wheref_p is the packing factor;L the average transmitting length through a sphere, [cm];

and D₀ the average transmitting length between spheres, [cm]. In Serpent Lepp¨anen [7, p. 3] uses an equation for the mean particle distance derived from (4.1):

D= 4r_pt 3

1 f_p −1

, (4.2)

whereD is the mean distance between fuel particles, [cm]; andr_pt the radius of a fuel particle, [cm]. HereL = 4r_pt/3 is the average transmitting length of a neutron through the fuel particle. It is more commonly known as the mean chord-length [13, pp. 1607–1608].

The mean particle distance is then used to sample the actual distances the neutrons travel with:

s=−Dlog(1−ξ) =−4r_pt 3

1 fp

−1

ln(1−ξ), (4.3)

where s is the optical distance to the nearest fuel particle, [cm]; and ξ a uniformly distributed random variable. To check the validity of this model and in the case of differing results to find the causes of found differences a MATLAB program is written for the simulation of a neutron traversing through a space filled with fuel particles.

The interaction mean free path simulation is performed by using a generated fuel particle distribution through which a neutron is moved. The neutron is randomly located within the fuel space. The goal of the simulation is then to determine the mean distances before collisions with TRISO particles a neutron undergoes while travelling in a constant direction. Determining the particles in the trajectory of the neutron as fast and as effectively as possible is the most important phase of the simulation. Possible methods include using all existing fuel particles simultaneously, dividing the space into subregions and using only those particles the subregion of which has been penetrated, and using all the particles

(19)

within the range of the radius of a fuel particle from the trajectory of the neutron.

4.1 Explicit distribution

In the explicit particle distribution the particles have predetermined positions inside the fuel pebble. They are randomly distributed with a uniform distribution. This makes the calculations simpler because the number of distributions can be, depending on the simulation, small. The use of a small number of distributions does, however, have a detrimental effect on the accuracy of the results as it loses the statistical randomness present with the use of infinite distributions, that is implicitly distributed particles.

Due to the predetermined distributions of particles in a space the simplest way of increasing the accuracy of a single distribution is, given a constant packing factor, by increasing the number of fuel particles and thus the total space. The larger the space the more uniformly distributed the particles will seem as the apparent local concentrations even out. The downside of increasing the number of particles is that it also increases the computation time significantly.

4.2 Linear collision detection

With linear collision detection, after the randomization of the positions of particles as well as the positions for the test neutrons and their directions, the collisions between neutrons and particles are determined. This is performed by solving the equation received from substituting the time-dependent coordinates of the neutrons for the variables in the equations of the fuel particles for a sphere. The coordinates of a neutron are:

r_n=r_n,0+tΩˆ_n =





 x_n,0 y_n,0 zn,0





+t





 Ωˆ_n,x Ωˆ_n,y Ωˆn,z





, (4.4)

wherer_n and r_n,0 are the current and initial position vectors of the neutron respectively, [cm];t a time parameter; andΩˆ_n the direction unit vector of the neutron, [cm]. Similarly the equation for the surface of a fuel particle is:

(x−x_pt)²+ (y−y_pt)² + (z−z_pt)² =r_pt² (4.5) (r−r_pt)·(r−r_pt) = r²_pt, (4.6)

wherexpt,ypt andzpt are the x-, y- and z-coordinates of a fuel particle respectively, [cm];

(20)

andr_ptthe position vector of a fuel particle, [cm]. The x-, y- and z-coordinates from (4.4) are then substituted for the respective variables in (4.6) to form a quadratic equation of the following form.

t²+ 2Ωˆ_n·(r_n,0−r_pt)t+|r_n,0−r_pt|²−r_pt² = 0 (4.7)

Equation (4.7) can be solved using the quadratic equation. The expressions for the variables a,b and c of the quadratic equation are:





 a= 1

b = 2Ωˆ_n·(r_n,0−r_pt).

c=|rn,0−rpt|²−r²_pt

(4.8)

Using these variables t can be solved. The resulting free path length is then:

l =

Ωˆ_n

t =t = −b±√

b²−4c

2 , (4.9)

where l is the free path length, [cm]. A collision between a neutron and a fuel particle then occurs only if the resulting free path length for a given neutron–particle combination is a real number, that isl /∈C.

4.3 Segmented collision detection

In segmented collision detection the fuel pebble is divided into cell segments with evenly spaced parallel planes in the x-, y- and z-directions with N_x, N_y and N_z planes, respectively, into (N_x−1)(N_y−1)(N_z−1) cells. After the division each fuel particle is then given an index value corresponding to the cell it is located in. Next the coordinates and direction vectors of the neutrons are randomized and they are indexed. After indexing the neutrons the points of intersection with the planes are calculated. This can be done using a modification of modular arithmetics. For example consider the x-axis. A neutron with coordinates (x_n, y_n, z_n) moving in the direction ( ˆΩ_x,Ωˆ_y,Ωˆ_z) will cross a yz-plane when:

x_n+ ˆΩ_xt≡0 (mod s_c), (4.10)

where s_c is the cell spacing, that is the distance between successive planes, [cm]. This equation gives, given the x-component of the directional vector is not zero and the pebble is centred at the origin:

(21)

t=− 1

Ωˆ_x(x−ns_c), (n =−N_x−1

2 , . . . ,N_x−1

2 ) (4.11)

This equation can be applied to the y- and z-directions by selecting the corresponding components from the position and direction vectors. After solving the time parameters for all the dimensions the coordinates of the points of intersection can be calculated by substituting the time parameter values into equation (4.4). Combining and sorting all the time parameter vectors the cells through which the neutrons moves can be determined using the original cell and time parameter vectors to judge which planes are crossed. This is done by comparing the time parameter from the combined vector to each individual time parameter vector to find a match. If a match is found in multiple vectors the point of intersection is either on an edge, if the number of matches is 2, or in a corner, if the number of matches is 3.

j+1

j

i i+1

(i,j+1)

(i,j)

(i+1,j+1)

(i+1,j)

Figure 4.1. A 2-D example of a neutron moving into a corner between cells.

In each of these cases the following procedure is to either increase or decrease the indices of the dimensions where a cell border is crossed depending on the signs of the corresponding direction vector components. In the example of figure 4.1 the neutron moves from the cell (i, j), that is the ith cell plane on the x-axis and thejth cell plane on y-axis, into the cell (i+ 1, j+ 1) with the signs of the direction vector components being positive. The indices of the new cell are then added to the index matrix of incident cells.

When all points of intersection have been considered and the cell indices have been gath-

(22)

ered, the fuel particles located in these cells are used simultaneously for the collision detection calculation instead of all the fuel particles inside the pebble. To ensure the inclusion of the particles in the chosen cells that are partially inside the cell even if their centres are in the adjoining cells, the particle indexing must be done using a virtual cell spacing of s⁰ = s+ 2r_pt where each side is extended by a length of r_pt. After choosing the particles the collision detection is performed, for example by using the linear collision detection introduced earlier. Using this method the collisions in both directions are solved by default.

The total number of cells must be chosen small enough for there to be a noticeable benefit for the computation speed but large enough to prevent the border cross check from taking too much time. If a particle radius of 0.046 cm is used the maximum number of cells, with enough room for one fuel particle, inside a graphite matrix with a radius of 2.5 cm is approximately 54 per each dimension. The smaller the spacing the larger the relative effect the addition of space for the virtual cell becomes, and thus the use of cells with the width of one diameter is essentially pointless as its volume is at most octuple its original volume. An illustration of the centre cells can be seen in figure 4.2.

Figure 4.2. An illustration of the cells -2, . . . ,2 on the x-, y- and z-axes.

(23)

4.4 Projectile collision detection

A combination of previous methods can be used to simplify the calculation procedure when compared to the segmented collision detection. The neutron can be thought to travel through an infinite cylinder, the axis of which lies on its direction vector. This cylinder can then be thought to have a radius equal to the radius of a fuel particle, r_pt. The fuel particles the neutron will hit are then those whose centres lie within the cylinder.

This is the case when the following equation holds true:

rpt−rn−

(rpt−rn)·Ωˆn

Ωˆn

2

< r_pt² . (4.12)

The equation (4.12) is received by taking the projection of the distance vector r_pt −r_n to the direction vector and using the received value to find the location of the neutron where it is closest to the fuel particle. The square of the distance is then calculated and compared to the square of the radius of a fuel particle. The particles found to lie within the cylinder are thus ones that are in the direction of the movement of the neutron and the free path length is then the distance to the surface of the nearest particle in the direction of movement.

Tracking the collisions of a neutron entering a region, filled with a given packing factor, with this method can be used to determine the packing factor experienced by the neutron. This can be performed by calculating the distances between every collision the neutron undergoes while passing though the region. The packing factor can then be estimated using the equation (4.1).

4.5 Projection cross section

This method is the simplest of those listed here. It uses randomly generated particle distributions to produce projections on a plane after which the area covered by the projections is determined. Using different numbers of fuel particles the effect of the differing number of fuel particles on the effective cross section, received by dividing the calculated area by the number of particles, is determined.

In the example of figure 4.3 the inside of the fuel pebble is coloured yellow with its border and the fuel particles coloured black. The plotted figures of the projections of particles are then saved as image files. A program able to distinguish between differently coloured pixels is then used to count the black pixels in a figure where only the projections are plotted and the pebble and axes are hidden. Using the average amount of pixels a single

(24)

Figure 4.3. A sample projection of particles on the xy-plane.

projection occupies as the cross section of a single particle, the cross sections of varying numbers of particles can be compared to it to track the decrease in the apparent cross section of the fuel particles as the amount of fuel particles increases.

This method can then be used to derive a correlation between the number of particles, radius of particles and the cross section. This is done by forming a function for the coefficient. Because this method uses the whole pebble to produce the average cross section it brings an error due to the significantly more limited width of the path of a neutron.

4.6 Chosen methods

The first chosen approach for the modelling is the linear collision detection method. It is chosen because it is simple to program due to the need of only solving a group of quadratic equations to find collisions. It can also perform the calculations fast due to the simplicity of the required programming. Some computation time can be wasted, however, on the detection of a collision with fuel particles that are not in the immediate vicinity of the path of a neutron.

The second approach is the segmented collision detection. The aim is to find a way

(25)

for the program to find the indices of the cells the particles are in and to find the cells that the neutron goes through without the cost of a significant increase in computation time. If succeeded this can reduce the time spent on solving the quadratic equation, where the calculation of the square root of the determinant takes the most time, by only checking the few particles in the vicinity.

Lastly the projectile collision detection is used. Due to its simplicity the possible sources of error in programming should be minimized when compared to the previous two. This is due to the speed at which the particles in the path of the neutron can be found, giving the ability to reduce the number of equations to be solved when determining the closest particles. The lack of a need to divide the space into cells also speeds the procedure without adding complexity.

4.7 Boundary conditions

If a neutron does not collide with a fuel particle during its travel inside the fuel pebble it will eventually arrive at the edge of the graphite matrix. After this point there are no fuel particles to be hit in the trajectory and the continued movement of the neutron must be taken into account with the use of one or more boundary conditions (BC).

In an ideal situation the generated space would near infinity in size and thus this could be ignored. Due to the limitations set by computer memory, however, a boundary must be set at a given radius to reduce the amount of capacity needed to perform the calculations.

4.7.1 Collision at boundary radius

The easiest and simplest way of resolving the movement of a neutron past the boundary of the graphite matrix is to terminate the travel of the neutron. Then the free path length of the neutron consists only of either the shortest optical distance to a fuel particle or the boundary radius.

This causes problems as the true lengths of the neutrons are not known, causing the resulting mean free path to appear shorter than it would be if the neutrons were to go on until a collision is encountered. If the fraction of neutrons arriving at the boundary radius is small, on the order of 1 pcm, the effect on the average is relatively small.

4.7.2 Regular or random periodic boundary

Another option for the BC is to have a periodic boundary, be it regular or random by nature. This means that when the neutron arrives at the boundary of the graphite matrix

(26)

it is relocated to another point on the surface of the fuel pebble. In the case of a regular periodic boundary the neutron is relocated on the opposite side of the pebble with the same direction. This can be performed only if a collision has been originally found in the opposite direction. Otherwise the neutron would end up in an infinite loop. The random periodic boundary samples a new location along the boundary and a new direction. Both of these options represent the travel of neutrons through multiple fuel pebbles where they exit one and enter another either with an identical orientation or a random rotation.

This approach does not take any scattering into account due to the graphite because the directions of neutrons stay constant. This is the case even with the resampling of the direction as a consequence of the use of random periodic boundary because the new direction of the neutron in the pebble can be thought to be relative to a static particle distribution common to all the pebbles, that is the sampled location and direction of a neutron are the point of entry and direction of the neutron from the point of view of the pebble.

4.7.3 Reflective boundary

The third option is to reflect the the neutron as it arrives at the boundary surface. The reflected neutron retains its position but its directional vector goes through a reflection against the normal of the local tangential plane. The equation for the reflected vector is:

Ωˆ⁰_n =Ωˆ_n−2

Ωˆ_n·nˆ_(x,y,z) ˆ

n_(x,y,z), (4.13)

whereΩˆ⁰_n is the reflected directional unit vector, [cm]; andnˆ_(x,y,z) the normal unit vector of the tangential plane at the coordinates ( x, y, z), [cm]. The tangential unit vector is received from:

ˆ

n_(x,y,z) =− r_n

|r_n|. (4.14)

If the normal is parallel to the original direction vector and there are no collisions in the opposite direction, the reflective BC cannot be used. In this case a better choice would be to use the random periodic boundary. The use of the reflective BC does reduce the randomness of the movements of the neutrons with the loss of the need to randomize their positions and directions, as opposed to the random periodic boundary. Thus it will not be used in the simulation.

(27)

5 PROGRAMMING

Four MATLAB scripts are written for the calculation of the free path lengths. The first one utilizes a simple linear collision detection using all the fuel particles to find the ones collided with. The second one first finds the particles near the path of a neutron before determining the collisions. The nature of these calculations and their differences means that while the first script can calculate all the collisions of the neutrons simultaneously the second one, as it is, is only able to calculate one neutron at a time. This slows the calculation significantly but also decreases the required amount of memory which can be many gigabytes with the first one.

The third script uses the aforementioned projectile collision detection method, utilizes an extended space and terminates the neutrons at the boundary. Finally the fourth script is used to determine the packing factors experienced by the neutrons passing through the fuel. It also utilizes the projectile collision method to determine the mean free paths of neutrons entering the space from the boundaries.

In this chapter the most important sections of the programming are explained. They include the distribution generation for the fuel particles and neutrons, determining and indexing the cells in which the fuel particles and neutrons reside, collision detection between the neutrons and fuel particles and the free path calculation.

5.1 Distribution generation

Two different methods for particle and neutron distribution methods are used. In the first method the particles are distributed randomly using therandfunction inside a cube with an edge length of 2r_pb after which their locations are checked to lie within the fuel pebble. A simpler functionsph2cartwas tried with the angles and radii randomized but it was found that the received distribution is not uniform but focused near the z-axis.

The distribution generation begins with an anchor particle generated within a reduced pebble radiusr⁰_pb =rpb -rpt. After this the following particles are given locations one by one, every time checking that they do not intersect with the preexisting fuel particles or lie partly outside of the boundary radius.

The aim of this function is to generate the coordinates as fast as possible to reduce the added runtime to the otherwise already time-consuming script. Thus the use of for

(28)

loops has been strictly avoided. Parallelization would speed the generation but due to the lack of parallel versions of the used commands the only ways of making it faster are for example reducing the number of redundant variables and trying to reduce the number of necessary iterations. Generating the variables directly in the spherical coordinates would remove the need to check if the particles are within the pebble. The coordinate transfor- mation to Cartesian can, depending on the method, take more time than that which is spared.

A similar procedure is applied to the generation of the test neutrons, but with a few differences. The maximum distance a neutron can be from the origin is the full radius of a fuel pebble instead of the reduced radius. The initial locations of the neutrons are generated outside of the particles and thus there is an additional check to reject the coordinates found to lie inside fuel particles.

A second method would be to generate a cube with a volume equal to that of a fuel pebble. A given number of particles, precalculated for a chosen packing factor, are then generated inside the cube. The program is written so that it produces ann×n×n space, where n is the number of cubes per dimension. The particles near the edges of a cube are listed in a cell matrix to avoid overlapping at the boundaries between cubes. The neutron generation is then done by choosing the centre cube and randomly picking the fuel particles in which the neutron is generated. The neutron is then given a random direction vector and moved in that direction to the surface of the fuel particle.

5.2 Cell indexing

The space is divided with 27 planes into 25 cells along each axis, creating a total of 15,625 cells. The cell spacing is calculated using the radius of a fuel pebble and the aforementioned axis-specific cell count,n_c:

sc= 2rpb

n_c , (5.1)

wheres_cis the cell spacing, [cm]. The cell numbers are received by dividing the coordinate values by the cell spacing and rounding the resulting value to the nearest integer. Thus the cell indices range from -13 to +13 with the cell centre coordinates being multiples of cell spacing. Using the equations (4.10) and (4.11) for all the axes the plane intersection coordinates can be solved. The time parameters are collected in one vector after which a frequency table is formed to find any duplicate time parameters. The presence of duplicate or triplicate parameters means a cell crossing has occurred either at an edge in the

(29)

case of a duplicate or in a corner in the case of a triplicate.

Following the check for duplicates and triplicates the neutron is tracked through the cell crossings by altering the current cell indices of the neutron. This is done by dividing the combined parameter vector into forward, which has the positive parameters, and backward, which has the negative parameters, vectors. Using these vectors, with their respective intersection coordinates, the indices of the new cell are determined with the coordinates of the crossing point, which are used to determined between a crossing at the face, edge or in the corner of the cubic cell, and the signs of the direction vector. In the case of backward movement the signs are reversed.

This method by itself calculates the cells for only one neutron at a time. This forces the script to loop at least once for every neutron with additional loops for every round without a collision. This causes a significant increase in the total runtime. This problem could possibly be avoided through the use of acellarray where the cell index matrices of neutrons are converted into text strings and then stored in thecell array. This, though requiring each neutron to be calculated separately in aforloop, could reduce the increase in runtime.

5.3 Collision detection

For the collision detection a number of test neutrons must be generated with the aforementioned distribution generation method. If cell indexing is not used, or thecellarray version is used, the amount of test neutrons is restricted by the available random access memory (RAM). For an example in the case of no cell indexing of any kind, using 15,000 fuel particles and 100,000 neutrons, the collision detection calculations of the first round require two matrices that take approximately 11.18 GB of RAM each. After the coordinates of the neutrons have been generated they are also given random unit direction vectors.

The next step is to calculate the determinants for the solution of equation (4.7). To simplify the calculation the coefficients are substituted as shown in equation (4.8). Us- ing these coefficients the determinants can be calculated for each neutron–particle pair to determine the particles that are hit. This can be used for calculations using a single neutron and those using a plethora of neutrons as it creates a determinant matrix with the required amount of columns for the determinant vectors of each neutron.

With the cell indexing and projectile collision detection it is not necessary to limit the

(30)

number of neutrons used in the calculation as the number of particles chosen for hit detection is significantly smaller than the total. It is the smallest with the projectile collision detection as only the ones actually hit are chosen. The downside of this approach is the need to simulate each neutron separately.

5.4 Free path calculation

The calculated determinants are used to find the least non-complex free path lengths for the neutrons from the determinant matrix using the equation (4.9). If only complex values are received the neutron does not collide with a fuel particle and has to be moved to the boundary. In the case of projectile collision detection only real determinants are received. The received mean free path lengths are saved in a free path matrix and the distances travelled by the non-colliding neutrons in a travelled free path matrix.

If enabled, the random periodic boundary condition is used after arriving at the boundary.

This simulates the different possible orientations the fuel pebbles can have in comparison to each other. If the distance travelled is zero, the coordinates of the neutron have been randomized during the previous iteration and the new direction has been out from the pebble. In this case the coordinates and direction are resampled.

5.5 Observed packing factor

In order to determine the packing factors observed by neutrons passing through the material, another script is to be written. This is done to help determine the effects the possible difference between the actual packing density and that experienced by a point object traversing through the region has on the mean free path and neutron multiplication factor.

The calculation is performed by randomly choosing a position of a virtual particle on the boundary of the generated space, sampling a random direction inwards and tracking the regions the virtual particle goes through. The regions are tracked by determining the intervals between the collisions with fuel particles, determined by solving the surface equations of the particles. After the intervals have been calculated, the observed packing factor is determined by calculating the ratio of the distance travelled through fuel particles and the total distance travelled.

(31)

6 RESULTS

This chapter includes the results received from the performed calculations. Reference measurements are performed using Serpent. Mean free path simulations are performed with and without BCs using the linear and segmented collision methods. Extended volume and projectile collision simulations as well as the simulations of the experienced packing factor are performed without a BC.

6.1 Serpent reference results

The reference results received from Serpent can be seen in table 6.1. Comparison calculations were run using packing factors ranging from 0.01 to 0.20 to compare the effects of the increase of packing factor. The listed results are calculated analogically. The corresponding absolute errors are listed in table 6.2.

Table 6.1. Neutron multiplication factors with f_p from 0.01 to 0.20.

Packing factor f_p Explicit Implicit Packing factor f_p Explicit Implicit

0.01 0.89828 0.89873 0.11 1.04223 1.04531

0.02 1.09350 1.09344 0.12 1.02139 1.02486

0.03 1.15314 1.15320 0.13 1.00120 1.00527

0.04 1.16695 1.16712 0.14 0.98174 0.98621

0.05 1.16132 1.16173 0.15 0.96299 0.96819

0.06 1.14688 1.14750 0.16 0.94517 0.95066

0.07 1.12796 1.12919 0.17 0.92771 0.93386

0.08 1.10680 1.10851 0.18 0.91125 0.91764

0.09 1.08529 1.08749 0.19 0.89528 0.90227

0.10 1.06366 1.06594 0.20 0.87999 0.88762

(32)

Table 6.2. Neutron multiplication factor errors with fp from 0.01 to 0.20.

Packing factor f_p Explicit Implicit Packing factor f_p Explicit Implicit 0.01 ±0.00015 ±0.00015 0.11 ±0.00013 ±0.00013 0.02 ±0.00013 ±0.00013 0.12 ±0.00014 ±0.00014 0.03 ±0.00012 ±0.00012 0.13 ±0.00014 ±0.00014 0.04 ±0.00012 ±0.00012 0.14 ±0.00014 ±0.00014 0.05 ±0.00012 ±0.00012 0.15 ±0.00014 ±0.00014 0.06 ±0.00012 ±0.00012 0.16 ±0.00015 ±0.00014 0.07 ±0.00012 ±0.00013 0.17 ±0.00015 ±0.00015 0.08 ±0.00013 ±0.00013 0.18 ±0.00015 ±0.00015 0.09 ±0.00013 ±0.00013 0.19 ±0.00015 ±0.00015 0.10 ±0.00013 ±0.00013 0.20 ±0.00015 ±0.00015

6.2 Linear and segmented collision simulations

The linear free path length script was used to generate results from one million neutrons using a packing factor f_p of 0.093. From these results the means, medians and modes are determined. These results are compared to those received from the reference model.

The confidence intervals are calculated at the significance of 4.75σ, that is the confidence interval covers approximately 99.9999 % of the population. The results are compiled in table 6.3.

Table 6.3. The results of the linear and segmented collision simulations with f_p of 0.093, compared with the expected results.

Description Symbol Unit Linear Segmented Expected Neutrons n_n [−] 1.0·10⁶ 1.0·10⁶ 1.0·10⁶

Mean free path ¯l [cm] 0.596 0.598 0.595

Error _a [cm] ±0.003 ±0.003 ±0.003

_r [%] ±0.499 ±0.499 0.475

Significance α [−] 4.750σ 4.750σ 4.750σ

Median ˜l [cm] 0.407 0.407 0.412

Min l_min [cm] 9.221·10⁻⁷ 3.255·10⁻⁷ 4.612·10⁻⁷

Max lmax [cm] 9.835 8.786 7.7356

Additionally the distributions for the previous results are plotted in figures 6.1 and 6.2

(33)

using histograms with a bin size of 0.02 cm.

Figure 6.1. The free path length distribution of the linear script.

Figure 6.2. The free path length distribution of the expected results.

(34)

The mean free paths of neutrons with varying packing factors was tested between packing factors of 0.01 and 0.20. The mean free paths were calculated using 200,000 results. These results are shown in table 6.4.

Table 6.4. The results from the linear program compared with the reference results.

Packing factor Reference Mean free path Error Median

f_p [-] l₀ [cm] ¯l [cm] [cm] ˜l [cm]

0.01 6.072 6.189 ±0.043 4.276

0.02 3.005 3.046 ±0.021 2.087

0.03 1.983 2.014 ±0.014 1.388

0.04 1.472 1.491 ±0.010 1.025

0.05 1.165 1.185 ±0.008 0.812

0.06 0.961 0.972 ±0.007 0.667

0.07 0.815 0.824 ±0.006 0.562

0.08 0.705 0.715 ±0.005 0.489

0.09 0.620 0.624 ±0.004 0.424

0.10 0.552 0.553 ±0.004 0.378

0.11 0.496 0.497 ±0.004 0.339

0.12 0.450 0.451 ±0.003 0.306

0.13 0.410 0.409 ±0.003 0.278

0.14 0.377 0.375 ±0.003 0.255

0.15 0.348 0.347 ±0.002 0.236

0.16 0.322 0.319 ±0.002 0.215

0.17 0.299 0.298 ±0.002 0.202

0.18 0.279 0.277 ±0.002 0.188

0.19 0.261 0.259 ±0.002 0.174

0.20 0.245 0.242 ±0.002 0.164

(35)

6.3 Projectile collision and packing factor simulations

One million neutrons were simulated, with a packing factor of 0.093, to compare the linear and segmented scripts with the projectile collision method. The results are listed in table 6.5. Additionally, the non-leaked free path length distribution can be seen from figure 6.3.

Table 6.5. Projectile collision results, using af_p of 0.093, compared to the reference model.

Description Symbol Unit Projectile collision Reference model

Neutrons n_n [−] 1.0·10⁶ 1.0·10⁶

Leaked neutrons n_l [−] 39 −

Mean free path

¯l [cm] 0.587

0.595

Leakless mean free path 0.587

Error _a [cm] ±0.003 ±0.003

_r [%] ±0.470 ±0.475

Leakless error _a [cm] ±0.003 ±0.003

r [%] ±0.468 ±0.475

Significance α [−] 4.750σ 4.750σ

Median

˜l [cm] 0.409

0.413

Leakless median [cm] 0.408

Min lmin [cm] 1.151·10⁻⁴

2.641·10⁻⁷

Leakless min 1.151·10⁻⁴

Max l_max [cm] 10.470

8.163

Leakless max 7.491

(36)

Figure 6.3. Projectile collision method free path length distribution with f_p of 0.093.

The projectile collision detection method was also tested with a packing factor of 17 % using 1,000,000 neutrons. These results were then compared to an equal amount of randomly generated mean free paths using equation (4.3). The results are shown in table 6.6. The comparison plot is shown in figure 6.4. The packing factor experienced by the neutron was sampled 400,000 times and these results are included in table 6.6.

Table 6.6. Projectile collision results, using af_p of 0.17, compared to the reference model.

Description Symbol Unit Projectile collision Reference model

Neutrons n_n [−] 1.4·10⁶ 1.4·10⁶

Mean free path ¯l [cm] 0.292 0.299

Error _a [cm] ±0.001 ±0.001

_r [%] ±0.205 ±0.210

Significance α [−] 4.750σ 4.750σ

Observed packing factor fp,o [−] 0.1729 0.1700

Median ˜l [cm] 0.204 0.207

Min l_min [cm] 3.509·10⁻⁷ 2.492·10⁻⁸

Max l_max [cm] 4.279 4.407

(37)

Figure 6.4. Projectile collision method and reference model results with a packing factor of 0.17.

In order to determine possible differences between the distributions of neutrons leaving from the surfaces of fuel particles and those generated in the graphite matrix, both cases were simulated 400,000 times. Based on these results, the combined mean free path was calculated. The results are listed in table 6.7 and they are plotted in figure 6.5.

Table 6.7. A comparison of results between generating the neutrons on the surface and off the surface of fuel particles, with a f_p of 0.017, using

the projectile collision method.

Description Symbol Unit On surface In matrix Combined

Neutrons n_n [−] 4·10⁵ 4·10⁵ 8·10⁵

Mean free path ¯l [cm] 0.293 0.284 0.288

Error _a [cm] ±0.002 ±0.002 ±0.002

_r [%] ±0.733 ±0.755 ±0.526

Significance α [−] 4.750σ 4.750σ 4.750σ

Median ˜l [cm] 0.205 0.196 0.201

Min l_min [cm] 1.753·10⁻⁴ 5.426·10⁻⁸ 5.426·10⁻⁸

Max l_max [cm] 3.883 3.679 3.883

(38)

Figure6.5.Projectilecollisionmethodfreepathlengthdistributioncomparisonbetweenneutronsgeneratedoff(top)oron(middle) thesurfaceofaparticle,andthecombineddistribution(bottom),withfpof0.17.Therelativeproportionsoftheoffandon distributionsisnotknownforarealcaseandtheycanthusbedifferent.

(39)

7 EVALUATION

The reference results received with Serpent, seen from the table 6.1, using the explicit model show an initial increase in the multiplication factor as a function of the packing factor. It peaks at a packing factor of 0.04, after which starts to constantly decrease. This peak is caused by the presence of an ideal mixture of fuel and moderator. Below 0.04 the amount of moderator present is larger than necessary, reducing the amount of interactions between neutrons and fuel. Above it the increased presence of U-238, combined with the reduction of the amount of moderator, causes more neutrons to be captured and less thermal neutrons to be produced and viable for fission.

A similar trend can be seen with the use of the implicit model. The initial increase of the neutron multiplication factor peaks at 4 %, after which it declines. The results differ noticeably from the explicit models from a packing factor of 7 % upwards. The error relative to the explicit results increase on average by 0.06 % for every 1 % of packing factor, rising from 0.11 % at a 7 % packing factor to 0.87 % at 20 %. These errors correspond to reactivity differences of 97 pcm and 977 pcm, respectively. The decreasing trend in the multiplication factors after the peak suggests that an increase in the packing factor causes a decrease in the neutron multiplication factor and, conversely, a decrease in the packing factor leads to an increase in the neutron multiplication factor. This implies that the implicit model observes a lower packing factor than is actually used. This would lead to a conclusion that the current model’s interaction mean free path, or the mean particle distance, given by the equation (4.2), is too large for packing factors above 7 %.

The results from the linear and segmented programs, seen from table 6.3 are greater than the reference model’s interaction mean free path but lie within the confidence interval. The neutrons were generated outside of the fuel particles and were found to follow an exponential distribution similar to the reference model, as can be seen from figures 6.1 and 6.2. This would suggest that the current model is still accurate enough around a packing factor of 0.093. The comparison between the reference results and the linear program, seen from table 6.4, suggests that the model produces too short interaction mean free paths up to a packing factor of 0.12, after which they are too large. This would, given the previously suggested explanation, mean that the neutron multiplication factors below a packing factor of 0.04 should be lower for the implicit method, due to the apparent lower packing density. Above 0.04 the multiplication factor would thus also be lower up to 0.12 due to an apparent increase in packing factor. Past 0.12, the neutron multiplication factor would be greater for the implicit model due to the lessened effect of the increase of

Simulation of the Interaction Mean Free Path between Neutrons and TRISO Particles