Photon flux with a layer of water and air

Hedi · June 20, 2024, 5:59pm

Hello everyone,
I am currently trying to simulate the photon flux through a layer of water and air to compare with my analytical calculation and check if my simulation is correct. I’m also trying to get to grips with the software and understand what I’m doing because I’ve only been using OpenMC for 1 month.
I’ve filled my box with nothing (empty) and just a layer of water and air but every time OpenMC measures the flux, it measures the flux in the empty space as if it couldn’t see the layer of water and air. And when I display my geometry and my universe, my layers are there. So I don’t understand, I’ve tried everything I can.
Can anyone help me?
Here is my code :

import openmc

import openmc.statepoint

import openmc.stats

import matplotlib.pyplot as plt

#################################################### Easy model (void) ####################################################################

##### Define material

air = openmc.Material(name='air')

air.add_element(element='N', percent=0.79, percent_type='ao')

air.add_element(element='O', percent=0.21, percent_type='ao')

air.set_density(units='g/cm3', density=1.293e-3)

water = openmc.Material(name='water')

water.add_elements_from_formula(formula='H2O')

#water.add_element(element='H', percent=2.0, percent_type='ao')

#water.add_element(element='O', percent=1.0, percent_type='ao')

water.set_density(units='g/cm3', density=1.0)

materials = openmc.Materials(materials=[air, water])

materials.export_to_xml()

##### Define the universe

pitch = 20.0

left = openmc.XPlane(x0=-pitch/2, boundary_type='vacuum')

right = openmc.XPlane(x0=pitch/2, boundary_type='vacuum')

bottom = openmc.YPlane(y0=-pitch/2, boundary_type='vacuum')

top = openmc.YPlane(y0=pitch/2, boundary_type='vacuum')

# Define a layer of water

water_left = openmc.XPlane(x0=0.)

water_right = openmc.XPlane(x0=2.0)

# Define a layer of air

air_left = openmc.XPlane(x0=2.0)

air_right = openmc.XPlane(x0=4.0)

##### Define the regions

water_box = +water_left & -water_right & +bottom & -top

air_box = +air_left & -air_right & +bottom & -top

box = +left & -right & +bottom & -top

##### Define the cell

water_cell = openmc.Cell(name='water_cell', region=water_box, fill=water)

air_cell = openmc.Cell(name='air_cell', region=air_box, fill=air)

box_cell = openmc.Cell(name='box_cell', region=box)

##### Define the universe

universe = openmc.Universe(name='universe', cells=[water_cell, air_cell, box_cell])

universe.plot(origin=(0.0, 0.0, 0.0), width=(10.0, 10.0), pixels=(1000, 1000), basis='xy', color_by='cell', colors={box_cell:'red', water_cell:'blue', air_cell:'green'},

legend=True, axis_units='cm')

#plt.show()

#quit()

##### Define the geometry

geometry = openmc.Geometry(root=universe)

geometry.export_to_xml()

##### Visualize the geometry

plot = openmc.Plot()

plot.filename = 'box_plot'

plot.width = (10.0, 10.0)

plot.pixels = (1000, 1000)

plot.color_by = 'cell'

plot.colors = {box_cell : 'red', water_cell : 'blue', air_cell : 'green'}

##### Visualize the 3D geometry

plot_3D = openmc.Plot()

plot_3D.filename = 'box_3D_plot'

plot_3D.type = 'voxel'

plot_3D.width = (10.0, 10.0, 10.0)

plot_3D.pixels = (100, 100, 100)

plots = openmc.Plots(plots=[plot, plot_3D])

plots.export_to_xml()

openmc.plot_geometry(output=True)

##### Define photon source

source_strenght = 10000 # Unit is source/cm/s in 2D and source/s in 3D

lower_source_coords = (-6.0, -1.0, 0.0)

upper_source_coords = (-4.0, 1.0, 0.0)

source_surface = (upper_source_coords[0]-lower_source_coords[0])*(upper_source_coords[1]-lower_source_coords[1]) # Unit is cm

source = openmc.IndependentSource()

source.space = openmc.stats.Box(lower_left=lower_source_coords, upper_right=upper_source_coords)

#source.space = openmc.stats.Point(xyz=(-5., 0., 0.))

source.angle = openmc.stats.Monodirectional(reference_uvw=[1, 0, 0])

source.energy = openmc.stats.Discrete(x=[30e3], p=[1.0])

source.time = openmc.stats.Uniform(a=0.0, b=1.0)

source.particle = 'photon'

source.strenght = source_strenght

##### Set the parameters of the simulation

settings = openmc.Settings()

settings.particles = 10000

settings.batches = 1

settings.source = source

settings.photon_transport = True # Activate the transport photon

settings.run_mode = 'fixed source' # Fixed source mode for photons

settings.export_to_xml()

##### Define a tally just in front of the source

# Define the mesh tally with a RegularMesh

mesh = openmc.RegularMesh(name='tally_mesh')

mesh_lower_coords = (5.0, -1.0)

mesh_upper_coords = (6.0, 1.0)

mesh.dimension = [100, 100] # Number of divisions in each direction

mesh.lower_left = mesh_lower_coords # Adjust the position at 10 cm of the source

mesh.upper_right = mesh_upper_coords # Adjust the position at 10 cm of the source

d = mesh_lower_coords[0] - upper_source_coords[0] # Tally-source distance

# Create a mesh filter

mesh_filter = openmc.MeshFilter(mesh)

##### Define a tally for photon flux

tally = openmc.Tally(name='photon_flux')

tally.scores = ['flux'] # Unit is particule.cm/source

tally.filters = [mesh_filter]

tallies = openmc.Tallies([tally])

tallies.export_to_xml()

##### Run the simulation

openmc.run()

# openmc.run(tracks=True) # Track every photons

##### Read the results

# Find the last statepoint file

statepoint_file = openmc.StatePoint('statepoint.1.h5')

##### Access the tally results

final_tally = statepoint_file.get_tally(name='photon_flux')

flux = final_tally.mean.flatten() # Unit is particule.cm/source

##### Print the total flux

total_flux = flux.sum() # Unit is particules.cm/source

print(f'Total photon flux: {(total_flux*source_strenght)} photons/source')

print(f'Total photon flux: {(total_flux*source_strenght)/(source_surface)} photons/cm2/source')

statepoint_file.close()

wahidluthfi · June 21, 2024, 8:40am

Hi Hedi,
I think the way you declare box region is not right.
I think you want to create a cell filled with nothing (vacuum) surrounding the water and air right? then to make sure that your box cell is not conflicted with the water and air cells, its region definition should look like

box = ((+left & -water_left) | (+air_right & -right)) & +bottom & -top

also, it is source.strength, not source.strenght, it might affect your source strength definition in your settings.xml.
then, regarding your mesh tally, since your system covers around -10 to 10 cm either in x and y axis, then your lower and upper coordinates should contain this whole system for better visualization, try

mesh_lower_coords = (-12.5, -12.5)
mesh_upper_coords = (12.5, 12.5)

I think you can check this example notebook about gamma detector from openmc, or maybe other examples that interest you.

github.com

openmc-dev/openmc-notebooks/blob/main/gamma-detector.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating the Signal of a Gamma Detector\n",
    "\n",
    "This IPython Notebook explains how to calculate the signal of a gamma detector using the pulse-height tally. \n",
    "It also shows how to create a realistic image of a detector simulation that takes into account the detector's intrinsic resolution.\n",
    "\n",
    "In particular, this Notebook illustrates the following features:\n",
    "* setting up a simulation of a gamma detector \n",
    "* defining a gamma source from a specific isotopic mixture\n",
    "* reprocessing the simulation output to obtain a realistic pulse-height spectrum"
   ]
  },
  {
   "cell_type": "code",

This file has been truncated. show original

Hedi · June 27, 2024, 2:59am

Hello wahidluthfi,
Thank you very much for your help.
I took a look at the gamma detector example you sent me and it gave me some ideas, but the way the tally is defined in the example is different from mine. As you can see, I’m using the RegularMesh class.
Which brings me to your comment, I understand that the tally must cover all the geometry for better visualisation but the value of my flux is different.
Let’s take the simple case where the box is filled with vacuum, in my code I’m simulating 10,000 particles so at the end I have to find my 10,000 particles, don’t I? If so, that’s not what I get. Knowing that the unit of flux is particles-cm/source, to find my original flux value do I have to multiply by the total size of the tally? Do I have to do the same if I then add my water and air layers?
I also had a question about the dimensions of RegularMesh because my aim is to display the flux as a function of x, so should I put [100, 1] in my dimensions, for example?

wahidluthfi · June 27, 2024, 7:10am

Hi Hedi,
I am not into photon transport since I am usually doing neutron transport for core modeling, so for my case (core modeling) the tally normalization factor can be found from the core thermal power.
But since photon transport tally outputs are dependent on source strength declared in simulation settings, then I think you could find other sources in discussion for photon tally normalization such as this

and some examples mentioned in this discussion might help you later

regarding this statement “Let’s take the simple case where the box is filled with vacuum, in my code I’m simulating 10,000 particles so at the end I have to find my 10,000 particles”
I think, since you talk about vacuum and monodirectional photons, then your photon source will not interact with any material, and it will move through space ‘to infinity and beyond’. hehe, I recommend you use air instead of vacuum since it feels more practical for me.

it will be different when you use a material since the interaction between photons and material will have various products: photoelectric, Compton scattering, and pair production which will change the total number of photons in the system. correct me if I’m wrong.

also, regarding the tally dimension, when using score: ‘flux’ on the mesh tally, I think the dimension is #-cm/source since when using mesh tally, openmc counts the number of particles in a specific mesh volume. so to bring back the flux into the typical flux dimension (#/cm2-sec), you will need to multiply this tally with the source and divide it by mesh volume, so #-cm/source. #/sec. 1/cm3. But that is what I understand when doing core calculations, I don’t know if it is also applicable to photon transport. somehow when tallying flux or power distribution, I only focus on the calculated trends or ratio, not the absolute value of each output.
documentation on tally: 8. Specifying Tallies — OpenMC Documentation

finally, regarding the mesh tally dimension, that’s right, and I think I could show how it could be done to plot along the x-axis with this notebook. but you will need to find other sources in other discussions regarding tally normalization for photon transport since I didn’t have any experience with it. Sorry Hedi
layerwater.ipynb (217.8 KB)

Wahid

Hedi · June 27, 2024, 2:30pm

Hello Wahid
Thank you very much for your reply.
It’s more or less what I was looking for but you’ve given me some interesting leads.
I’m going to try a few things and I’ll come back to you if I get stuck.
For the moment I’m just trying to get the same value and the same graph as this analytical formula: I = I0exp(-µx).
In order to see if my simulation is correct or not but for the moment it’s not there yet but I think the main problem comes from how I’ve defined my tally I think.

wahidluthfi · June 27, 2024, 9:03pm

If you want to look at how exp(-µx) works dramatically, then you are gonna need to use a thick layer of material with a high µ value for photons, like water, concrete, or lead since dry air does not attenuate photons properly. For modeling, You can use material composition provided by PNNL at
https://www.pnnl.gov/publications/compendium-material-composition-data-radiation-transport-modeling-1
Then I hope openmc calculation could be consistent enough to compare with an analytical approach like when using photon mass attenuation from NIST data

Since you are using 2cm of water on your script, then it only drops 1-exp(-µ×2cm) of the incoming photon with the energy you describe in input as µ also depends on photon energy.

Hedi · June 27, 2024, 9:30pm

Yes, you’re right, I’ve stopped using air and I use water, graphite, lead, etc…
As you said, I used the NIST site to get the µ values for each material, but every time I compare my theoretical calculations with the results of my simulation, I get totally different values.
On the other hand, what’s reassuring is that I’ve got the right curve shape but not the right values.
For example, I don’t have the right value at the beginning. If I simulate 10,000 particles, I should start at 10,000 but that’s not the case.
By trying to simulate my geometry in a vacuum and then using other materials, I’ve been able to understand certain things.
For example, by simulating 10,000 particles with a tally of dimension [100, 1] and size 20 cm, OpenMC calculates the flux on each cell as follows. It takes the 10,000 particles and divides them into 100, so that the tally counts 100 particles per cm, which it then multiplies by the size of the tally to give a total of 2,000 particles-cm/source. I don’t know if my analysis is right or not, but from what I’ve seen, that’s it.

wahidluthfi · June 28, 2024, 10:04am

Hi Hedi,
since you want to find the 10,000 particles that you declare as source strength, I am making a simple case that I think could help us to understand what openmc did in photon transport.
I am making a simple 3D geometry consisting of a small cube with a size 1x1x1 cm in front of a shield cell, also a cube 1x1x1 on its back. A photon source was placed at -5cm (point source monodirectional monoenergetic) with a source strength of 10,000 as you want.

when all material being used by all cells (cell 01, shield, cube 01, and cube 02) is vacuum or dry air, the photon distribution will look like this

I am using a mesh tally with a size of 1 cm in the x and y directions because I noticed that the 10,000 particle strength was reported at this size, so I think the source strength was defined as 10,000 per cc for this source definition. then, because there is almost no interaction in dry air (vacuum), each mesh in front of the source will report that the particle is 10,000. that’s why the total sum of particles on this mesh tally is 250,000 because there is 25 mesh in front of the source and each reported 10,000 photons. you also noticed that the photon rays are slightly shifted to y negative, that’s because the point source is counted along the mesh on the x-axis. If you move this source to be centered on the x-axis, then your 10,000 particles will be divided into 2, or 5,000 because it was counted in separate mesh bins, but the total photon counted will still be the same ~250,000.

then the 2 cube cells (1x1x1 cm) were used for tallying, so in a vacuum, everything is perfect on 10,000 as you want. if you change the size of the cube, i.e. become 2x2x2 cm, then this cell will count 2x10,000 since it covers 2 cm on the x-axis.

then the plot of photons along the x-axis also shows a 10,000 as you want

and when it comes to changing the 5cm cube of shield cell’s material, let’s say we want it to become water, then the plot will change, and the peak count also drops from 10,000 on cube 01 to around ~1000 on cube 02 on the back of the water shield.

I think there is something interesting about the absorption rate plotted here, but if I activate heating (photon heating), any tally in a vacuum will drop around 0, so yeah, I hope the expert on this topic could help to explain this behavior because I notice this when using vacuum. it is different when I use dry air as shown below.

then this is when we used dry air instead of a vacuum surrounding the shield

and last one is when we use lead which terminates almost all photons counted on cube 02 cell behind the lead shield

I think this notebook could be used if anybody wants to try it. I am also open for any correction
void.ipynb (262.8 KB)

Wahid

Topic		Replies	Views
Flux tally resulting in 0 for outer layers User Support	0	22	July 24, 2024
Need to understand photon interactions User Support	4	212	September 7, 2023
Photon flux discrepancies with EGSnrc, MCNP, and Serpent Validation and Verification	10	615	January 20, 2023
Photon flux with NaN value User Support	13	1465	June 7, 2021
Debugging a simple photon transport model generated with ChatGPT User Support	4	682	May 18, 2023

Photon flux with a layer of water and air

Related topics