I. M. Sobol – The Monte Carlo Method – Section 07 – Computation of Neutron Transmission Through a Plate

The probabilistic laws of interaction of an individual elementary particle (neutron, photon, meson, etc.) with matter are known. Usually it is needed to find the macroscopic characteristics of processes in which the number of participating particles is enormous, such as densities, fluxes and so on.

This situation is quite similar to that encountered in Section 5 and Section 6, and offers a very convenient object for the Monte Carlo calculations. The neutron physics is probably the field where the Monte Carlo method is used more frequently than anywhere else. We shall consider the simplest version of the problem of neutron transmission through a plate.

Statement of the Problem: Let the flux of neutrons with energy Eo be incident on a homogeneous infinite plate 0 ≤ x ≤ h. The angle of incidence is 900. Upon collision with atoms of the plate material neutrons may be either elastically scattered or captured. Let us assume for simplicity that energy of neutrons is not changed in scattering and that any direction of “recoil” of a neutron from an atom is equally probable (this is sometimes the case in neutron collisions with heavy atoms).

This would mean that the neutron could be absorbed or reflected or it simply passes through the plate, depending on the thickness of it. Given that the neutron that enters could perpetually stay inside the place, we define here that if the neutron happens to stay inside the plate for the 10 successive collisions, it is absorbed. If it is on the left side of the place (the direction that it originally entered through) before it’s tenth collision, and if it is on the right side the plate before it completes 10 collisions, then it is regarded as passed.


Simulation

The simulation of this problem was done in two phases. All the code has been presented at the end.

  • Phase:1 involved the simulation of the neutron paths using Python
  • Phase:2 involved the plotting of the phase 01 data using Wolfram Mathematica

Reading the images:

  • There are three sets of images shown below
    • the first one with the lines in blue shows the neutrons that have passed the plate
    • the second one with the lines in red shows the neutrons that got reflected
    • the third one with the lines in green shows the neutrons that got absorbed in the plate
  • Each image that you see has two parts to it which essentially represent the same data
    • the first part on the top is simply a straight line that passes through all the point that the neutron has visited
    • the second part of the bottom is the B-Spline Curve that uses the same set of points as the plot that is above it. This is just a smoothed out representation of the path of the neutrons and does not represent the actual paths

Visualization

 

Statistics

  • Each of the following charts below shows the same data in three different formats.
  • We see the statistics for the Neutrons that have
    • passed the barrier
    • reflected from the barrier
    • absorbed from the barrier
  • There are 1000 data points in each of the chart and the number of neutrons per data point is shown on the top of the chart
    • example 10000 Neutrons means that each point has 10000 neutrons

 


Code

Python Code: This runs the simulation

from mathematica.random_functions import *
from mathematica.monitoring import *
import csv

for t in range(0, 1000):
    _filePath = "D:\Mathematica Files 4K\probability_problems\i_m_sobol_the_monte_carlo_method" \
                "\i_m_sobol_the_monte_carlo_method_section_7"
    _exportFileName = _filePath + "\\i_m_sobol_the_monte_carlo_method_section_7_temp_" + str(t) + ".csv"

    TimeTagMessage('Opening file ' + _exportFileName)
    __file = open(_exportFileName, 'w')
    TimeTagMessage('Opening writer object')
    __writer = writer = csv.writer(__file, delimiter=',')
    TimeTagMessage('Writing rows to the file')

    for i in range(0, 1000):
        _seed = RandomReal(0, 1)
        _seedList = [_seed]
        while 0 <= _seed <= 2 and _seedList.__len__() <= 10:
            _seed += RandomReal(-1, 1)
            _seedList.append(_seed)
        writer.writerow(_seedList)
        if (i + 1) % 1000 == 0:
            TimeTagMessage("\t\tWriting row # " + str(i + 1) + " to the file")

    TimeTagMessage('Closing the file ' + _exportFileName)
    __file.close()

Mathematica Code: This plots the python exported data

ClearAll[yRandomizer];
yRandomizer[list_List] :=
    Module[{yRand = RandomReal[{-.02, .02}]}, {{-1, yRand}} ~ Join ~
        MapThread[{#1, #2} &, {list, {yRand} ~ Join ~
            RandomReal[{-1, 1}, Length@list - 1]}]]

Module[{data =
    DeleteCases[
      Import[StringReplace[NotebookFileName[], ".nb" -> ".csv"]], {}],
  collisions, last, reflected, absorbed, passed, plotDataReflected,
  plotDataPassed, plotDataAbsorbed},
  collisions = {Length[#], Last[#]} & /@ data;
  last = Last /@ data;
  reflected = Select[data, (Last[#] < 0) &];
  absorbed = Select[data, (0 <= Last@# <= 2) &];
  passed = Select[data, (Last@# > 2) &];
  plotDataReflected = yRandomizer /@ reflected;
  plotDataPassed = yRandomizer /@ passed;
  plotDataAbsorbed = yRandomizer /@ absorbed;

  Echo[Length /@ {passed, reflected, absorbed}];

  Framed@Column[Graphics[
    {
      {Red, Opacity[0.1], Arrowheads@0.02,
        Thickness@0, #[plotDataReflected]},
      {Gray, Opacity[0.04], Arrowheads@0.02,
        Thickness@0, #[plotDataPassed]},
      {Gray, Opacity[0.04], Arrowheads@0.02,
        Thickness@0, #[plotDataAbsorbed]},
      {InfiniteLine[{{0, -1}, {0, 1}}]},
      {InfiniteLine[{{2, -1}, {2, 1}}]}
    }
    , PlotRange -> All, ImageSize -> 788,
    PlotLabel ->
        "Reflected Neutrons: " <> ToString@# <>
            " representation"] & /@ {Arrow, BSplineCurve}]
]

(*Need to run python code before running this one*)
Module[{multiSimData,
  printTemporary = PrintTemporary["Beginning analysis"]},
  multiSimData =
      Module[{data = DeleteCases[Import[#], {}], reflected, absorbed,
        passed},
        reflected = Select[data, (Last[#] < 0) &];
        absorbed = Select[data, (0 <= Last@# <= 2) &];
        passed = Select[data, (Last@# > 2) &];
        NotebookDelete[printTemporary];
        printTemporary =
            PrintTemporary[Style["Analyzing the file " <> #, Darker@Green]];
        Length /@ {passed, reflected, absorbed}
      ] & /@
          SortBy[FileNames["*temp*", NotebookDirectory[]],
            ToExpression@Last[StringSplit[StringDelete[#, ".csv"], "_"]] &];
  Export[StringReplace[NotebookFileName[], ".nb" -> "_1000.xlsx"],
    multiSimData]
]

Column@Module[{data =
    Import[StringReplace[NotebookFileName[], ".nb" -> "_1000.xlsx"]][[
        1]]}, data =
    Association@
        MapThread[#1 -> #2 &, {{"Passed", "Reflected", "Absorbed"},
          Transpose[data]}];
DistributionChart[data, ImageSize -> 500,
  ChartLabels ->
      Placed[{"Mean = " <> ToString@# & /@ (Mean /@ Values@data),
        Keys@data}, {Above, Below}], ChartElementFunction -> #,
  PlotLabel ->
      "1000 Neutron Statistics for Barrier Width = 2 units & \
Absorption number = 10"] & /@ {"PointDensity", "SmoothDensity",
  "Density"}
]

End of the post 😉


.

.

.

.

.

.

.

.

.

.

.

.

.

.