Error in the normalization of thermal elastic data

[marquezj]

I am moving the discussion we started in the forum with @yrrepy and @paulromano. The bug is that the elastic scattering cross section of several moderators (including CH2 and Lucite, but also other hydrocarbons in JEFF) are wrong in the OpenMC libraries.

To be clear, this on itself is not an OpenMC bug, but more likely proof that we (the thermal scattering evaluators) are horrible, horrible people, but let's try to fix the issue in OpenMC.

The physical fact is that for an incoherent scatterer there are two limits that need to be satisfied: for low energy the elastic cross section should go to the bound atom cross section value ($\sigma_b$), and for epithermal energy the inelastic cross section should go to the free atom cross section:

$$\sigma_{free} = \frac{\sigma_b}{(1+1/AWR)^2}$$

where AWR is the mass of the atom in neutron mass units.

For instance, for hydrogen the bound scattering cross section is 82 b, and the free atom cross section is 82 b / (1 + 1/0.99917)^2 = 20.5 b.

Now, thermal elastic and thermal inelastic processes are stored separately in the ENDF-6 files, and have independent (and sometimes not consistent) normalization cross sections.

For incoherent inelastic the value is stored in the first position of the B array, as a product of the free atom cross section times the number of principal atoms (which sometimes, but not always, is the number of that type of atoms in the molecule). The number of atoms is stored in the sixth position of the B array.

For incoherent elastic the value is stored in the "SB" scalar.

So, here is the part we are horrible people. The standard is:

$\sigma_{free} = B(1)/B(6)$

$\sigma_b = SB$

but in many evaluations $\sigma_b = SB/B(6)$, because that is the way NJOY treats the value in THERMR.

And that leads to the following: NJOY processing gives the right value but the (better) OpenMC incoherent elastic model does not:

E = np.geomspace(1e-4,5,1000)
data1 = openmc.data.ThermalScattering.from_njoy('n-001_H_001.endf', 'tsl_0039_H(Lucite).dat', use_endf_data=False)
data2 = openmc.data.ThermalScattering.from_njoy('n-001_H_001.endf', 'tsl_0039_H(Lucite).dat', use_endf_data=True)

xs = data1.elastic.xs['300K'](E)
plt.loglog(E,xs, 'k-', label='use_endf_data=False')

xs = data2.elastic.xs['300K'](E)
plt.loglog(E,xs, 'r--', label='use_endf_data=True')

plt.xlabel('Neutron energy [eV]')
plt.ylabel('Thermal elastic XS [b/atom]')
plt.legend()

One possible solution is to add a flag for this non standard data and process it accordingly.