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
`
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
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
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())
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 ` 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]:
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():
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))
`