Thermal Neutron diffusion: Monte Carlo approach

Here some easy program for the calculation of the flux of thermal neutron in medium like water, emitted from point and planar source. The procedure consist in following several neutron sampling randomly the distance of the next collision from the actual neutron position. Then all the collision coordinates are recorded and used to evaluate the flux.

Sampling

in analog monte carlo method the idea is to sample the distance of the next collision will happen from the following distribution: pdf = \Sigma e^{-\Sigma d} if we integrate the pdf from zero to infinity we obtain the cumulative distribution function CDF = 1 - e^{-\Sigma d} the CDF have the property to be uniformly distributed, so it’s possible to extract a random value RAND and obtain d from it.

In this case i need an uniform scattering, so the angle \phi is sampled uniformly from zero to 2*Pi as the cosine of \theta so doing the emission direction result to be uniformly distributed. Below you can see the Matlab iplementation of the position updating:

            phi = rand*2*pi;
            mu = 2*rand -1;             %cosine of theta
            theta = acos(mu);
            d = - 1/SIGMA*log(1-rand);  %sampling the distance
            %update the position
            X = X + d*sin(theta)*cos(phi);
            Y = Y + d*sin(theta)*sin(phi);
            Z = Z + d*mu;

In case of multiple type of interaction ( i.e. scattering, absorption, fission ..) or homogeneus mixture of materials you have to evaluate which interaction take place, a script can be found in iterated source method page. Anyway the idea is to compare a random number with the ratio of the macroscopic cross section of interaction “a” over the total macroscopic cross section \Sigma_{tot} = \sum_a^N \Sigma_a, then if RAND < \frac{\Sigma_a}{\Sigma_{tot}} interaction “a” take place. If not we check RAND < \frac{\Sigma_a + \Sigma_{a+1}}{\Sigma_{tot}} and so on until we find the correct event.

Point source

Show matlab code

in these simulations i have take in account elastic scattering and absorption and i have neglected inelastic scattering. in the comparison that you can see below i have used the following analytical solution for the equation of diffusion:

\phi(r) = {S e^{-R/L} \over 4\pi D R}

S: source strength

R: radius

L: diffusion length

D: diffusion parameter

Planar source

Show matlab code

In this case the analytical solution is:

\Phi(x) = {{SL \over 2D}e^{-x/L}}