Hi. I wanted to simulate a neutron collimator and find out how many particles it could collimate but it did not work. May I know what I should correct in the below 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) # Deuterium
highvacuum.set_density(‘kg/m3’, 1e-8) # Very low density
lowvacuum = openmc.Material(material_id=4, name=‘low vacuum’)
lowvacuum.add_nuclide(‘H2’, 2.0)
lowvacuum.set_density(‘kg/m3’, 1e-5)
feed_gas = openmc.Material(material_id=5, name=‘deuterium’)
feed_gas.add_nuclide(‘H2’, 2.0) # Deuterium
feed_gas.set_density(‘kg/m3’, 0.0003207) # Density of deuterium gas at 4 mbar
materials = openmc.Materials([hdpe, lead, highvacuum, lowvacuum, 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=‘vacuum’)
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,)
chamber_wall_right = openmc.YPlane(surface_id=7, y0=chamber_width)
chamber_wall_left = openmc.YPlane(surface_id=8, y0=-chamber_width)
chamber_front = openmc.XPlane(surface_id=9, x0=-chamber_length)
chamber_back = openmc.XPlane(surface_id=10, x0=chamber_length)
Define the regions
collimator_region = -collimator_outer_surface & +collimator_inner_surface & -collimator_front_surface & +collimator_back_surface
aperture_region = -collimator_inner_surface
entry_region = +collimator_front_surface & -chamber_back
exit_region = -collimator_back_surface & +scintillator_position
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, empty_cell, lead_cell]))
Create the settings
settings = openmc.Settings()
settings.run_mode = “fixed source”
settings.particles = 1000000
settings.batches = 20
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(entry_cell), openmc.CellFilter(aperture_cell), openmc.CellFilter(collimator_cell), openmc.CellFilter(exit_cell), openmc.CellFilter(empty_cell), openmc.CellFilter(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”, lead: (165, 42, 42), highvacuum: ( 47, 79, 79), lowvacuum: (176, 196, 222), feed_gas: ‘blue’}
plots = openmc.Plots([plot])
Run the simulation
model = openmc.model.Model(geometry, materials, settings, tallies, plots)
sp_filename = model. Run()
sp = openmc.StatePoint(sp_filename)
flux_tally = sp.get_tally(name=‘neutron_collimation’)
df = flux_tally.get_pandas_dataframe()
print(df)
Thank you. Regards.