Negative MT=301 Macroscopic cross-sections for Inconel

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:
inconel_301_xs

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

2 Likes

This is, unfortunately, a well-known problem with nuclear data evaluations. I would recommend reading this article which discusses the issue in general as well as many specific cases (including the Mo isotopes in JEFF 3.3).

1 Like

Very interesting article, though the proposed solutions seem kind of out of reach, at least for me.

Is there of way of telling OpenMC to ignore negative kermas? Otherwise I think I’ll just set all negative values to zero in the *.h5 data files for now.

I think the OpenMC data repo has a script for generating CENDL data in the openmc format, according to the article you posted CENDL-3.2 provides correct data for Mo kerma. I guess a better approach so as not to underestimate inconel heating would be to replace the faulty MT301 data in the JEFF-3.3 files with the CENDL values.

Thanks!

At present, no. So as you suggested, the easiest fix is to set the negatives to zero directly in the .h5 files. I’m definitely open to adding an option to set negatives to zero so users don’t have to muck around with the data, although I’d first like to get a sense of what other codes do (and why).

1 Like