On the Deviation of Results in Using OpenMC to Produce Multi group Cross Sections for Openmoc in the Presence of Vacuum Boundaries

[wind zhou]

I need to use OpenMC to create a cross-section library of a single three-dimensional fuel rod for use by OpenMOC. When all boundary conditions are reflection, I can use the cross-section library produced to calculate good results in OpenMOC. However, when the upper and lower surface boundary conditions are vacuum, the error reaches several thousand pcm. Can someone tell me the reason for this.Here is the code,
'

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-dark')

import openmoc

import openmc
import openmc.mgxs as mgxs
import openmc.model as model
import openmc.data
from openmc.openmoc_compatible import get_openmoc_geometry


# create a model object to tie together geometry, materials, settings, and tallies
#model = openmc.model()
model = model.Model()


# 1.6% enriched fuel
fuel = openmc.Material(name='3.1% Fuel')
fuel.set_density('g/cm3', 10.257)
fuel.add_nuclide('U234', 6.11864E-06)
fuel.add_nuclide('U235', 7.18132E-04)
fuel.add_nuclide('U236', 3.29861E-06)
fuel.add_nuclide('U238', 2.21546E-02)
fuel.add_nuclide('O16', 4.57642E-02)

# borated water
water = openmc.Material(name='Water')
water.set_density('g/cm3', 0.743)
water.add_nuclide('H1', 4.96224E-02)
water.add_nuclide('B10', 1.07070E-05)
water.add_nuclide('B11', 4.30971E-05)
water.add_nuclide('O16', 2.48112E-2)

# zircaloy
gap = openmc.Material(name='gap')
gap.set_density('g/cm3', 0.0001786)
gap.add_nuclide('He4', 2.68714E-05)

# zircaloy
zircaloy = openmc.Material(name='clad')
zircaloy.set_density('g/cm3', 6.56)
zircaloy.add_nuclide('Zr90', 2.18865E-02)
zircaloy.add_nuclide('Zr91', 4.77292E-03)
zircaloy.add_nuclide('Zr92', 7.29551E-03)
zircaloy.add_nuclide('Zr94', 7.39335E-03)
zircaloy.add_nuclide('Zr96', 1.19110E-03)
zircaloy.add_nuclide('Sn112', 4.68066E-06)
zircaloy.add_nuclide('Sn114', 3.18478E-06)
zircaloy.add_nuclide('Sn115', 1.64064E-06)
zircaloy.add_nuclide('Sn116', 7.01616E-05)
zircaloy.add_nuclide('Sn117', 3.70592E-05)
zircaloy.add_nuclide('Sn118', 1.16872E-04)
zircaloy.add_nuclide('Sn119', 4.14504E-05)
zircaloy.add_nuclide('Sn120', 1.57212E-04)
zircaloy.add_nuclide('Sn122', 2.23417E-05)
zircaloy.add_nuclide('Sn124', 2.79392E-05)
zircaloy.add_nuclide('Fe54', 8.68307E-06)
zircaloy.add_nuclide('Fe56', 1.36306E-04)
zircaloy.add_nuclide('Fe57', 3.14789E-06)
zircaloy.add_nuclide('Fe58', 4.18926E-07)
zircaloy.add_nuclide('Cr50', 3.30121E-06)
zircaloy.add_nuclide('Cr52', 6.36606E-05)
zircaloy.add_nuclide('Cr53', 7.21860E-06)
zircaloy.add_nuclide('Cr54', 1.79686E-06)
zircaloy.add_nuclide('Hf174', 3.54138E-09)
zircaloy.add_nuclide('Hf176', 1.16423E-07)
zircaloy.add_nuclide('Hf177', 4.11686E-07)
zircaloy.add_nuclide('Hf178', 6.03806E-07)
zircaloy.add_nuclide('Hf179', 3.01460E-07)
zircaloy.add_nuclide('Hf180', 7.76449E-07)

# Instantiate a Materials collection
model.materials = openmc.Materials([fuel, gap ,water, zircaloy])

# Create cylinders for the fuel and clad
fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.4096)
gap_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.418)
clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.475)

# Create box to surround the geometry
#box = openmc.model.rectangular_prism(1.26, 1.26, boundary_type='reflective')
min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')
max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')
min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')
max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')
min_z = openmc.ZPlane(z0=-182.88, boundary_type='vacuum')
max_z = openmc.ZPlane(z0=+182.88, boundary_type='vacuum')
 
# Create a Universe to encapsulate a fuel pin
pin_cell_universe = openmc.Universe(name='1.6% Fuel Pin')

# Create fuel Cell
fuel_cell = openmc.Cell(name='1.6% Fuel')
fuel_cell.fill = fuel
fuel_cell.region = -fuel_outer_radius
pin_cell_universe.add_cell(fuel_cell)

# Create a gap Cell
gap_cell = openmc.Cell(name='1.6% Clad')
gap_cell.fill = gap
gap_cell.region = +fuel_outer_radius & -gap_outer_radius
pin_cell_universe.add_cell(gap_cell)
# Create a clad Cell
clad_cell = openmc.Cell(name='1.6% Clad')
clad_cell.fill = zircaloy
clad_cell.region = +gap_outer_radius & -clad_outer_radius
pin_cell_universe.add_cell(clad_cell)

# Create a moderator Cell
moderator_cell = openmc.Cell(name='1.6% Moderator')
moderator_cell.fill = water
moderator_cell.region = +clad_outer_radius
pin_cell_universe.add_cell(moderator_cell)

# Create root Cell
root_cell = openmc.Cell(name='root cell')
root_cell.fill = pin_cell_universe

# Add boundary planes
root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z

# Create root Universe
root_universe = openmc.Universe(universe_id=0, name='root universe')
root_universe.add_cell(root_cell)
# Create Geometry and set root Universe
model.geometry = openmc.Geometry(root_universe)


# OpenMC simulation parameters
batches = 500
inactive = 50
particles = 10000

# Instantiate a Settings object
settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.output = {'tallies': False}

# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-1.26, -1.26, -1.26, 1.26, 1.26, 1.26]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)

# Activate tally precision triggers
settings.trigger_active = True
settings.trigger_max_batches = settings.batches * 4

model.settings = settings

# Instantiate a "fine" 8-group EnergyGroups object
fine_groups = mgxs.EnergyGroups([0., 0.058, 0.14, 0.28,
                                 0.625, 4.0, 5.53e3, 821.0e3, 20.0e6])



# Extract all Cells filled by Materials
openmc_cells = model.geometry.get_all_material_cells().values()

# Create dictionary to store multi-group cross sections for all cells
xs_library = {}

# Instantiate 8-group cross sections for each cell
for cell in openmc_cells:
    xs_library[cell.id] = {}
    xs_library[cell.id]['total']  = mgxs.TotalXS(groups=fine_groups)
    xs_library[cell.id]['transport']  = mgxs.TransportXS(groups=fine_groups)
    xs_library[cell.id]['fission'] = mgxs.FissionXS(groups=fine_groups)
    xs_library[cell.id]['nu-fission'] = mgxs.FissionXS(groups=fine_groups, nu=True)
    xs_library[cell.id]['nu-scatter'] = mgxs.ScatterMatrixXS(groups=fine_groups, nu=True)
    xs_library[cell.id]['chi'] = mgxs.Chi(groups=fine_groups)


# Create a tally trigger for +/- 0.01 on each tally used to compute the multi-group cross sections
tally_trigger = openmc.Trigger('std_dev', 1e-2)

# Add the tally trigger to each of the multi-group cross section tallies
for cell in openmc_cells:
    for mgxs_type in xs_library[cell.id]:
        xs_library[cell.id][mgxs_type].tally_trigger = tally_trigger


# Instantiate an empty Tallies object
tallies = openmc.Tallies()

# Iterate over all cells and cross section types
for cell in openmc_cells:
    for rxn_type in xs_library[cell.id]:

        # Set the cross sections domain to the cell
        xs_library[cell.id][rxn_type].domain = cell
       
        # Tally cross sections by nuclide
        xs_library[cell.id][rxn_type].by_nuclide = True
               
        # Add OpenMC tallies to the tallies file for XML generation
        for tally in xs_library[cell.id][rxn_type].tallies.values():
            tallies.append(tally, merge=True)
           
model.tallies = tallies

# Run OpenMC
sp_file = model.run()

# Load the last statepoint file
sp = openmc.StatePoint(sp_file)
#sp = openmc.StatePoint('statepoint.2000.h5')

# Iterate over all cells and cross section types
for cell in openmc_cells:
    for rxn_type in xs_library[cell.id]:
        xs_library[cell.id][rxn_type].load_from_statepoint(sp)


# Create an OpenMOC Geometry from the OpenMC Geometry
openmoc_geometry = get_openmoc_geometry(sp.summary.geometry)

# Get all OpenMOC cells in the gometry
openmoc_cells = openmoc_geometry.getRootUniverse().getAllCells()

# Inject multi-group cross sections into OpenMOC Materials
for cell_id, cell in openmoc_cells.items():
   
    # Ignore the root cell
    if cell.getName() == 'root cell':
        continue
   
    # Get a reference to the Material filling this Cell
    openmoc_material = cell.getFillMaterial()
   
    # Set the number of energy groups for the Material
    openmoc_material.setNumEnergyGroups(fine_groups.num_groups)
   
    # Extract the appropriate cross section objects for this cell
    transport = xs_library[cell_id]['transport']
    nufission = xs_library[cell_id]['nu-fission']
    nuscatter = xs_library[cell_id]['nu-scatter']
    chi = xs_library[cell_id]['chi']
   
    # Inject NumPy arrays of cross section data into the Material
    # NOTE: Sum across nuclides to get macro cross sections needed by OpenMOC
    openmoc_material.setSigmaT(transport.get_xs(nuclides='sum').flatten())
    openmoc_material.setNuSigmaF(nufission.get_xs(nuclides='sum').flatten())
    openmoc_material.setSigmaS0(nuscatter.get_xs(nuclides='sum').flatten())
    openmoc_material.setChi(chi.get_xs(nuclides='sum').flatten())

# Generate tracks for OpenMOC
track_generator = openmoc.TrackGenerator3D(openmoc_geometry, num_azim=16, num_polar=6,
                                           azim_spacing=0.1, z_spacing=0.75)
track_generator.generateTracks()

# Run OpenMOC
solver = openmoc.CPUSolver(track_generator)
solver.setNumThreads(12)
solver.computeEigenvalue()
solver.setConvergenceThreshold(1e-6)

# Print report of keff and bias with OpenMC
openmoc_keff = solver.getKeff()
#openmc_keff = sp.keff.n
openmc_keff = sp.k_combined.nominal_value
bias = (openmoc_keff - openmc_keff) * 1e5

print('openmc keff = {0:1.6f}'.format(openmc_keff))
print('openmoc keff = {0:1.6f}'.format(openmoc_keff))
print('bias [pcm]: {0:1.1f}'.format(bias))
'
The results of running the above code,
openmc_keff = 1.198678
openmoc_keff = 1.127340
bias [pcm]: -7133.8

[Guillaume Giudicelli]

Hello

One of the first things that differ between MC and MOC is the representation of neutron scattering. MOC in your input file is isotropic, MC is anisotropic considering Hydrogen in the moderator.

To get better agreement you can set OpenMC to run with isotropic scattering then compare again.

Then you can let OpenMC be anisotropic again but use the basic inscatter transport correction implemented in MGXS.

And finally you can use custom transport correction factors derived from infinite medium water calculations in the moderator.

This should improve that source of discrepancy.

Then you will have to look at the discretization as well, and finally resonance self-shielding effects.

I ll refer you to Dr Boyd's thesis (isotropic scattering, 2D), Dr Gunow's (3D), and my (3D + equivalence) thesis (available on MIT dspace) for some of these comparisons and convergence studies

Best,

Guillaume