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

Validation and Verification
1

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