Theory and Methodology
1. General Program Flow
The software conducts a collision-driven Monte Carlo simulation, where it simulates each particle individually. Prior to running the simulation, it must undergo initialization, which involves the following steps:
- Interpolating the total, absorbtion, and scattering cross section data.
- Interpolating data of the initial energy distribution of the neutrons.
- Sampling the initial neutron energies from the energy distribution.
- Setting the initial position and energies of the neutrons.
- Setting the position and size of the tank.
The initial spectrum is discussed in Initial Energy Spectrum, the method of sampling from a distribution is described in Monte Carlo Sampling, and the interpolation of distributions in Interpolation. After the initializing, the actual simulation of the diffusing neutrons can start. Each neutron's energy and position is updated consecutively which consists of the following steps:
- Give the neutron a random direction.
- Check if the neutron is inside the tank, if false end the loop else continue.
- Get the total cross section corresponding to the current energy of the neutron.
- Determine the nuclide with which the neutron collides.
- Check if the the neutron has a thermal energy.
- Compute the mean free path of the neutron and sample the distance to the next collision.
- Determine the type of reaction. If it is absorption, end the loop; if it is elastic scattering, proceed.
- Sample a direction from the angular distribution and compute the associated energy loss fraction.
- Update the Neutron's position and energy.
- Repeat steps 2-9 until the neutron either escaped or is absorbed.
Each step is discussed in more detail in the following sections:
| Steps | Section |
|---|---|
| 2 | Geometry |
| 3-4 | Cross Sections |
| 5 | Thermal Motion |
| 6 | Distance To Next Collision |
| 7 | Neutron Interactions |
| 8-9 | Handling Scattering Interaction |
In the simulation, positions and energies of the neutrons are stored. Several properties can be computed and measured from this data post-simulation. The procedure of measuring these properties is described in Measuring Quantities.
2. Neutron Physics
2.1 Initial Energy Spectrum
The neutron source used is a mixture of americium-241 and beryllium-9. The following reactions take place in the source:
The first reaction involves the \( \alpha- \) decay of the americium. The \( \alpha \) particles have a probability of about \( 10^{-4}\) to undergo fusion with the berylium, described by the second reaction. This reaction leads to the emmission of neutrons with energies between 3 and 10 MeV. We use the Am-Be energy spectrum (Figure 1) from the ISO 8529-1 standard for sampling the initial energies of the neutrons.
The neutron spectrum is interpolated to serve as a probability distribution of energies \( P(\text{Energy}) \). This interpolation and subsequent sampling follow the methodology described in Interpolation and Monte Carlo Sampling respectively. The interpolated initial neutron energy spectrum from the Am-Be source is shown in Figure 2.
2.2 Cross Sections
The cross section indicates the likelihood of a particular interaction occurring between two particles. It is influenced by two main factors: the types of particles involved and their energies. In our scenario, there are two possible combinations of particles in an interaction: a neutron with an H-nucleus and a neutron with an O-nucleus. Unfortunately, it's not possible to calculate the cross section as a function of energy analytically. However, there is an abundance of freely accessible nuclear data available to everyone. In the program, cross section data from ENDF/B is used an interpolated following the procedure outlined in Interpolation.
In various reactions, different cross sections are involved. To proceed, we employ the following notations:
| Notation | Description |
|---|---|
| \(\sigma_{H, tot}, \sigma_{O, tot}\) | Total cross section for the specified nuclide. |
| \(\sigma_{H, s}, \sigma_{O, s}\) | Scattering cross section for the specified nuclide. |
| \(\sigma_{H, a}, \sigma_{O, a}\) | Absorption cross section for the specified nuclide. |
| \(\Sigma_{H_2O}, \Sigma_{H}, \Sigma_{O}\) | Total macroscopic for the specified nuclide or molecule. |
The theory behind these reactions is described in detail in Neutron Interactions. For each type of cross section and nuclide there are separate datafiles. Total macroscopic cross section for a molecule can be computed as follows:
where \(N_i\) is the atomic density of atoms of type \(i\):
where \(\rho\) is the density of the compound, \(M\) the molecular weight, \(N_a\) Avogrado's constant, and \(n_i\) the number of atoms of type \(i\) in the molecule. Thus, the total macrocscopic cross section for \(H_2O\) is given by:
The cross sections in barns (\(10^{-28} m^2\)) as a function of neutron energy for the O and H nuclide are shown in Figure 3 and 4.
2.3 Distance to Next Collision
The distance to the next collision \(l\) follows a distribution function given by:
where \(\Sigma\) is the total macroscopic cross section of an arbitrary material. Integrating both sides to some distance \(l\) gives:
As the probability distribution is normalized, we can set the this cummulative distribution to a uniformly distributed random number \(\xi\) on the interval \( [0, 1) \). Isolating \(l \) and using that in this case \(\Sigma \equiv \Sigma_{H_2O} \), gives us a formula for sampling the distance to the next collision:
where \(\lambda \) is the mean-free-path which is given by
The mean-free-path as a function of neutron energy is shown in Figure 5.
2.4 Neutron Interactions
To determine the type of reaction involving a nuclide, we initially identify which type of nuclide the neutron interacts with, either \( H\) or \( O\). The probability of the neutron interacting with H is determined by:
If \(\xi < P(H)\), where \(\xi\) is another uniformly distributed random number on the interval \( [0, 1) \), the neutron interacts with an \( H \) nuclide; otherwise, it interacts with an \( O \) nuclide. Since the cross sections depend on the energy of the neutron, \(P(H) \) also varies accordingly. The dependence of \( P(H) \) on neutron energy is depicted in Figure 6.
Once the nuclide of interaction is determined, it is necessary to consider the reactions that can take place. Firstly we check absorbtion. The probability for absorbtion is given by:
where \( \sigma_{i, tot} \) and \( \sigma_{i, a} \) are the total and absorbtion cross section of the previously identified nuclide respectively. Again, if \(\xi < P(\text{absorbtion})\), then the neutron is absorbed by the nuclide. Absorption results in the complete dissipation of the neutron's energy and we proceed to simulate the next neutron. The absorption probability \(P(\text{absorbtion})\) as a function of neutron energy for H and O are shown in Figure 7 and 8 respectively.
If the criteria mentioned above is not met, we consider a scattering interaction. This interaction leads to a reduction in the neutron's energy, determined by the outgoing angle. This angle is sampled from a distribution specific to the nuclide in question. The process of sampling the angle and computing the neutron's new direction and energy is elaborated upon in the Section Handling Scattering Interaction.
2.5 Handling Scattering Interaction
A scattering interaction involving a neutron and nuclide results in an energy loss, which depends on the outgoing angle of the neutron. At lower (thermal) neutron energies, there is notable relative motion between the neutron and the nuclide. The treatment of this scenario is discussed in the Section Thermal Motion. However, at higher energies, it is a reasonable approximation to consider the nuclides as stationary. We can extract the outgoing angles from datasets containing angular distributions sourced from ENDF/B. These distributions, which vary with increasing neutron energies, are defined by Legendre polynomials and are provided in the center-of-mass frame. The representation of these angular distributions using Legendre polynomials is detailed in the ENDF-6 Formats Manual, which we will briefly explain in Angular Distributions. With this angle the new energy and direction of the neutron can be computed as described in Transforming Coordinates.
2.5.1 Angular Distributions
The angular distributions depend on the energy \(E\) in the labframe and the cosine of the outgoing angle in the center-of-mass frame \(\mu\). The distributions are normalized:
We assume that the angular distributions of the scattered neutrons have azimuthal symmetry. This implies that we can uniformly select the azimuthal angle \(\phi\) from the range \( [0,2\pi) \). The distributions \( f(\mu, E) \) are represented as Legendre polynomials:
where \( l\) is the order of the Legendre polynomial, \(a_l \) the \(l^{th}\) Legendre polynomial coefficient, and \( NL \) the number of Legendre polynomials that represent the distribution. Below in Figure 9 and 10 the angular distributions are plotted for varying neutron energy.
The neutron energy distributions are discretely spaced for each incoming neutron energy \(E_i \), requiring a selection between the distribution corresponding to the energy immediately above or below the current neutron energy. This is done by defining an interpolation factor:
where \( E \) is the incoming energy of the neutron. We designate \(k\) as the selected distribution index. If \( \xi > f \), then \(k=i\); otherwise, \(k=i+1\), where \( \xi \) is a uniform distributed random number on the interval \( [0, 1) \). The cosine of the scattering angle \(\mu \) in the center-of-mass (CM) frame is sampled from the angular distribution corresponding to \(E_k\).
2.5.2 Transforming Coordinates
From the cosine of the scattering angle in the CM frame \(\mu \), the secondary velocity of the neutron can be calculated analytically. We start by considering the velocity of the center-of-mass which is given by:
where A is the mass of the target nuclide, \( \mathbf{v}_n \) the initial velocity of the neutron, and \( \mathbf{v}_t \) is the velocity of the target nuclide. For energies significantly surpassing thermal energy, \( \mathbf{v}_t \) can be considered negligble. The velocity of the neutron can now be transformed to the CM frame:
where the capital letter \( \mathbf{V} \) denotes the CM frame. If we express the post-collision direction in the CM frame as \( \mathbf{\Omega}' = (u', v', w') \), the components are given by:
where \( \phi \) is the azimuthal angle and \( \theta = \cos^{-1}(\mu) \). We assume azimuthal symmetry, so \( \phi \) can be sampled from a random uniformly distributed interval \( [0, 2\pi) \). Finally, transforming back to the lab frame results in the new velocity of the neutron:
The relation between the post- and pre-collision neutron energy in the labframe can be computed as
and the relation between the old and new position is given by:
where \( l \) is obtained following the methodology outlined in Distance to Next Collision.
2.6 Thermal Motion
When the energy of the neutron approaches thermal energies, the target nuclide cannot be assumed to be at rest. The target nucleus will have motion associated with its thermal vibration. Thermal motion has an effect on the cross sections and the secondary angular distribution. In The OpenMC Monte Carlo Code a method is described for sampling the velocity of the target nucleus to be used in the elastic scattering kinematic equations described in Handling Scattering Interaction.
For the activation process we are mainly interested in the high energy neutrons. For simplicity in the simulation, thermal neutron energies are drawn from a Maxwell-Boltzmann distribution. Furthermore, the resulting direction of a neutron from a scattering reaction is chosen to be random. A neutron is considered thermal when it reaches an energy below the threshold of \( 10 K_B T \approx 0.252 \ eV\), where \( K_B \) is Boltzmann's constant and \(T \) is the temperature of the water. The Maxwell-Boltzmann distribution for \( T = 293 K \) is used in the simulation and shown in Figure 11.
3 Geometry
The geometry of the experiment is shown in Figure 12. In the simulation the neutron source is placed at the origin. The water tank of radius \( R\) and height \( h \) is placed at position \( \mathbf{r}_t\) such that the neutrons diffuse from the correct initial position.
A Neutron at position \( \mathbf{r}_n \) is inside the tank if the following conditions are true:
where the notation \( || \mathbf{r} ||_{xy} = \sqrt{r_x^2 + r_y^2} \) is used. If these conditions are not met, the neutron is flagged as escaped and the simulation proceeds with simulating the next neutron.
4 Monte Carlo Sampling
In the initilization the energies of the neutrons are sampled from the distribution shown in Figure 2. The secondary angles of the neutrons in a scattering interaction are sampled from angular distributions on the fly. Some examples of these angular distributions are shown in Figure 10. Finally, when a neutron reaches the thermal energy threshold, its energies are sampled from the Boltzmann distribution shown in Figure 11. The Monte Carlo method used in the simulation is rejection sampling. For simplicity, the proposal distribution is chosen to be a uniform distribution on the interval \([0, 1) \) for each of the sampling processes. The following steps are performed to sample a point from a PDF \(f\):
- Sample y uniformly on the interval \([0, 1) \) and \( x\) between the bounds of the \( f \).
- Compute \( P = f(x) \).
- If \( y \leq P \) stop the process and return \(x\), else go back to step 1.
Figure 13 illustrates a large sample obtained from an angular distribution using rejection sampling. This figure is meant to demonstrate the sampling method, as in the simulation, only one secondary angle is sampled per elastic collision.
5 Interpolation
In the initialization of the simulation, data for the cross sections and initial energy distribution is interpolated. This means we estimate values that lie between known data points. It's important for cross sections to be continuous, ensuring that for every neutron energy, there's a corresponding cross section. Likewise, we need the initial energy distribution to be interpolated so we can generate an initial spectrum of neutron energies using the sampling method described in Monte Carlo Sampling.
We accomplish this interpolation using the Akima1DInterpolator from scipy.interpolate. This tool fits smooth curves between the data points by using cubic polynomials, ensuring a continuous representation of the data.
6 Measuring Quantities
For neutron activation, we are interested in the flux and the energy spectrum of the neutrons at varying radius in the tank. These quantities can be computed from the positions and energies of the neutrons at each collisions site in the tank, which are stored after running the simulation.
6.1 Flux
To compute the flux at a specific distance \(r\) from the origin, we consider a spherical shell of radius \(r\) and denote the position at the \(i_{th}\) collision site of a neutron as \(\mathbf{r}_{n_i}\). If the line between two consecutive positions \(\mathbf{r}_{n_i}\) and \(\mathbf{r}_{n_{i+1}}\) passes through the spherical shell, a counter is increased by one. Note that the neutron can pass through the shell multiple times. If we repeat this process for all neutrons, the flux \(\Phi\) at a distance \(r\) from the origin is given by:
where \(N\) is the number of times the neutrons pass through the spherical shell.
6.2 Energy Spectrum
For the energy spectrum at a distance \(r\) from the origin, we consider the same spherical shell. If the line between \(\mathbf{r}_{n_i}\) and \(\mathbf{r}_{n_{i+1}}\) passes through the spherical shell, we save the corresponding energy \(E_i\) of the neutron in a list. We repeat this process for all neutrons, appending the energies to the same list. Note that neutron can contribute multiple energy values to this list if it passes through the spherical shell multiple times. The resulting list is the energy spectrum at a distance \(r\) from the origin.