Particle crossed surface but could not be located in any cell

User Support
1

Hi everyone,
I am new to OpenMC. I am attempting to simulate a HDPE neutron collimator with diameter 5cm, 18.5cm long, and aperture diameter 3.8cm. The neutron source energy is 2.45MeV (D-D reactions). However, I ran into a warning stating that after particle #number crossed surface 4 it could not be located in any cell and it did not leak. May I know how can I solve this problem? Below is my code. Thank you. Regards.

import openmc

Define the materials

hdpe = openmc.Material(name=“HDPE”)
hdpe.add_element(“C”, 2)
hdpe.add_element(“H”, 4)
hdpe.set_density(“g/cm3”, 0.95)
materials = openmc.Materials([hdpe])

Define the geometry

geometry = openmc.Geometry()

Define the collimator dimensions

outer_radius = 2.5 # cm
inner_radius = 1.9 # cm
length = 18.5 # cm
distance = 1.0 # cm

Create the collimator geometry

collimator_outer_surface = openmc.YCylinder(surface_id=1, r=outer_radius, boundary_type=‘vacuum’)
collimator_inner_surface = openmc.YCylinder(surface_id=2, r=inner_radius)
collimator_right_surface = openmc.YPlane(surface_id=3, y0=length / 2)
collimator_left_surface = openmc.YPlane(surface_id=4, y0=-length / 2)

collimator_region = -collimator_outer_surface & +collimator_inner_surface & -collimator_right_surface & +collimator_left_surface

collimator_cell = openmc.Cell(cell_id=1, region=collimator_region, fill=hdpe)
aperture_cell = openmc.Cell(cell_id=2, region=-collimator_inner_surface)
entry_cell = openmc.Cell(cell_id=3, region=-collimator_left_surface)
exit_cell = openmc.Cell(cell_id=4, region=+collimator_right_surface)

Add the cells to the geometry

geometry.root_universe = openmc.Universe(cells=[collimator_cell, aperture_cell, entry_cell, exit_cell])

Create the settings

settings = openmc.Settings()
settings.run_mode = “fixed source”
settings.particles = 1000000
settings.batches = 10
settings.inactive = 0
settings.source = openmc.Source(space=openmc.stats.Point(xyz=(0,-(length / 2) - distance,0)), energy=openmc.stats.Discrete([2.45e6], [1]))

Create the tallies

tally = openmc.Tally(name=“neutron_collimation”)
tally.scores = [“flux”]
tally.filters = [openmc.CellFilter(collimator_cell)]
tallies = openmc.Tallies([tally])

Run the simulation

model = openmc.model.Model(geometry, materials, settings, tallies)
model.run()

2

Hello. I recommend you add the max_lost_particle specification to your settings, with a limit of 1. Also, there is a particle hdf5 file created when a particle is lost. I do not recall if this is dependent on using the max_lost_particle spec. You can view data in the particle hdf5 file from the command line using: h5dump your_data_file.h5 | less

3

Hi,
@luluq279, I join @AlexandreTrottier about the issue you got there. It is defintely a problem with your geometry, mainly in the way you define boundary type of the other surfaces, since openmc take by default ‘transmission’ surface, unless you define another type. Usually, adjacent cell (filled with material) should exist, and I think it is not the case in your example:

collimator_outer_surface = openmc.YCylinder(surface_id=1, r=outer_radius, boundary_type=‘vacuum’)
collimator_inner_surface = openmc.YCylinder(surface_id=2, r=inner_radius) #<<<< here
collimator_right_surface = openmc.YPlane(surface_id=3, y0=length / 2) #<<<<< and here
collimator_left_surface = openmc.YPlane(surface_id=4, y0=-length / 2) # <<<<< and here

Besides that, I think that the definition of the following cells, should be reconsidered, since they are infinite cells and it is more likely causing scoring particles issue:

aperture_cell = openmc.Cell(cell_id=2, region=-collimator_inner_surface)
entry_cell = openmc.Cell(cell_id=3, region=-collimator_left_surface)
exit_cell = openmc.Cell(cell_id=4, region=+collimator_right_surface)

6

Hi sir. I have added more to my code but I ran into a warning: could not find the cell containing particle #some number, and also a runtime error: maximum number of lost particles has been reached. Would you mind sharing your idea on how to deal with this problem? Below is my code:

Define the materials

hdpe = openmc.Material(material_id=1, name=“HDPE”)

hdpe.add_element(“C”, 2)

hdpe.add_element(“H”, 4)

hdpe.set_density(“g/cm3”, 0.96)

lead = openmc.Material(material_id=2, name=“Pb”)

lead.add_element(‘Pb’, 1)

lead.set_density(‘g/cm3’, 11.34)

highvacuum = openmc.Material(material_id=3, name=‘high vacuum’)

highvacuum.add_nuclide(‘H2’, 2.0)

highvacuum.set_density(‘kg/m3’, 1e-7)

lowvacuum = openmc.Material(material_id=4, name=‘low vacuum’)

lowvacuum.add_nuclide(‘H2’, 2.0)

lowvacuum.set_density(‘kg/m3’, 1e-6)

feed_gas = openmc.Material(material_id=5, name=‘deuterium’)

feed_gas.add_nuclide(‘H2’, 2.0) # Deuterium

feed_gas.set_density(‘kg/m3’, 3.2e-4) # Density of deuterium gas at 4 mbar

materials = openmc.Materials([hdpe, lowvacuum, lead, highvacuum, feed_gas])

Define the collimator dimensions

outer_radius = 2.45 # cm

inner_radius = 1.9 # cm

length = 18.5 # cm

source_distance = 1.0 # cm

scintillator_distance = 3.4 # cm

chamber_width = 3.5 # cm

chamber_length = 15.0 # cm

Create the collimator geometry

collimator_outer_surface = openmc.XCylinder(surface_id=1, r=outer_radius, boundary_type=‘vacuum’)

collimator_inner_surface = openmc.XCylinder(surface_id=2, r=inner_radius, boundary_type=‘reflective’)

collimator_front_surface = openmc.XPlane(surface_id=3, x0=length / 2, )

collimator_back_surface = openmc.XPlane(surface_id=4, x0=-length / 2, )

scintillator_position = openmc.XPlane(surface_id=5, x0=-(length / 2) - scintillator_distance, boundary_type=‘vacuum’)

source_position = openmc.XPlane(surface_id=6, x0=length / 2 + source_distance, boundary_type=‘vacuum’)

chamber_wall_right = openmc.YPlane(surface_id=7, y0=chamber_width, boundary_type=‘vacuum’)

chamber_wall_left = openmc.YPlane(surface_id=8, y0=-chamber_width, boundary_type=‘vacuum’)

chamber_front = openmc.XPlane(surface_id=9, x0=-chamber_length, boundary_type=‘vacuum’)

chamber_back = openmc.XPlane(surface_id=10, x0=chamber_length, boundary_type=‘vacuum’)

Define the regions

collimator_region = -collimator_outer_surface & +collimator_inner_surface & -collimator_front_surface & +collimator_back_surface

aperture_region = -collimator_inner_surface & -collimator_front_surface & +collimator_back_surface

entry_region = +collimator_front_surface & -source_position & -collimator_outer_surface

exit_region = -collimator_back_surface & +scintillator_position & -collimator_outer_surface

empty_region = +collimator_outer_surface & -chamber_wall_right & +chamber_wall_left & -chamber_back & +chamber_front

lead_region = -collimator_outer_surface & -collimator_back_surface & +chamber_front & ~exit_region

Define the cells

collimator_cell = openmc.Cell(cell_id=1, region=collimator_region, fill=hdpe)

aperture_cell = openmc.Cell(cell_id=2, region=aperture_region, fill=lowvacuum)

entry_cell = openmc.Cell(cell_id=3, region=entry_region, fill=feed_gas)

exit_cell = openmc.Cell(cell_id=4, region=exit_region, fill=highvacuum)

empty_cell = openmc.Cell(cell_id=5, region=empty_region, fill=lowvacuum)

lead_cell = openmc.Cell(cell_id=6, region=lead_region, fill=lead)

Define the geometry and add the cells to the geometry

geometry = openmc.Geometry(openmc.Universe(cells=[collimator_cell, aperture_cell, entry_cell, exit_cell, lead_cell, empty_cell]))

Create the settings

settings = openmc.Settings()

settings.run_mode = “fixed source”

settings.particles = 1000000

settings.batches = 10

settings.inactive = 0

settings.source = openmc.Source(space=openmc.stats.Point((length / 2 + source_distance,0,0)), energy=openmc.stats.Discrete([2.5e6], [1]))

Create the tallies

tally = openmc.Tally(name=“neutron_collimation”)

tally.scores = [‘flux’, ‘absorption’, ‘elastic’, ‘scatter’, ‘total’, ‘events’]

tally.filters = [openmc.CellFilter([collimator_cell, aperture_cell, entry_cell, exit_cell, empty_cell, lead_cell])]

tallies = openmc.Tallies([tally])

Create the plot

plot = openmc.Plot()

plot.filename = “collimator”

plot.width = (35.0, 35.0)

plot.pixels = (500, 500)

plot.color_by = “material”

plot.colors = {hdpe: “green”, lowvacuum: (176, 196, 222), highvacuum: ( 47, 79, 79), lead: (165, 42, 42), feed_gas: ‘blue’}

plots = openmc.Plots([plot])

Run the simulation

model = openmc.model.Model(geometry, materials, settings, tallies, plots)
model.run()

Thank you very much sir. Regards.

7

HI @luluq279 ;

I took a look in your new model, and I tried some modifications mainly on the boundary type of the used surfaces (don’t define “vacuum” surface for cells which are adjacent to another material cells), and also some cells, like the empty_cell, which create an undefined volume on the right of the source surface.

What I did and it works is to put the source in the middle of the collimator (0, 0, 0). Of course, this is not the purpose of the model, but I did just to see if the issue came from the definition of the source, since there are too much lost particles.

I will let you know if I will find something new ;), unless you or someone else will find the origin of the problem !

Good luck

8

Hi sir. Thank you very much for your help. Now I am building a new model but without any collimator, only the detector.

Setup2
Setup21000×1000 2.25 KB

Attached below is the code:
setup_2.py (10.1 KB)

However, the result that I obtained is 0 flux at the detector (scintillator). If you run the code and save the result for mesh in excel format, you will notice starting at radius 6, there is no flux detected. May I know if there is anything wrong with the geometry defined?

9

Hi @luluq279 ,

First, please you can call me Salah and no need to call me Sir :wink:
Secondly, I looked in your python file and I got the same results as you : only the room chamber is registering neutrons and nothing elsewhere.
That why, I would to suggest you to do something simple in the first step … for example only your cylindrical source and detector, and put the whole system in spherical volume (with vacuum spherical boundary). Then, you can see if your detector cell is recording any of travelling neutrons or not. At this stage it will be much easier to find the issue with zero-tally value on your detector and to fix it.

regards

Salah

10

Dear Salah,
Thank you for the response. I have tried your suggestion, and I have finally figured out the problem. The boundary_type of the YCylinder, some YPlanes, and the body of the scintillator need to be changed to transmission so that flux reading will be non-zero. Thank you very much for your idea.

Regards,
Luqman

11

Dear Luqman;

It is very nice to hear that you could solve the issue, you are very welcome

Best regards

Salah