In this notebook, we demonstrate how to normalize tally values in OpenMC assuming you have a reactor operating at a known power. We'll begin by creating a very simple "reactor" model consisting of a sphere of U235 surrounded by water.
import openmc
u235 = openmc.Material()
u235.add_nuclide('U235', 1.0)
u235.set_density('g/cm3', 10.0)
water = openmc.Material()
water.add_components({'H': 2.0, 'O': 1.0})
water.add_s_alpha_beta('c_H_in_H2O')
fuel_sph = openmc.Sphere(r=6.5)
outer_sph = openmc.Sphere(r=15.0, boundary_type='vacuum')
fuel = openmc.Cell(fill=u235, region=-fuel_sph)
moderator = openmc.Cell(fill=water, region=+fuel_sph & -outer_sph)
geometry = openmc.Geometry([fuel, moderator])
model = openmc.Model(geometry=geometry)
model.settings.batches = 100
model.settings.inactive = 10
model.settings.particles = 1000
model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())
Let's run OpenMC to see what we get for $k_\text{eff}$:
sp_path = model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
keff = sp.keff
print(f'k-effective = {keff}')
k-effective = 0.9943+/-0.0030
Now let's say we want to know the fission reaction rate in [fissions/sec] in the fuel. We'll create a tally for this using the "fission" score:
fission_tally = openmc.Tally()
fission_tally.filters = [openmc.CellFilter([fuel])]
fission_tally.scores = ['fission']
model.tallies = [fission_tally]
We then run OpenMC and get the result from the statepoint:
sp_path = model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
tally = sp.tallies[fission_tally.id]
fission_rate = tally.mean.ravel()[0]
print(f'Fission rate = {fission_rate} [reactions/source]')
Fission rate = 0.3955335940378975 [reactions/source]
So, the number we get here is not quite what we wanted—it's the fission reaction rate per source particle. We want the reaction rate per second. In order to determine that, we need to figure out the source rate, i.e., the number of source particles emitted per second. The method for doing this is described in the OpenMC user's guide but here we'll give a practical demonstration.
We need two things in order to determine our desired normalization factor in units of [source/sec]. First, we need to know the power in [W] or [J/sec]. This is not something OpenMC or any other neutronics code will provide for you—that is something you have to provide as a user. In addition to the power, we need the observed heating rate from the simulation in [J/source]. To get that, we can add another tally for the heating rate. We'll use the fission-q-recoverable score but there are some pitfalls as to which score you choose (see further discussion below).
heating_tally = openmc.Tally()
heating_tally.scores = ['fission-q-recoverable']
model.tallies = [fission_tally, heating_tally]
sp_path = model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
tally = sp.tallies[heating_tally.id]
heating_rate_ev = tally.mean.ravel()[0]
print(f'Heating rate = {heating_rate_ev} [eV/source]')
Heating rate = 76235017.82874997 [eV/source]
The OpenMC tally gives us units of [eV/source] and we need to convert this to [J/source]:
joule_per_ev = 1.60218e-19
heating_rate = heating_rate_ev * joule_per_ev
print(f'Heating rate = {heating_rate} [J/source]')
Heating rate = 1.2214222086486662e-11 [J/source]
At this point, let's assume that our spherical reactor has a power of 30 kW. That is:
power = 30e3 # [J/s]
The tally normalization factor is simply the power in [J/s] divided by the heating rate in [J/source]:
source_per_sec = power / heating_rate
Whenever we want our tallied quantities in units of "per sec" rather than "per source neutron", all we have to do is multiply by the normalization factor. For example, our fission rate would be:
fission_rate_sec = fission_rate * source_per_sec
print(f'Fission rate = {fission_rate_sec:.3e} [reactions/sec]')
Fission rate = 9.715e+14 [reactions/sec]
As a quick sanity check, we can multiply by 200 MeV/reaction to see how much heating we'd expect from this many fission reactions:
expected_power = fission_rate_sec * 200e6 * joule_per_ev
print(f'Expected power = {expected_power:.1f} [W]')
Expected power = 31130.1 [W]
The value is fairly close to the expected 30 kW but not quite exact since we used a crude estimate of 200 MeV/reaction. The actual amount of heat being deposited per reaction can be determined by dividing the heating rate by the fission rate.
mev_per_reaction = heating_rate_ev / fission_rate * 1e-6
print(f'Heat per reaction = {mev_per_reaction:.1f} [MeV]')
Heat per reaction = 192.7 [MeV]
How much heat you observe will depend on which score you choose to estimate the heating. OpenMC has a variety of scores that allow you to tally heating that differ in subtle ways, so you should exercise some caution when picking a score. A few notes about the available scores:
"heating" — This score gives you the heating from all reactions, not just fission, so it gives you a way of accounting for heat from (n,$\gamma$) reactions. It will also include photon energy deposition but only if you are running a coupled neutron–photon calculation; for a neutron-only calculation, you should use the "heating-local" score to account for photon energy deposition."heating-local" — This score is the same as "heating" except that it assumes photon energy is deposited locally and thus can be used in a neutron-only calculation. This score requires special cross section data (MT=901) that is not present in all data libraries, particularly those that were converted straight from ACE files. The "official" data libraries from https://openmc.org do have the necessary data included."fission-q-recoverable" — This score represents the recoverable energy release from fission, which is basically the total fission energy release less neutrinos. Do note that it does not account for heating from non-fission reactions."kappa-fission" — This older score is basically equivalent to "fission-q-recoverable" but does not include an energy-dependence on the incident neutron energy.