Hi everyone,
while trying to compute energy deposition in a rod a Inconel, I noticed that the “heating” tally sometimes ended up being negative. This shows up in the openmc-plotter as the tally not being fully drawn:
To investigate this phenomenon I plotted the macroscopic MT=301 cross section of inconel using the following script:
import openmc
import matplotlib.pyplot as plt
import h5py
import numpy as np
data_path = "/data/nuclear_data/jeff33/HDF5/"
def get_301_xs(isot, temperature=900):
with h5py.File(f"{data_path}/{isot}.h5") as f:
energy = f[f"{isot}/energy/{temperature}K"][...]
xs = f[f"{isot}/reactions/reaction_301/{temperature}K/xs"][...]
return energy, xs
def compute_301_macro(mat, filter="all"):
energies = {}
xs = {}
for nuclide in mat:
e_, xs_= get_301_xs(nuclide)
if filter=="all":
energies[nuclide], xs[nuclide] = get_301_xs(nuclide)
elif filter=="noMo":
if nuclide[:2] != "Mo":
energies[nuclide], xs[nuclide] = get_301_xs(nuclide)
elif filter=="Mo":
if nuclide[:2] == "Mo":
energies[nuclide], xs[nuclide] = get_301_xs(nuclide)
E = np.unique(np.concatenate(list(energies.values())))
maxi = min([l[-1] for l in energies.values()])
E = E[E<maxi]
XS = sum([np.interp(E, energies[isot], xs[isot]) for isot in energies])
return E, XS
if __name__=="__main__":
mass_fraction = {"Ni": 0.58, "Cr": 0.215, "Mo": 0.09, "Fe": 0.05, "Nb": 0.0365,
"Co": 0.01, "Si": 0.005, "Mn": 0.005, "Ti": 0.004, "Ag": 0.004,
"C": 0.001, "S": 0.01, "P": 0.001,}
inconel = openmc.Material()
inconel.set_density = 8.44
for k, v in mass_fraction.items():
inconel.add_element(k, v*100, percent_type="wo")
inconel_cc = inconel.get_nuclide_atom_densities()
fig, ax = plt.subplots(1, 1)
ax.plot(*compute_301_macro(inconel_cc, filter="all"), label=r"$\kappa_{301}(\mathrm{Inconel})$")
ax.plot(*compute_301_macro(inconel_cc, filter="Mo"), label=r"$\kappa_{301}(\mathrm{Mo})$")
ax.plot(*compute_301_macro(inconel_cc, filter="noMo"), label=r"$\kappa_{301}(\mathrm{Inconel} \ \backslash \ \{\mathrm{Mo}\})$")
ax.set_xscale("log")
ax.set_yscale("symlog")
ax.set_title("MT=301 Macroscopic cross-sections of Inconel625")
ax.set_xlabel("Energy [eV]")
ax.set_ylabel(r"MT=301 Cross section [$\mathrm{barn}\cdot \mathrm{eV}$]")
ax.grid()
ax.legend()
plt.savefig("MT301_inconel.png")
This script plots the following image:
It seems that some nuclides (in particular isotopes of molybdenum) in the inconel material have negative MT=301 heating cross sections, and that these negative cross-sections are so large that the overall macroscopic heating cross-section of inconel becomes negative for energies below ~100KeV.
I am having some troubles interpreting these negative cross-sections and the negative heating tally results we are getting. Surely this doesn’t mean that inconel cools down when subject to a neutron flux? Should we set the heating tally of molybdenum to zero and only keep the positive contribution to heating?
Any help would be very welcome!
Thanks!
Nicolas