fusion-energy / neutronics-workshop

A workshop covering a range of fusion relevant analysis and simulations with OpenMC, DAGMC, Paramak and other open source fusion neutronics tools
MIT License
117 stars 51 forks source link

should we have an example for photonuclear neutron production from photons #229

Open shimwell opened 1 year ago

shimwell commented 1 year ago

the code is ease enough if this is useful we can add an example

# coupled photon transport can be enabled for a neutron transport
# with the settings.photon_transport = True
# that is fairly standard and we need to enable photon production and transport
# for tallies like heat, dose and it is important to transport the photon
# as it will not deposit its energy locally.
# Charged particles will deposit their energy locally so they don't need to be transported

# However this example is somewhat different
# This examples shows that a 4MeV photon source creates neutrons
# We do this by creating a photon source and tallying neutron flux

# You could try changing the material to a lower Z atom like Be or Li and the
# flux should go up partly because these shield photons to a lesser extent.

import openmc

# MATERIALS
material = openmc.Material()
material.add_element('Fe', 1)
material.set_density('g/cm3', 7)
my_materials = openmc.Materials([material])

# GEOMETRY
# surfaces
vessel_inner = openmc.Sphere(r=500, boundary_type='vacuum')
# cells
inner_vessel_region = -vessel_inner
inner_vessel_cell = openmc.Cell(region=inner_vessel_region)
inner_vessel_cell.fill = material
my_geometry = openmc.Geometry([inner_vessel_cell])

# SIMULATION SETTINGS
# Instantiate a Settings object
my_settings = openmc.Settings()
my_settings.batches = 10
my_settings.inactive = 0
my_settings.particles = 500
my_settings.run_mode = 'fixed source'
# Create a photon point source
my_source = openmc.Source()
my_source.space = openmc.stats.Point((0, 0, 0))
my_source.angle = openmc.stats.Isotropic()
my_source.particle = 'photon'  # default would otherwise be neutron
my_source.energy = openmc.stats.Discrete([4e6], [1])
my_settings.source = my_source
# added a cell tally for neutron flux production
cell_filter = openmc.CellFilter(inner_vessel_cell)
flux_tally = openmc.Tally(name='flux')
flux_tally.filters = [cell_filter]
flux_tally.scores = ['flux']
my_tallies = openmc.Tallies([flux_tally])

# Run OpenMC!
model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies)
sp_filename = model.run()

# open the results file
sp = openmc.StatePoint(sp_filename)

# access the tally using pandas dataframes
flux_tally = sp.get_tally(name='flux')

# printing the tally shows that we get a neutron flux as the tally is not 0
print(flux_tally.mean, flux_tally.std_dev)
# printed results are ...
# >>> [[[14.42136056]]] [[[0.14143424]]]
py1sl commented 1 year ago

Not certain but i think this should give a zero value, the photonuclear cross section threshold for Iron is above 10 MeV, Lead or Be are typical low threshold materials for photonuclear reactions. does the tally need a particle filter?