Polyethylene (HDPE) cross section for neutron moderation

I have an experiment with a He3 detector to detect D-D neutrons. The detector is encased in a HDPE block that I’m trying to model. Is it valid to use the following to define the material?

HDPE_material.add_elements_from_formula(“C2H4”, percent_type=“ao”)
HDPE_material.set_density(“g/cm3”, 0.95)

I guess as a broader statement more to do with the structure of materials, one wouldn’t use plain ole C12 to model a material made from graphite. So since the cross section data in openmc doesn’t include “polyethylene,” is it valid to use it’s chemical formula instead?

Yep that would make a material from the specified ratio of elements. The formula method you have there is the most concise way to do it. You could also use the material.add_element method to build up a material but I think in this case the formula is more convenient. You get the same end result as

HDPE_material.add_element(‘C’, 2)
HDPE_material.add_element(‘H’, 4)
HDPE_material.set_density(“g/cm3”, 0.95)

I looked to see if s alpha beta cross sections exist for polyethylene but couldn’t see any in ENDF/B8.0. So I think what you have is the best you can do

Here are some examples of making a material from elements if you want another option

Here is a collection of material definitions which includes a polyethylene neutronics_material_maker/pnnl_materials.json at 55fa45b6811c19a34874ab9f501c253a56df6051 · fusion-energy/neutronics_material_maker · GitHub

That’s very helpful thank you. I’ve also installed the neutronics material maker – what a great resource! I do see there’s thermal scattering data c_H_in_CH2 which could be what I’m looking for. Is there a way in openmc to view the thermal scattering cross sections of .h5 files to verify it is indeed correct?

I would reference this paper to confirm its accuracy.

I quite like this calculate_cexs function for getting the cross section data from materials.

The plot_xs function in that same file is also handy if you don’t want the data and just want a plot.

1 Like

Hi, if you are using the standard ENDF/B-VIII.0 nuclear data files, there is a normalization error in the file for polyethylene. But actually it is really easy to fix it using openmc.data:

ch2_path = '/home/student/tsl_software/endfb-viii.0-hdf5/neutron/c_H_in_CH2.h5'
data_h5 = openmc.data.ThermalScattering.from_hdf5(ch2_path)

for T in data_h5.temperatures:
    sigma_bound = data_h5.elastic.xs[T].bound_xs
    data_h5.elastic.xs[T].bound_xs = sigma_bound / 2.0
    print(f'Original sigma_bound: {sigma_bound} -> updated to: {(sigma_bound/2.0)}')

ch2_path_bak = ch2_path + '.bak'
os.system(f'mv {ch2_path} {ch2_path_bak}')
data_h5.export_to_hdf5(ch2_path)

This comes from one of the examples from the thermal scattering school at ESS:

1 Like

@Shimwell that function is very helpful – thank you!
@Ignacio wow what a great find – I wonder if anyone else has found similar discrepancies between thermal cross sections of NCrystals and ENDF? Are we sure that it’s not the NCrystal results that are off by a factor of two in the other direction?

@nickschw: because the bound scattering cross section of H1 is ~82 b. That is a know fact, but you can also see it here:

https://www.ncnr.nist.gov/resources/n-lengths/list.html

Aside of that, it is not just a comparison between codes: I am comparing with experimental data. There I am comparing with EXFOR entry # 30748002, measured by Dawidowski and Granada.

This is not an unusual normalization problem with these materials. I need to check if this is a problem in the source data used to generate the library, or in the processing scripts.

1 Like

As far as I can tell, the data that shows up in the HDF5 file (at least for the bound cross section) is exactly reproduced from the ENDF file itself, so it does not appear to be a processing issue (unless you’re talking about LEAPR→ENDF?)

Hi Paul. Sorry for hijacking the thread, maybe we should move the discussion somewhere else. I dig a bit and the problem is that some of the TSL evaluations with incoherent elastic component are not standard, and sigma_bound is multiplied by natom.

In particular, for ENDF/B-VIII.0 this affects polyethylene and lucite. A simple workaround is to check that the free atom scattering cross section derived from the bound atom scattering cross section in the elastic component matches the free atom scattering cross section from the inelastic component (or, if you have to divide by natom).

import openmc.data
import os
from io import StringIO

def get_file(url):
    user_agent = '--user-agent="Mozilla/5.0"'
    zipfile = 'file.zip'
    file = url.split('/')[-1].replace('zip', 'dat')
    os.system(f'wget --quiet -O {zipfile} {user_agent} "{url}"')
    os.system(f'unzip -o {zipfile}')
    return file

def read_file(file):
    data_endf = openmc.data.endf.Evaluation(file)

    file_obj = StringIO(data_endf.section[7, 4])
    items = openmc.data.endf.get_head_record(file_obj)
    items, values = openmc.data.endf.get_list_record(file_obj)

    msigma = values[0]
    awr = values[2]
    energy_max = values[3]
    natom = int(values[5])

    sigma_free = msigma/natom
    
    file_obj = StringIO(data_endf.section[7, 2])
    items = openmc.data.endf.get_head_record(file_obj)
    items, values = openmc.data.endf.get_tab1_record(file_obj)

    sigma_bound = items[0]
    
    return sigma_free, sigma_bound, awr, natom

endf8_urls = ['https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0005_H(YH2).zip',
              'https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0007_H(ZrH).zip',
              'https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0010_H(ice-Ih).zip',
              'https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0037_H(CH2).zip',
              'https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0039_H(Lucite).zip',
              'https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/tsl/tsl_0034_s-ch4.zip']

files = [get_file(url) for url in endf8_urls]

tol = 0.5 # tolerance in barns

print('\n Checking files:\n')

for fn in files:
    sigma_free, sigma_bound, awr, natom = read_file(fn)
    if (abs(sigma_free - sigma_bound/(1+1/awr)**2) < tol):
        print(f'{fn} is standard: {sigma_free} b aprox. eq. {(sigma_bound/(1+1/awr)**2)} b with natom = {natom}')
    elif (abs(sigma_free - sigma_bound/natom*1.0/(1+1/awr)**2) < tol):
        print(f'{fn} is non standard: {sigma_free} b aprox. eq. {(sigma_bound/(1+1/awr)**2)}/{natom} b')
    else:
        print(f'Something is wrong with {fn}: sigma_free= {sigma_free} b, sigma_bound = {sigma_bound} b, awr = {awr}, natom = {natom} ')

which gives:

tsl_0005_H(YH2).dat is standard: 20.43634 b aprox. eq. 20.43634000403785 b with natom = 1
tsl_0007_H(ZrH).dat is standard: 20.43634 b aprox. eq. 20.47800060779665 b with natom = 1
tsl_0010_H(ice-Ih).dat is standard: 20.505 b aprox. eq. 20.065000186502324 b with natom = 1
tsl_0037_H(CH2).dat is non standard: 20.43634 b aprox. eq. 40.872679858200634/2 b
tsl_0039_H(Lucite).dat is non standard: 20.4363375 b aprox. eq. 163.49071943280254/8 b
tsl_0034_s-ch4.dat is standard: 20.43634 b aprox. eq. 20.41903460041051 b with natom = 4