Iterated source method – only thermal neutrons

With this method is possible to evaluate the moltiplication factor k of a reactor. The idea is to start with a source of neutron uniformly distributed, for simplicty in this work the particles are monoenergetic (thermal) these neutrons will collide, escape, be absorbed or cause fission. when a fission event happen the position is collected for the next generation of neutrons.

Each generation consist in the same number of neutrons, this is equivalent to force criticality even if the system is sub/super-critical.

To force criticality we must change the average number of neutron emitted in fission event \bar \nu , so in each generation n \bar \nu_n = \bar \nu / k_n , where k_n is the moltiplication factor in generation n. Since the number of neutron in the system is constant whe have that \bar \nu_1 l_1 = \bar \nu_n l_n where l_1 and l_nare the number of fission event in generation 1 and n. It follw that k_n = k_1 l_n /l_1. And that’s all, what remain is the sampling procedure and the collision events.

In the following example i refer to an homogeneus reactor Uranium + Moderator (Graphite) with cylindrical shape.

SAMPLING UNIFORMLY IN A CYLINDER

Since the first generation is distributed uniformly (you can choose better distribution more similar to the real one) here a little script to perform this task (Matlab):

function [x,y,z] = sampling_cylinder(H,R,N)

z = rand([N, 1])*H-H/2;
theta = rand([N, 1])*2*pi;
r = sqrt(rand([N, 1]))*R;
y = r.*cos(theta);
x = r.*sin(theta);

end

Given in input the height, the radius of the cylinder and the number of neutron it returns a set of coordinates x,y,z for each neutron.

Note that z and \theta are sampled uniformly, r instead is r = R\sqrt{RAND}

CHECK IF NEUTRON IS INSIDE THE REACTOR

There are several way to do this, i decided to use an external function that given the neutron coordinate return a flag equal to 1 if the particle is out of the cylinder:

function [flag_out] = check_if_inside(pos_x,pos_y,pos_z,R,H)

%cylindrical coordinates

D = sqrt(pos_x^2 + pos_y^2);

if(abs(pos_z) >= H/2 || D >= R)         %if out from the system
    flag_out = 1;
else
    flag_out = 0;
end

end

The flag will be used in the main while cycle as an “exit condition”

FISSION AND ABSORPTION

The possible reaction are elastic scattering absorption and fission with U-235 U-238 and mooderator, but a neutron will be followed until it’s absorbed or cause fission, so we have to check for these events:

function [flag_fiss,flag_abs] = interaction()

flag_fiss = 0;
flag_abs = 0;

data                                       %import data for cross section

random = rand;
random1 = rand;
if random < macro_sigmaU5/macro_sigma       %interaction with U5
    if random1 < sigmaU5(1)/sigmaU5_tot     %fission 
        flag_fiss = 1;
    elseif random1 < (sigmaU5(1) + sigmaU5(2))/sigmaU5_tot  %absorption
        flag_abs = 1;
    end      
elseif random < (macro_sigmaU5 + macro_sigmaU8)/macro_sigma       %interaction with U8    
    if random1 < sigmaU8(1)/sigmaU8_tot     %fission 
        flag_fiss = 1;
    elseif random1 < (sigmaU8(1) + sigmaU8(2))/sigmaU8_tot  %absorption
        flag_abs = 1;
    end
else                                        %interaction with moderator    
    if random1 < sigmaM(1)/sigmaM_tot       %fission 
        flag_fiss = 1;
    elseif random1 < (sigmaM(1) + sigmaM(2))/sigmaM_tot  %absorption
        flag_abs = 1;
    end      
end

end

This part of the code is a little bit tricky but conceptually is very simple: “random” is a random number (from uniform distribution) used to select the collision nuclide (U-235 U-238 or C, must compare with macroscopic cross section) then random1 is used to select the type of collision. Take for exemple the variable sigmaU5: it’s a 3 component vector, it contain in order fission, absorption and scattering cross section sigmaU5_tot, instead, is the sum of the three (you can see it in data file). So for exemple knowing that the neutron will collide with U-235 if random1<sigmaU5(1)/sigmaU5_tot then the nuclide will fission.

In the end the function return two flags : flag_fiss = 1 mean that fission event happened, flag_abs = 1 mean that neutron is absorbed. Both flags equal to zero means elastic scattering.

You can find some easy examples of scattering in monte carlo here

SAMPLING FROM NEW FISSION POSITIONS

For the generations after the first one the neutron distribution is not uniform anymore, we have to sample from the fission coordinates of the previous generation. To do so we can create a vector of coordinater and then sample uniformly from it:

function [ID] = sampling_from_fission_vector(N,fission_number)

%sample N position of the fission vector
ID = round(rand([N,1])*(fission_number-1))+1;

end

so given the number of neutron N and the number of fission event we can sample uniformly the IDs in the fission coordinates vector, then the positions of the new particle will be something like (in the main cycle):

[ID] = sampling_from_fission_vector(N,fission_number);

    x    = fiss_pos_x(ID);
    y    = fiss_pos_y(ID);
    z    = fiss_pos_z(ID);

RESULTS

The average value of k is 1.07 against the theoretical one 1.09.

You can see that this script do not predict well the space distribution since some clustering is present. Big cores need a lot of particles, this clustering disappear with smaller geometries.

SOME FILES

Most of the feature of the code are described in this page, nevertheless could be useful to take a look at the complete set of files.

download (iterated_source.rar)