crp/ecp-benchmarks
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Compare commits

base: crp:main
Branches Tags
crp:main
crp:fom-model
crp:fullcore-long
crp:fullcore-short
crp:ecp-assembly-short
crp:zlattice
crp:rings
..
compare: crp:zlattice
Branches Tags
crp:fom-model
crp:main
crp:fullcore-long
crp:fullcore-short
crp:ecp-assembly-short
crp:zlattice
crp:rings

4 commits
main ... zlattice

Author SHA1 Message Date
Paul Romano
2b175b4411 Fix subdivided fuel (position of lower z-plane, add outer universe) 2018-05-02 22:04:10 -05:00
Paul Romano
4af7da2ea7 Use shallow copy of material for differentiation 2018-05-02 22:03:40 -05:00
Paul Romano
7dbc5489d1 Improve z-lattice by using same universe in each level 2018-04-10 15:04:38 -05:00
Paul Romano
38bcb13900 Use lattice in Z to create universe for fuel pins 2018-04-10 14:45:31 -05:00
14 changed files with 725 additions and 963 deletions
Download patch file Download diff file Expand all files Collapse all files

186
smr/build-assembly-long.py
View file

@ -1,186 +0,0 @@
#!/usr/bin/env python3
import argparse
from math import pi, isclose
from pathlib import Path
import numpy as np
from tqdm import tqdm
import openmc
from smr.materials import materials, clone
from smr.surfaces import surfs, lattice_pitch, pin_pitch, bottom_fuel_stack, \
top_active_core, pellet_OR, active_fuel_length
from smr.pins import pin_universes, make_stack
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', dest='multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('--clone', action='store_true',
help='Clone materials for each cell instance')
parser.add_argument('--no-clone', dest='clone', action='store_false',
help='Do not clone materials for each cell instance')
parser.add_argument('-a', '--axial', type=int, default=100,
help='Number of axial subdivisions in fuel')
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(clone=False, multipole=True)
args = parser.parse_args()
# Make directory for inputs
if args.output_dir is None:
if args.depleted:
directory = Path('assembly-long-depleted')
else:
directory = Path('assembly-long-fresh')
else:
directory = args.output_dir
directory.mkdir(exist_ok=True)
rings = [0.1*pin_pitch, 0.2*pin_pitch]
assembly_long_surfs = [
surfs['bottom FR'],
surfs['bot active core'],
surfs['top active core'],
surfs['top pin plenum'],
surfs['top FR'],
surfs['bot upper nozzle'],
surfs['top upper nozzle']
]
univs = pin_universes(rings, args.axial, args.depleted)
fuel_univ = make_stack(
'Fuel (3.1%) stack no grid',
surfaces=assembly_long_surfs,
universes=[
univs['water pin'],
univs['end plug'],
univs['Fuel pin (3.1%) no grid'],
univs['pin plenum'],
univs['end plug'],
univs['water pin']
]
)
# Define the NumPy array indices for assembly locations where there
# may be CR guide tubes, instrument tubes and burnable absorbers
nonfuel_y = np.array([2,2,2,3,3,5,5,5,5,5,8,8,8,8,8,11,11,11,11,11,13,13,14,14,14])
nonfuel_x = np.array([5,8,11,3,13,2,5,8,11,14,2,5,8,11,14,2,5,8,11,14,3,13,5,8,11])
universes = np.full((17,17), fuel_univ)
universes[nonfuel_y, nonfuel_x] = univs['GT empty']
# Instantiate the lattice
lattice = openmc.RectLattice(name='Pin lattice')
lattice.lower_left = (-17.*pin_pitch/2., -17.*pin_pitch/2.)
lattice.pitch = (pin_pitch, pin_pitch)
lattice.universes = universes
# Add lattice to bounding cell
root_universe = openmc.Universe(name='Root universe')
cell = openmc.Cell(name='Lattice cell')
cell.fill = lattice
z_bounds = +surfs['bottom FR'] & -surfs['top FR']
cell.region = surfs['lat grid box inner'] & z_bounds
root_universe.add_cell(cell)
# Apply reflective boundaries on sides and vacuum on bottom/top
surfs['bottom FR'].boundary_type = 'vacuum'
surfs['top FR'].boundary_type = 'vacuum'
for halfspace in surfs['lat grid box inner']:
halfspace.surface.boundary_type = 'reflective'
# Define geometry with a single assembly
geometry = openmc.Geometry(root_universe)
h = active_fuel_length / args.axial
fuel_mats = {}
# Count the number of instances for each cell and material
if args.clone:
geometry.determine_paths(instances_only=True)
for cell in tqdm(geometry.get_all_material_cells().values(),
desc='Differentiating materials / assigning volume'):
if cell.fill in materials:
# Determine if this material is fuel
is_fuel = 'UO2 Fuel' in cell.fill.name
# Fill cell with list of "differentiated" materials if requested
if args.clone:
cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
# Determine volume of each fuel material
if is_fuel:
upper_right = cell.region.bounding_box[1]
if isclose(upper_right[0], rings[0]):
ri, ro = 0.0, rings[0]
elif isclose(upper_right[0], rings[1]):
ri, ro = rings[0], rings[1]
else:
ri, ro = rings[1], pellet_OR
if args.clone:
for mat in cell.fill:
mat.volume = pi * (ro*ro - ri*ri) * h
else:
# In non-clone mode, we still need to create a copy of the
# material for each ring since they get different volumes
if ri not in fuel_mats:
cell.fill = cell.fill.clone()
cell.fill.volume = pi * (ro*ro - ri*ri) * h
fuel_mats[ri] = cell.fill
else:
cell.fill = fuel_mats[ri]
else:
if args.clone:
for mat in cell.fill:
mat.volume = 1.0
else:
cell.fill.volume = 1.0
#### Create OpenMC "materials.xml" file
print('Getting materials...')
all_materials = geometry.get_all_materials()
print('Creating materials collection...')
materials = openmc.Materials(all_materials.values())
print('Exporting materials to XML...')
materials.export_to_xml(str(directory / 'materials.xml'))
#### Create OpenMC "geometry.xml" file
geometry.export_to_xml(str(directory / 'geometry.xml'))
#### Create OpenMC "settings.xml" file
# Construct uniform initial source distribution over fissionable zones
lower_left = (-lattice_pitch/2, -lattice_pitch/2, bottom_fuel_stack)
upper_right = (lattice_pitch/2, lattice_pitch/2, top_active_core)
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings = openmc.Settings()
settings.batches = 200
settings.inactive = 100
settings.particles = 10000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
settings.temperature = {
'default': 531.5,
'method': 'interpolation',
'range': (500.0, 1300.0)
}
if args.multipole:
settings.temperature['multipole'] = True
settings.temperature['tolerance'] = 1000
settings.export_to_xml(str(directory / 'settings.xml'))

192
smr/build-assembly-short.py
View file

@ -1,192 +0,0 @@
#!/usr/bin/env python3
import argparse
from math import pi, isclose
from pathlib import Path
import numpy as np
from tqdm import tqdm
import openmc
from smr.materials import materials, clone
from smr.surfaces import surfs, lattice_pitch, pin_pitch, bottom_fuel_stack, \
top_active_core, pellet_OR, clad_OR, clad_IR, guide_tube_IR, guide_tube_OR
from smr.pins import pin_universes
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='mat',
help='Whether to use distribmats or distribcells for tallies')
parser.add_argument('-a', '--axial', type=int, default=10,
help='Number of axial subdivisions in fuel')
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(multipole=True)
args = parser.parse_args()
# Make directory for inputs
if args.output_dir is None:
if args.depleted:
directory = Path('assembly-short-depleted')
else:
directory = Path('assembly-short-fresh')
else:
directory = args.output_dir
directory.mkdir(exist_ok=True)
rings = [0.1*pin_pitch, 0.2*pin_pitch]
# Define the NumPy array indices for assembly locations where there
# may be CR guide tubes, instrument tubes and burnable absorbers
nonfuel_y = np.array([2,2,2,3,3,5,5,5,5,5,8,8,8,8,8,11,11,11,11,11,13,13,14,14,14])
nonfuel_x = np.array([5,8,11,3,13,2,5,8,11,14,2,5,8,11,14,2,5,8,11,14,3,13,5,8,11])
# NO BURNABLE ABSORBERS
pins = pin_universes(rings, args.axial, args.depleted)
gtu = pins['GT empty']
#gti = pins['GT empty instr']
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel pin (3.1%) no grid']
universes[nonfuel_y, nonfuel_x] = [ gtu, gtu, gtu,
gtu, gtu,
gtu, gtu, gtu, gtu, gtu,
gtu, gtu, gtu, gtu, gtu,
gtu, gtu, gtu, gtu, gtu,
gtu, gtu,
gtu, gtu, gtu ]
# Instantiate the lattice
lattice = openmc.RectLattice(name='Pin lattice')
lattice.lower_left = (-17.*pin_pitch/2., -17.*pin_pitch/2.)
lattice.pitch = (pin_pitch, pin_pitch)
lattice.universes = universes
# Add lattice to bounding cell
root_universe = openmc.Universe(name='Root universe')
cell = openmc.Cell(name='Lattice cell')
cell.fill = lattice
z_bounds = +surfs['bot active core'] & -surfs['top active core']
cell.region = surfs['lat grid box inner'] & z_bounds
root_universe.add_cell(cell)
# Apply reflective boundaries
surfs['bot active core'].boundary_type = 'reflective'
surfs['top active core'].boundary_type = 'reflective'
for halfspace in surfs['lat grid box inner']:
halfspace.surface.boundary_type = 'reflective'
# Define geometry with a single assembly
geometry = openmc.Geometry(root_universe)
#### "Differentiate" the geometry if using distribmats
h = 10.0*pin_pitch / args.axial
if args.tallies == 'mat':
# Count the number of instances for each cell and material
geometry.determine_paths(instances_only=True)
for cell in tqdm(geometry.get_all_material_cells().values(),
desc='Differentiating materials'):
if cell.fill in materials:
# Fill cell with list of "differentiated" materials
cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
# Determine volume of each fuel material
if 'UO2 Fuel' in cell.fill[0].name:
lower_left, _ = cell.region.bounding_box
if isclose(lower_left[0], rings[0]):
ri, ro = 0.0, rings[0]
elif isclose(lower_left[0], rings[1]):
ri, ro = rings[0], rings[1]
else:
ri, ro = rings[1], pellet_OR
for mat in cell.fill:
mat.volume = pi * (ro*ro - ri*ri) * h
elif cell.fill[0].name == 'Borated Water':
for mat in cell.fill:
mat.volume = pin_pitch**2 - pi*clad_OR**2 * h
elif cell.fill[0].name == 'Helium':
for mat in cell.fill:
mat.volume = pi * (clad_IR**2 - pellet_OR**2) * h
elif cell.fill[0].name == 'M5':
for mat in cell.fill:
mat.volume = pi * (clad_OR**2 - clad_IR**2) * h
elif cell.fill[0].name == 'Zircaloy-4':
for mat in cell.fill:
mat.volume = pi * (guide_tube_OR**2 - guide_tube_IR**2) * h
#### Create OpenMC "materials.xml" file
print('Getting materials...')
all_materials = geometry.get_all_materials()
print('Creating materials collection...')
materials = openmc.Materials(all_materials.values())
print('Exporting materials to XML...')
materials.export_to_xml(str(directory / 'materials.xml'))
#### Create OpenMC "geometry.xml" file
geometry.export_to_xml(str(directory / 'geometry.xml'))
#### Create OpenMC "settings.xml" file
# Construct uniform initial source distribution over fissionable zones
lower_left = (-lattice_pitch/2, -lattice_pitch/2, bottom_fuel_stack)
upper_right = (lattice_pitch/2, lattice_pitch/2, top_active_core)
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings = openmc.Settings()
settings.batches = 200
settings.inactive = 100
settings.particles = 10000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
if args.multipole:
settings.temperature = {
'multipole': True,
'tolerance': 1000,
'default': 531.5,
'method': 'interpolation',
'range': (500.0, 1300.0)
}
settings.export_to_xml(str(directory / 'settings.xml'))
#### Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
# Extract all fuel materials
materials = geometry.get_materials_by_name(name='Fuel', matching=False)
# If using distribcells, create distribcell tally needed for depletion
if args.tallies == 'cell':
# Extract all cells filled by a fuel material
fuel_cells = []
for cell in geometry.get_all_cells().values():
if cell.fill in materials:
tally = openmc.Tally(name='depletion tally')
tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
'fission', '(n,2n)', '(n,3n)', '(n,4n)']
tally.nuclides = cell.fill.get_nuclides()
tally.filters.append(openmc.DistribcellFilter([cell]))
tallies.append(tally)
# If using distribmats, create material tally needed for depletion
elif args.tallies == 'mat':
tally = openmc.Tally(name='depletion tally')
tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
'fission', '(n,2n)', '(n,3n)', '(n,4n)']
tally.nuclides = materials[0].get_nuclides()
tally.filters = [openmc.MaterialFilter(materials)]
tallies.append(tally)
tallies.export_to_xml(str(directory / 'tallies.xml'))

55
smr/build-assembly.py
View file

@ -1,28 +1,23 @@
#!/usr/bin/env python3
import argparse
import copy
from pathlib import Path
import numpy as np
from tqdm import tqdm
import openmc
from smr.materials import materials, clone
from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core, pellet_OR
from smr.materials import materials
from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core
from smr.assemblies import assembly_universes
from smr import inlet_temperature
from smr.plots import assembly_plots
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', dest='multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('--clone', action='store_true',
help='Clone materials for each cell instance')
parser.add_argument('--no-clone', dest='clone', action='store_false',
help='Do not clone materials for each cell instance')
parser.add_argument('-m', '--multipole', action='store_true',
help='Whether to use multipole cross sections')
parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='mat',
help='Whether to use distribmats or distribcells for tallies')
parser.add_argument('-r', '--rings', type=int, default=10,
@ -32,7 +27,6 @@ parser.add_argument('-a', '--axial', type=int, default=196,
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(clone=False, multipole=True)
args = parser.parse_args()
# Make directory for inputs
@ -46,22 +40,26 @@ else:
directory.mkdir(exist_ok=True)
# Define geometry with a single assembly
if args.rings > 1:
ring_radii = np.sqrt(np.arange(1, args.rings)*pellet_OR**2 / args.rings)
else:
ring_radii = None
assembly = assembly_universes(ring_radii, args.axial, args.depleted)
lattice_sides = openmc.model.rectangular_prism(lattice_pitch, lattice_pitch,
boundary_type='reflective')
assembly = assembly_universes(args.rings, args.axial, args.depleted)
lattice_sides = openmc.model.get_rectangular_prism(lattice_pitch, lattice_pitch,
boundary_type='reflective')
main_cell = openmc.Cell(
fill=assembly['Assembly (3.1%)'],
fill=assembly['Assembly (3.1%) 16BA'],
region=lattice_sides & +surfs['lower bound'] & -surfs['upper bound']
)
root_univ = openmc.Universe(cells=[main_cell])
geometry = openmc.Geometry(root_univ)
def clone(material):
"""Perform copy of material but share nuclide densities"""
shared_mat = copy.copy(material)
shared_mat.id = None
return shared_mat
#### "Differentiate" the geometry if using distribmats
if args.clone:
if args.tallies == 'mat':
# Count the number of instances for each cell and material
geometry.determine_paths(instances_only=True)
@ -101,15 +99,10 @@ settings.inactive = 100
settings.particles = 10000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
settings.temperature = {
'default': inlet_temperature,
'method': 'interpolation',
'range': (300.0, 1500.0),
}
settings.sourcepoint_write = False
if args.multipole:
settings.temperature['multipole'] = True
settings.temperature['tolerance'] = 1000
settings.temperature = {'multipole': True, 'tolerance': 1000}
settings.export_to_xml(str(directory / 'settings.xml'))
@ -143,3 +136,7 @@ elif args.tallies == 'mat':
tallies.append(tally)
tallies.export_to_xml(str(directory / 'tallies.xml'))
# Create plots
plots = assembly_plots(main_cell.fill)
plots.export_to_xml(str(directory / 'plots.xml'))

125
smr/build-core-fresh.py
View file

@ -1,30 +1,33 @@
#!/usr/bin/env python3
import os
import shutil
import copy
import argparse
from math import pi
from pathlib import Path
import numpy as np
import openmc
from tqdm import tqdm
from smr.materials import materials, clone
from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core, \
pellet_OR, active_fuel_length
import openmc
from smr.materials import materials
from smr.plots import core_plots
from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core
from smr.core import core_geometry
from smr import inlet_temperature
def clone(mat):
"""Make a shallow copy of a material (sharing compositions)."""
mat_copy = copy.copy(mat)
mat_copy.id = None
return mat_copy
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', dest='multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('--clone', action='store_true',
help='Clone materials for each cell instance')
parser.add_argument('--no-clone', dest='clone', action='store_false',
help='Do not clone materials for each cell instance')
parser.add_argument('-m', '--multipole', action='store_true',
help='Whether to use multipole cross sections')
parser.add_argument('-t', '--tallies', choices=('cell', 'mat'), default='cell',
help='Whether to use distribmats or distribcells for tallies')
parser.add_argument('-r', '--rings', type=int, default=10,
help='Number of annular regions in fuel')
parser.add_argument('-a', '--axial', type=int, default=196,
@ -32,7 +35,6 @@ parser.add_argument('-a', '--axial', type=int, default=196,
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(clone=False, multipole=True)
args = parser.parse_args()
# Make directory for inputs
@ -45,44 +47,20 @@ else:
directory = args.output_dir
directory.mkdir(exist_ok=True)
if args.rings > 1:
ring_radii = np.sqrt(np.arange(1, args.rings)*pellet_OR**2 / args.rings)
else:
ring_radii = None
geometry = core_geometry(ring_radii, args.axial, args.depleted)
geometry = core_geometry(args.rings, args.axial, args.depleted)
h = active_fuel_length / args.axial
fuel_mats = {}
# Count the number of instances for each cell and material
if args.clone:
#### "Differentiate" the geometry if using distribmats
if args.tallies == 'mat':
# Count the number of instances for each cell and material
geometry.determine_paths(instances_only=True)
fuel_volume = pi * pellet_OR**2 * h / args.rings
for cell in tqdm(geometry.get_all_cells().values(),
desc='Differentiating materials / assigning volume'):
if cell.fill in materials:
# Determine if this material is fuel
name = cell.fill.name
is_fuel = 'UO2 Fuel' in name
# Extract all cells filled by a fuel material
fuel_mats = {m for m in materials if 'UO2 Fuel' in m.name}
# Determine volume of each fuel material
if is_fuel:
if args.clone:
# Fill cell with list of "differentiated" materials if requested
cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
for mat in cell.fill:
mat.volume = fuel_volume
else:
r_o = cell.region.bounding_box[1][0]
if (name, r_o) not in fuel_mats:
cell.fill = cell.fill.clone()
cell.fill.volume = fuel_volume
fuel_mats[name, r_o] = cell.fill
else:
cell.fill = fuel_mats[name, r_o]
else:
cell.fill.volume = 1.0
for cell in geometry.get_all_cells().values():
if cell.fill in fuel_mats:
# Fill cell with list of "differentiated" materials
cell.fill = [clone(cell.fill) for i in range(cell.num_instances)]
#### Create OpenMC "materials.xml" file
all_materials = geometry.get_all_materials()
@ -108,14 +86,45 @@ settings.inactive = 100
settings.particles = 10000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
settings.temperature = {
'default': inlet_temperature,
'method': 'interpolation',
'range': (300.0, 1500.0),
}
settings.sourcepoint_write = False
if args.multipole:
settings.temperature['multipole'] = True
settings.temperature['tolerance'] = 1000
settings.temperature = {'multipole': True, 'tolerance': 1000}
settings.export_to_xml(str(directory / 'settings.xml'))
#### Create OpenMC "plots.xml" file
plots = core_plots()
plots.export_to_xml(str(directory / 'plots.xml'))
#### Create OpenMC "tallies.xml" file
tallies = openmc.Tallies()
# Extract all fuel materials
materials = geometry.get_materials_by_name(name='Fuel', matching=False)
# If using distribcells, create distribcell tally needed for depletion
if args.tallies == 'cell':
# Extract all cells filled by a fuel material
fuel_cells = []
for cell in geometry.get_all_cells().values():
if cell.fill in materials:
tally = openmc.Tally(name='depletion tally')
tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
'fission', '(n,2n)', '(n,3n)', '(n,4n)']
tally.nuclides = cell.fill.get_nuclides()
tally.filters.append(openmc.DistribcellFilter([cell]))
tallies.append(tally)
# If using distribmats, create material tally needed for depletion
elif args.tallies == 'mat':
tally = openmc.Tally(name='depletion tally')
tally.scores = ['(n,p)', '(n,a)', '(n,gamma)',
'fission', '(n,2n)', '(n,3n)', '(n,4n)']
tally.nuclides = materials[0].get_nuclides()
tally.filters = [openmc.MaterialFilter(materials)]
tallies.append(tally)
tallies.export_to_xml(str(directory / 'tallies.xml'))

118
smr/build-core-long.py
View file

@ -1,118 +0,0 @@
#!/usr/bin/env python3
import argparse
from math import pi, isclose
from pathlib import Path
import openmc
from smr.materials import materials
from smr.surfaces import bottom_fuel_stack, top_active_core, \
pellet_OR, pin_pitch, clad_IR, clad_OR, active_fuel_length
from smr.core import core_geometry
from smr import inlet_temperature
import smr.surfaces
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', dest='multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('-a', '--axial', type=int, default=100,
help='Number of axial subdivisions in fuel')
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(multipole=True)
args = parser.parse_args()
# Make directory for inputs
if args.output_dir is None:
if args.depleted:
directory = Path('core-long-depleted')
else:
directory = Path('core-long-fresh')
else:
directory = args.output_dir
directory.mkdir(exist_ok=True)
# Modify lattice pitch
smr.surfaces.lattice_pitch = lattice_pitch = 17*smr.surfaces.pin_pitch
ring_radii = [0.1*pin_pitch, 0.2*pin_pitch]
geometry = core_geometry(ring_radii, args.axial, args.depleted)
h = active_fuel_length / args.axial
fuel_mats = {}
for cell in geometry.get_all_cells().values():
if cell.fill in materials:
# Determine volume of each fuel material
name = cell.fill.name
if 'UO2 Fuel' in name:
upper_right = cell.region.bounding_box[1][0]
if isclose(upper_right, ring_radii[0]):
ri, ro = 0.0, ring_radii[0]
elif isclose(upper_right, ring_radii[1]):
ri, ro = ring_radii[0], ring_radii[1]
else:
ri, ro = ring_radii[1], pellet_OR
if (name, ri) not in fuel_mats:
cell.fill = cell.fill.clone()
cell.fill.volume = pi * (ro*ro - ri*ri) * h
fuel_mats[name, ri] = cell.fill
else:
cell.fill = fuel_mats[name, ri]
elif name == 'Helium':
cell.fill.volume = pi * (clad_IR**2 - pellet_OR**2) * h
elif name == 'M5':
# Clad is not subdivided
cell.fill.volume = pi * (clad_OR**2 - clad_IR**2) * active_fuel_length
else:
cell.fill.volume = 1.0
#### Create OpenMC "materials.xml" file
all_materials = geometry.get_all_materials()
materials = openmc.Materials(all_materials.values())
materials.export_to_xml(str(directory / 'materials.xml'))
#### Create OpenMC "geometry.xml" file
geometry.export_to_xml(str(directory / 'geometry.xml'))
#### Create OpenMC "settings.xml" file
# Construct uniform initial source distribution over fissionable zones
lower_left = [-7.*lattice_pitch/2., -7.*lattice_pitch/2., bottom_fuel_stack]
upper_right = [+7.*lattice_pitch/2., +7.*lattice_pitch/2., top_active_core]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings = openmc.Settings()
settings.batches = 200
settings.inactive = 100
settings.particles = 20_000_000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
settings.temperature = {
'default': inlet_temperature,
'method': 'interpolation',
'range': (300.0, 1500.0),
}
if args.multipole:
settings.temperature['multipole'] = True
settings.temperature['tolerance'] = 1000
settings.export_to_xml(str(directory / 'settings.xml'))
# Check assembly power distribution
core_lattice = geometry.get_cells_by_fill_name('Main core')[0].fill
mesh = openmc.RegularMesh.from_rect_lattice(core_lattice)
assembly_power = openmc.Tally()
assembly_power.filters = [openmc.MeshFilter(mesh)]
assembly_power.scores = ['nu-fission']
tallies = openmc.Tallies([assembly_power])
tallies.export_to_xml(directory / 'tallies.xml')

131
smr/build-core-short.py
View file

@ -1,131 +0,0 @@
#!/usr/bin/env python3
import argparse
from math import pi, isclose
from pathlib import Path
import openmc
from smr.materials import materials
from smr.surfaces import bottom_fuel_stack, top_active_core, \
pellet_OR, surfs, pin_pitch, clad_IR, clad_OR
import smr.surfaces
import smr.pins
from smr.core import core_geometry
from smr import inlet_temperature
# Define command-line options
parser = argparse.ArgumentParser()
parser.add_argument('--multipole', action='store_true',
help='Use multipole cross sections')
parser.add_argument('--no-multipole', action='store_false',
help='Do not use multipole cross sections')
parser.add_argument('-a', '--axial', type=int, default=3,
help='Number of axial subdivisions in fuel')
parser.add_argument('-d', '--depleted', action='store_true',
help='Whether UO2 compositions should represent depleted fuel')
parser.add_argument('-o', '--output-dir', type=Path, default=None)
parser.set_defaults(multipole=True)
args = parser.parse_args()
# Make directory for inputs
if args.output_dir is None:
if args.depleted:
directory = Path('core-short-depleted')
else:
directory = Path('core-short-fresh')
else:
directory = args.output_dir
directory.mkdir(exist_ok=True)
# Modify lattice pitch
smr.surfaces.lattice_pitch = lattice_pitch = 17*smr.surfaces.pin_pitch
# Modify fuel length
length = 3. * pin_pitch
smr.surfaces.active_fuel_length = length
smr.pins.top_active_core = length
surfs['top active core'].z0 = length
# Change top and bottom of model to contain only fuel
surfs['lower bound'].z0 = 0.0
surfs['lower bound'].boundary_type = 'reflective'
surfs['upper bound'].z0 = length
surfs['upper bound'].boundary_type = 'reflective'
ring_radii = [0.1*pin_pitch, 0.2*pin_pitch]
geometry = core_geometry(ring_radii, args.axial, args.depleted)
h = length / args.axial
fuel_mats = {}
for cell in geometry.get_all_cells().values():
if cell.fill in materials:
# Determine volume of each fuel material
name = cell.fill.name
if 'UO2 Fuel' in name:
upper_right = cell.region.bounding_box[1][0]
if isclose(upper_right, ring_radii[0]):
ri, ro = 0.0, ring_radii[0]
elif isclose(upper_right, ring_radii[1]):
ri, ro = ring_radii[0], ring_radii[1]
else:
ri, ro = ring_radii[1], pellet_OR
if (name, ri) not in fuel_mats:
cell.fill = cell.fill.clone()
cell.fill.volume = pi * (ro*ro - ri*ri) * h
fuel_mats[name, ri] = cell.fill
else:
cell.fill = fuel_mats[name, ri]
elif name == 'Helium':
cell.fill.volume = pi * (clad_IR**2 - pellet_OR**2) * h
elif name == 'M5':
# Clad is not subdivided
cell.fill.volume = pi * (clad_OR**2 - clad_IR**2) * length
else:
cell.fill.volume = 1.0
#### Create OpenMC "materials.xml" file
all_materials = geometry.get_all_materials()
materials = openmc.Materials(all_materials.values())
materials.export_to_xml(str(directory / 'materials.xml'))
#### Create OpenMC "geometry.xml" file
geometry.export_to_xml(str(directory / 'geometry.xml'))
#### Create OpenMC "settings.xml" file
# Construct uniform initial source distribution over fissionable zones
lower_left = [-7.*lattice_pitch/2., -7.*lattice_pitch/2., bottom_fuel_stack]
upper_right = [+7.*lattice_pitch/2., +7.*lattice_pitch/2., top_active_core]
source = openmc.source.Source(space=openmc.stats.Box(lower_left, upper_right))
source.space.only_fissionable = True
settings = openmc.Settings()
settings.batches = 200
settings.inactive = 100
settings.particles = 10000
settings.output = {'tallies': False, 'summary': False}
settings.source = source
settings.sourcepoint = {'write': False}
settings.temperature = {
'default': inlet_temperature,
'method': 'interpolation',
'range': (300.0, 1500.0),
}
if args.multipole:
settings.temperature['multipole'] = True
settings.temperature['tolerance'] = 1000
settings.export_to_xml(str(directory / 'settings.xml'))
# Check assembly power distribution
core_lattice = geometry.get_cells_by_fill_name('Main core')[0].fill
mesh = openmc.RegularMesh.from_rect_lattice(core_lattice)
assembly_power = openmc.Tally()
assembly_power.filters = [openmc.MeshFilter(mesh)]
assembly_power.scores = ['nu-fission']
tallies = openmc.Tallies([assembly_power])
tallies.export_to_xml(directory / 'tallies.xml')

20
smr/milestones.md
View file

@ -1,20 +0,0 @@
# Milestone Models
- **AD-SE-08-61, Coupled Multiphysics Driver Implementation** --- This milestone
used the singlerod short/long problems. Generating the model was done with the
script `tests/singlerod/make_openmc_model.py` from the ENRICO repository
(there is a `--short` command line option to generate the short version)
- **AD-SE-08-66, Coupled Assembly Analysis** --- This milestone used the
assembly short and long (v2) problems. Generating the model was done with
`smr/build-assembly-long.py -a 100 --clone` on the ecp-benchmarks repository
(git commit `631fefe`, after pull request #14). In these models, materials are
fully differentiated across each fuel ring/axial segment.
- **AD-SE-08-73, Full core coupled-physics simulation** --- This milestone used
the core-short and core-long (90 layer) models. Generating the models was done
with `smr/build-core-short.py` and `smr/build-assembly-long.py -a 90` on the
ecp-benchmarks repository (git commit `c7b89db`, after pull request #16). In
these models, materials are not differentiated and no grid spacers are
present. The lattice pitch is modified to be exactly 17 times the pin pitch
(slightly different than NuScale specification).

6
smr/smr/__init__.py
View file

@ -1,6 +0,0 @@
_MPA_PER_PSIA = 0.0068947572931683625
# Values from ML17013A274, Table 4.1-1
inlet_temperature = (497 - 32)*5/9 + 273.15
core_average_temperature = (543 - 32)*5/9 + 273.15
system_pressure = 1850*_MPA_PER_PSIA

170
smr/smr/assemblies.py
View file

@ -120,14 +120,13 @@ def make_assembly(name, universes):
return universe
def assembly_universes(ring_radii, num_axial, depleted):
def assembly_universes(num_rings, num_axial, depleted):
"""Generate universes for SMR fuel assemblies.
Parameters
----------
ring_radii : iterable of float
Radii of rings in fuel (note that this doesn't need to include the
full fuel pin radius)
num_rings : int
Number of annual regions in fuel
num_axial : int
Number of axial subdivisions in fuel
depleted : bool
@ -139,14 +138,15 @@ def assembly_universes(ring_radii, num_axial, depleted):
Dictionary mapping a universe name to a openmc.Universe object
"""
pins = pin_universes(ring_radii, num_axial, depleted)
pins = pin_universes(num_rings, num_axial, depleted)
# Create dictionary to store assembly universes
univs = {}
# commonly needed universes
gtu = pins['GT empty stack']
gtu = pins['GT empty']
gti = pins['GT empty instr']
bas = pins['BA stack']
ins = pins['IT stack']
crA = pins['GT CR bank A']
crB = pins['GT CR bank B']
@ -215,7 +215,7 @@ def assembly_universes(ring_radii, num_axial, depleted):
gtu, gtu, gtu, gtu, gtu,
gtu, gtu,
gtu, gtu, gtu ]
univs['Assembly (2.4%)' + comment] = \
univs['Assembly (2.4%) no BAs' + comment] = \
make_assembly('Assembly (2.4%) no BAs' + comment, universes)
# WITH CONTROL ROD D BANK
@ -231,6 +231,33 @@ def assembly_universes(ring_radii, num_axial, depleted):
univs['Assembly (2.4%) CR D' + comment] = \
make_assembly('Assembly (2.4%) CR D' + comment, universes)
# WITH 12 BURNABLE ABSORBERS
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (2.4%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, gtu, bas,
bas, bas,
bas, gtu, gtu, gtu, bas,
gtu, gtu, cent, gtu, gtu,
bas, gtu, gtu, gtu, bas,
bas, bas,
bas, gtu, bas ]
univs['Assembly (2.4%) 12BA' + comment] = \
make_assembly('Assembly (2.4%) 12BA' + comment, universes)
# WITH 16 BURNABLE ABSORBERS
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (2.4%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, bas, bas,
bas, bas,
bas, gtu, gtu, gtu, bas,
bas, gtu, cent, gtu, bas,
bas, gtu, gtu, gtu, bas,
bas, bas,
bas, bas, bas ]
univs['Assembly (2.4%) 16BA' + comment] = \
make_assembly('Assembly (2.4%) 16BA' + comment, universes)
#### 3.1% ENRICHED ASSEMBLIES
for cent, comment in [(gti, ''), (ins, ' instr')]:
@ -261,5 +288,134 @@ def assembly_universes(ring_radii, num_axial, depleted):
univs['Assembly (3.1%) CR SA' + comment] = \
make_assembly('Assembly (3.1%) CR SA' + comment, universes)
# WITH 20 BURNABLE ABSORBERS
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, bas, bas,
bas, bas,
bas, bas, gtu, bas, bas,
bas, gtu, cent, gtu, bas,
bas, bas, gtu, bas, bas,
bas, bas,
bas, bas, bas ]
univs['Assembly (3.1%) 20BA' + comment] = \
make_assembly('Assembly (3.1%) 20BA' + comment, universes)
# WITH 16 BURNABLE ABSORBERS
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, bas, bas,
bas, bas,
bas, gtu, gtu, gtu, bas,
bas, gtu, cent, gtu, bas,
bas, gtu, gtu, gtu, bas,
bas, bas,
bas, bas, bas ]
univs['Assembly (3.1%) 16BA' + comment] = \
make_assembly('Assembly (3.1%) 16BA' + comment, universes)
# WITH 15 BURNABLE ABSORBERS NW
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ gtu, gtu, gtu,
gtu, gtu,
gtu, bas, bas, bas, bas,
gtu, bas, cent, bas, bas,
gtu, bas, bas, bas, bas,
gtu, bas,
bas, bas, bas ]
univs['Assembly (3.1%) 15BANW' + comment] = \
make_assembly('Assembly (3.1%) 15BANW' + comment, universes)
# WITH 15 BURNABLE ABSORBERS NE
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ gtu, gtu, gtu,
gtu, gtu,
bas, bas, bas, bas, gtu,
bas, bas, cent, bas, gtu,
bas, bas, bas, bas, gtu,
bas, gtu,
bas, bas, bas ]
univs['Assembly (3.1%) 15BANE' + comment] = \
make_assembly('Assembly (3.1%) 15BANE' + comment, universes)
# WITH 15 BURNABLE ABSORBERS SW
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, bas, bas,
gtu, bas,
gtu, bas, bas, bas, bas,
gtu, bas, cent, bas, bas,
gtu, bas, bas, bas, bas,
gtu, gtu,
gtu, gtu, gtu ]
univs['Assembly (3.1%) 15BASW' + comment] = \
make_assembly('Assembly (3.1%) 15BASW' + comment, universes)
# WITH 15 BURNABLE ABSORBERS SE
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, bas, bas,
bas, gtu,
bas, bas, bas, bas, gtu,
bas, bas, cent, bas, gtu,
bas, bas, bas, bas, gtu,
gtu, gtu,
gtu, gtu, gtu ]
univs['Assembly (3.1%) 15BASE' + comment] = \
make_assembly('Assembly (3.1%) 15BASE' + comment, universes)
# WITH 6 BURNABLE ABSORBERS N
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ gtu, gtu, gtu,
gtu, gtu,
gtu, gtu, gtu, gtu, gtu,
gtu, gtu, cent, gtu, gtu,
bas, gtu, gtu, gtu, bas,
bas, bas,
bas, gtu, bas ]
univs['Assembly (3.1%) 6BAN' + comment] = \
make_assembly('Assembly (3.1%) 6BAN' + comment, universes)
# WITH 6 BURNABLE ABSORBERS S
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, gtu, bas,
bas, bas,
bas, gtu, gtu, gtu, bas,
gtu, gtu, cent, gtu, gtu,
gtu, gtu, gtu, gtu, gtu,
gtu, gtu,
gtu, gtu, gtu ]
univs['Assembly (3.1%) 6BAS' + comment] = \
make_assembly('Assembly (3.1%) 6BAS' + comment, universes)
# WITH 6 BURNABLE ABSORBERS W
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ gtu, gtu, bas,
gtu, bas,
gtu, gtu, gtu, gtu, bas,
gtu, gtu, cent, gtu, gtu,
gtu, gtu, gtu, gtu, bas,
gtu, bas,
gtu, gtu, bas ]
univs['Assembly (3.1%) 6BAW' + comment] = \
make_assembly('Assembly (3.1%) 6BAW' + comment, universes)
# WITH 6 BURNABLE ABSORBERS E
universes = np.empty((17,17), dtype=openmc.Universe)
universes[:,:] = pins['Fuel (3.1%) stack']
universes[nonfuel_y, nonfuel_x] = [ bas, gtu, gtu,
bas, gtu,
bas, gtu, gtu, gtu, gtu,
gtu, gtu, cent, gtu, gtu,
bas, gtu, gtu, gtu, gtu,
bas, gtu,
bas, gtu, gtu ]
univs['Assembly (3.1%) 6BAE' + comment] = \
make_assembly('Assembly (3.1%) 6BAE' + comment, universes)
return univs

154
smr/smr/core.py
View file

@ -5,19 +5,18 @@ import numpy as np
import openmc
from .materials import mats
from .surfaces import surfs, lattice_pitch
from .reflector import reflector_universes
from .assemblies import assembly_universes
from smr import surfaces
def core_geometry(ring_radii, num_axial, depleted):
def core_geometry(num_rings, num_axial, depleted):
"""Generate full core SMR geometry.
Parameters
----------
ring_radii : iterable of float
Radii of rings in fuel (note that this doesn't need to include the
full fuel pin radius)
num_rings : int
Number of annual regions in fuel
num_axial : int
Number of axial subdivisions in fuel
depleted : bool
@ -29,12 +28,11 @@ def core_geometry(ring_radii, num_axial, depleted):
SMR full core geometry
"""
assembly = assembly_universes(ring_radii, num_axial, depleted)
assembly = assembly_universes(num_rings, num_axial, depleted)
reflector = reflector_universes()
# Construct main core lattice
core = openmc.RectLattice(name='Main core')
lattice_pitch = surfaces.lattice_pitch
core.lower_left = (-9*lattice_pitch/2, -9*lattice_pitch/2)
core.pitch = (lattice_pitch, lattice_pitch)
universes = np.tile(reflector['solid'], (9, 9))
@ -47,67 +45,67 @@ def core_geometry(ring_radii, num_axial, depleted):
universes[1, 1] = reflector['1,1']
universes[1, 2] = reflector['NW']
universes[1, 3] = assembly['Assembly (3.1%)']
universes[1, 4] = assembly['Assembly (3.1%)']
universes[1, 5] = assembly['Assembly (3.1%)']
universes[1, 3] = assembly['Assembly (3.1%) instr']
universes[1, 4] = assembly['Assembly (2.4%) CR D']
universes[1, 5] = assembly['Assembly (3.1%) instr']
universes[1, 6] = reflector['NE']
universes[1, 7] = reflector['1,7']
universes[2, 0] = reflector['2,0']
universes[2, 1] = reflector['NW']
universes[2, 2] = assembly['Assembly (3.1%)']
universes[2, 3] = assembly['Assembly (2.4%)']
universes[2, 4] = assembly['Assembly (1.6%)']
universes[2, 5] = assembly['Assembly (2.4%)']
universes[2, 6] = assembly['Assembly (3.1%)']
universes[2, 2] = assembly['Assembly (3.1%) instr']
universes[2, 3] = assembly['Assembly (2.4%) CR D']
universes[2, 4] = assembly['Assembly (3.1%) 16BA']
universes[2, 5] = assembly['Assembly (2.4%) CR D']
universes[2, 6] = assembly['Assembly (3.1%) instr']
universes[2, 7] = reflector['NE']
universes[2, 8] = reflector['2,8']
universes[3, 0] = reflector['3,0']
universes[3, 1] = assembly['Assembly (3.1%)']
universes[3, 2] = assembly['Assembly (2.4%)']
universes[3, 3] = assembly['Assembly (1.6%)']
universes[3, 4] = assembly['Assembly (1.6%)']
universes[3, 5] = assembly['Assembly (1.6%)']
universes[3, 6] = assembly['Assembly (2.4%)']
universes[3, 7] = assembly['Assembly (3.1%)']
universes[3, 1] = assembly['Assembly (3.1%) instr']
universes[3, 2] = assembly['Assembly (2.4%) CR D']
universes[3, 3] = assembly['Assembly (3.1%) 16BA']
universes[3, 4] = assembly['Assembly (2.4%) CR D']
universes[3, 5] = assembly['Assembly (3.1%) 16BA']
universes[3, 6] = assembly['Assembly (2.4%) CR D']
universes[3, 7] = assembly['Assembly (3.1%) instr']
universes[3, 8] = reflector['3,8']
universes[4, 0] = reflector['4,0']
universes[4, 1] = assembly['Assembly (3.1%)']
universes[4, 2] = assembly['Assembly (1.6%)']
universes[4, 3] = assembly['Assembly (1.6%)']
universes[4, 4] = assembly['Assembly (2.4%)']
universes[4, 5] = assembly['Assembly (1.6%)']
universes[4, 6] = assembly['Assembly (1.6%)']
universes[4, 7] = assembly['Assembly (3.1%)']
universes[4, 1] = assembly['Assembly (2.4%) CR D']
universes[4, 2] = assembly['Assembly (3.1%) 16BA']
universes[4, 3] = assembly['Assembly (2.4%) CR D']
universes[4, 4] = assembly['Assembly (1.6%) instr']
universes[4, 5] = assembly['Assembly (2.4%) CR D']
universes[4, 6] = assembly['Assembly (3.1%) 16BA']
universes[4, 7] = assembly['Assembly (2.4%) CR D']
universes[4, 8] = reflector['4,8']
universes[5, 0] = reflector['5,0']
universes[5, 1] = assembly['Assembly (3.1%)']
universes[5, 2] = assembly['Assembly (2.4%)']
universes[5, 3] = assembly['Assembly (1.6%)']
universes[5, 4] = assembly['Assembly (1.6%)']
universes[5, 5] = assembly['Assembly (1.6%)']
universes[5, 6] = assembly['Assembly (2.4%)']
universes[5, 7] = assembly['Assembly (3.1%)']
universes[5, 1] = assembly['Assembly (3.1%) instr']
universes[5, 2] = assembly['Assembly (2.4%) CR D']
universes[5, 3] = assembly['Assembly (3.1%) 16BA']
universes[5, 4] = assembly['Assembly (2.4%) CR D']
universes[5, 5] = assembly['Assembly (3.1%) 16BA']
universes[5, 6] = assembly['Assembly (2.4%) CR D']
universes[5, 7] = assembly['Assembly (3.1%) instr']
universes[5, 8] = reflector['5,8']
universes[6, 0] = reflector['6,0']
universes[6, 1] = reflector['SW']
universes[6, 2] = assembly['Assembly (3.1%)']
universes[6, 3] = assembly['Assembly (2.4%)']
universes[6, 4] = assembly['Assembly (1.6%)']
universes[6, 5] = assembly['Assembly (2.4%)']
universes[6, 6] = assembly['Assembly (3.1%)']
universes[6, 2] = assembly['Assembly (3.1%) instr']
universes[6, 3] = assembly['Assembly (2.4%) CR D']
universes[6, 4] = assembly['Assembly (3.1%) 16BA']
universes[6, 5] = assembly['Assembly (2.4%) CR D']
universes[6, 6] = assembly['Assembly (3.1%) instr']
universes[6, 7] = reflector['SE']
universes[6, 8] = reflector['6,8']
universes[7, 1] = reflector['7,1']
universes[7, 2] = reflector['SW']
universes[7, 3] = assembly['Assembly (3.1%)']
universes[7, 4] = assembly['Assembly (3.1%)']
universes[7, 5] = assembly['Assembly (3.1%)']
universes[7, 3] = assembly['Assembly (3.1%) instr']
universes[7, 4] = assembly['Assembly (2.4%) CR D']
universes[7, 5] = assembly['Assembly (3.1%) instr']
universes[7, 6] = reflector['SE']
universes[7, 7] = reflector['7,7']
@ -120,7 +118,6 @@ def core_geometry(ring_radii, num_axial, depleted):
core.universes = universes
root_univ = openmc.Universe(universe_id=0, name='root universe')
surfs = surfaces.surfs
# Cylinder filled with core lattice
cell = openmc.Cell(name='Main core')
@ -136,10 +133,75 @@ def core_geometry(ring_radii, num_axial, depleted):
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
# Neutron shield panels
cell = openmc.Cell(name='neutron shield panel NW')
cell.fill = mats['SS']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
+surfs['neutron shield NWbot SEtop'] &
-surfs['neutron shield NWtop SEbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel N')
cell.fill = mats['H2O']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
+surfs['neutron shield NWtop SEbot'] &
-surfs['neutron shield NEtop SWbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel SE')
cell.fill = mats['SS']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-surfs['neutron shield NWbot SEtop'] &
+surfs['neutron shield NWtop SEbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel E')
cell.fill = mats['H2O']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
+surfs['neutron shield NWbot SEtop'] &
+surfs['neutron shield NEbot SWtop'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel NE')
cell.fill = mats['SS']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
+surfs['neutron shield NEbot SWtop'] &
-surfs['neutron shield NEtop SWbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel S')
cell.fill = mats['H2O']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-surfs['neutron shield NWtop SEbot'] &
+surfs['neutron shield NEtop SWbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel SW')
cell.fill = mats['SS']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-surfs['neutron shield NEbot SWtop'] &
+surfs['neutron shield NEtop SWbot'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
cell = openmc.Cell(name='neutron shield panel W')
cell.fill = mats['H2O']
cell.region = (+surfs['core barrel OR'] & -surfs['neutron shield OR'] &
-surfs['neutron shield NWbot SEtop'] &
-surfs['neutron shield NEbot SWtop'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)
# Downcomer
cell = openmc.Cell(name='downcomer')
cell.fill = mats['H2O']
cell.region = (+surfs['core barrel OR'] & -surfs['RPV IR'] &
cell.region = (+surfs['neutron shield OR'] & -surfs['RPV IR'] &
+surfs['lower bound'] & -surfs['upper bound'])
root_univ.add_cell(cell)

54
smr/smr/materials.py
View file

@ -1,11 +1,7 @@
"""Instantiate the OpenMC Materials needed by the core model."""
import copy
import openmc
from openmc.data import atomic_weight, atomic_mass, water_density
from . import system_pressure, core_average_temperature
from openmc.data import atomic_weight, atomic_mass
_DEPLETION_NUCLIDES = [
@ -45,11 +41,13 @@ mats = {}
# Create He gas material for fuel pin gap
mats['He'] = openmc.Material(name='Helium')
mats['He'].temperature = 300
mats['He'].set_density('g/cc', 0.0015981)
mats['He'].add_element('He', 1.0, 'ao')
# Create air material for instrument tubes
mats['Air'] = openmc.Material(name='Air')
mats['Air'].temperature = 300
mats['Air'].set_density('g/cc', 0.00616)
mats['Air'].add_element('O', 0.2095, 'ao')
mats['Air'].add_element('N', 0.7809, 'ao')
@ -58,6 +56,7 @@ mats['Air'].add_element('C', 0.00027, 'ao')
# Create inconel 718 material
mats['In'] = openmc.Material(name='Inconel')
mats['In'].temperature = 300
mats['In'].set_density('g/cc', 8.2)
mats['In'].add_element('Si', 0.0035, 'wo')
mats['In'].add_element('Cr', 0.1896, 'wo')
@ -65,17 +64,9 @@ mats['In'].add_element('Mn', 0.0087, 'wo')
mats['In'].add_element('Fe', 0.2863, 'wo')
mats['In'].add_element('Ni', 0.5119, 'wo')
# Create stainless steel 302
mats['SS302'] = openmc.Material(name='SS302')
mats['SS302'].set_density('g/cm3', 7.86)
mats['SS302'].add_element('Si', 0.01, 'wo')
mats['SS302'].add_element('Cr', 0.18, 'wo')
mats['SS302'].add_element('Mn', 0.02, 'wo')
mats['SS302'].add_element('Fe', 0.70, 'wo')
mats['SS302'].add_element('Ni', 0.09, 'wo')
# Create stainless steel material
mats['SS'] = openmc.Material(name='SS304')
mats['SS'].temperature = 300
mats['SS'].set_density('g/cc', 8.03)
mats['SS'].add_element('Si', 0.0060, 'wo')
mats['SS'].add_element('Cr', 0.1900, 'wo')
@ -85,6 +76,7 @@ mats['SS'].add_element('Ni', 0.1000, 'wo')
# Create carbon steel material
mats['CS'] = openmc.Material(name='Carbon Steel')
mats['CS'].temperature = 300
mats['CS'].set_density('g/cc', 7.8)
mats['CS'].add_element('C', 0.00270, 'wo')
mats['CS'].add_element('Mn', 0.00750, 'wo')
@ -105,6 +97,7 @@ mats['CS'].add_element('Fe', 0.96487, 'wo')
# Create zircaloy 4 material
mats['Zr'] = openmc.Material(name='Zircaloy-4')
mats['Zr'].temperature = 300
mats['Zr'].set_density('g/cc', 6.55)
mats['Zr'].add_element('O', 0.00125, 'wo')
mats['Zr'].add_element('Cr', 0.0010, 'wo')
@ -112,18 +105,9 @@ mats['Zr'].add_element('Fe', 0.0021, 'wo')
mats['Zr'].add_element('Zr', 0.98115, 'wo')
mats['Zr'].add_element('Sn', 0.0145, 'wo')
# Create M5 alloy material
m5_niobium = 0.01 # http://publications.jrc.ec.europa.eu/repository/bitstream/JRC100644/lcna28366enn.pdf
m5_oxygen = 0.00135 # http://publications.jrc.ec.europa.eu/repository/bitstream/JRC100644/lcna28366enn.pdf
m5_density = 6.494 # 10.1039/C5DT03403E
mats['M5'] = openmc.Material(name='M5')
mats['M5'].add_element('Zr', 1.0 - m5_niobium - m5_oxygen)
mats['M5'].add_element('Nb', m5_niobium)
mats['M5'].add_element('O', m5_oxygen)
mats['M5'].set_density('g/cm3', m5_density)
# Create Ag-In-Cd control rod material
mats['AIC'] = openmc.Material(name='Ag-In-Cd')
mats['AIC'].temperature = 300
mats['AIC'].set_density('g/cc', 10.16)
mats['AIC'].add_element('Ag', 0.80, 'wo')
mats['AIC'].add_element('In', 0.15, 'wo')
@ -132,11 +116,10 @@ mats['AIC'].add_element('Cd', 0.05, 'wo')
#### Borated Water
# Concentration of boron at beginning of equilibrium cycle
boron_ppm = 1240 # ML17013A274, Figure 4.3-17
boron_ppm = 975
# Density of water
h2o_dens = water_density(core_average_temperature, system_pressure)
# Density of clean water at 2250 psia T=560F NIST
h2o_dens = 0.73986
# Weight percent of natural boron in borated water
wB_Bh2o = boron_ppm * 1.0e-6
@ -163,6 +146,7 @@ aho_Bh2o = ah2o_Bh2o
# Create borated water for coolant / moderator
mats['H2O'] = openmc.Material(name='Borated Water')
mats['H2O'].temperature = 300
mats['H2O'].set_density('g/cc', rho_Bh2o)
mats['H2O'].add_element('B', aB_Bh2o, 'ao')
mats['H2O'].add_element('H', ah_Bh2o, 'ao')
@ -194,6 +178,7 @@ aB_bsg = aB10_bsg + aB11_bsg
# Create borosilicate glass material
mats['BSG'] = openmc.Material(name='Borosilicate Glass')
mats['BSG'].temperature = 300
mats['BSG'].set_density('g/cc', 2.26)
mats['BSG'].add_element('O', aO_bsg, 'ao')
mats['BSG'].add_element('Si', aSi_bsg, 'ao')
@ -206,6 +191,7 @@ mats['BSG'].add_nuclide('B11', aB11_bsg, 'ao')
# Create 1.6% enriched UO2 fuel material
mat = openmc.Material(name='1.6% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.31341)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=1.61006)
@ -213,6 +199,7 @@ mats['UO2 1.6 fresh'] = mat
# Create 2.4% enriched UO2 fuel material
mat = openmc.Material(name='2.4% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.29748)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=2.39993)
@ -220,6 +207,7 @@ mats['UO2 2.4 fresh'] = mat
# Create 3.1% enriched UO2 fuel material
mat = openmc.Material(name='3.1% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.30166)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=3.10221)
@ -227,6 +215,7 @@ mats['UO2 3.1 fresh'] = mat
# Depleted versions of 1.6%, 2.4%, 3.1% fuel
mat = openmc.Material(name='2.4% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.29748)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=2.39993)
@ -235,6 +224,7 @@ for nuc in _DEPLETION_NUCLIDES:
mats['UO2 2.4 depleted'] = mat
mat = openmc.Material(name='1.6% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.31341)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=1.61006)
@ -243,6 +233,7 @@ for nuc in _DEPLETION_NUCLIDES:
mats['UO2 1.6 depleted'] = mat
mat = openmc.Material(name='3.1% Enr. UO2 Fuel')
mat.temperature = 300
mat.set_density('g/cc', 10.30166)
mat.add_element('O', 2., 'ao')
mat.add_element('U', 1., 'ao', enrichment=3.10221)
@ -254,10 +245,3 @@ mats['UO2 3.1 depleted'] = mat
# Construct a collection of Materials to export to XML
materials = openmc.Materials(mats.values())
def clone(material):
"""Perform copy of material but share nuclide densities"""
shared_mat = copy.copy(material)
shared_mat.id = None
return shared_mat

295
smr/smr/pins.py
View file

@ -142,14 +142,13 @@ def make_pin_stack(name, zsurfaces, universes, boundary, fuel_fill):
return universe
def pin_universes(ring_radii=None, num_axial=196, depleted=False):
def pin_universes(num_rings=10, num_axial=196, depleted=False):
"""Generate universes for SMR fuel pins.
Parameters
----------
ring_radii : iterable of float
Radii of rings in fuel (note that this doesn't need to include the
full fuel pin radius)
num_rings : int
Number of annual regions in fuel
num_axial : int
Number of axial subdivisions in fuel
depleted : bool
@ -231,8 +230,8 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
surfs['top upper nozzle']
]
univs['GT empty stack'] = make_stack(
'GT empty stack', surfaces=stack_surfs,
univs['GT empty'] = make_stack(
'GT empty', surfaces=stack_surfs,
universes=[univs['water pin'],
univs['water pin'],
univs['water pin'],
@ -473,6 +472,199 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
univs['GT CR bank {} dummy bare'.format(b)]])
#### BURNABLE ABSORBER PIN CELLS
univs['BA'] = make_pin(
'BA',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['BA IR 7'],
surfs['BA IR 8']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']])
univs['BA grid (bottom)'] = make_pin(
'BA grid (bottom)',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['BA IR 7'],
surfs['BA IR 8']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']],
grid='bottom')
univs['BA grid (intermediate)'] = make_pin(
'BA grid (intermediate)',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['BA IR 7'],
surfs['BA IR 8']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']],
grid='intermediate')
univs['BA dashpot'] = make_pin(
'BA dashpot',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['GT dashpot IR'],
surfs['GT dashpot OR']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']])
univs['BA dashpot grid (bottom)'] = make_pin(
'BA dashpot grid (bottom)',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['GT dashpot IR'],
surfs['GT dashpot OR']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']],
grid='bottom')
univs['BA dashpot grid (intermediate)'] = make_pin(
'BA dashpot grid (intermediate)',
surfaces=[surfs['BA IR 1'],
surfs['BA IR 2'],
surfs['BA IR 3'],
surfs['BA IR 4'],
surfs['BA IR 5'],
surfs['BA IR 6'],
surfs['GT dashpot IR'],
surfs['GT dashpot OR']],
materials=[mats['Air'],
mats['SS'],
mats['Air'],
mats['BSG'],
mats['Air'],
mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']],
grid='intermediate')
univs['BA blank SS'] = make_pin(
'BA blank SS',
surfaces=[surfs['BA IR 6'],
surfs['BA IR 7'],
surfs['BA IR 8']],
materials=[mats['SS'],
mats['H2O'],
mats['Zr'],
mats['H2O']])
univs['BA blank SS bare'] = make_pin(
'BA blank SS bare',
surfaces=[surfs['BA IR 6']],
materials=[mats['SS'],
mats['H2O']])
stack_surfs_BA = [
surfs['bot support plate'],
surfs['top support plate'],
surfs['top lower nozzle'],
surfs['top lower thimble'],
surfs['grid1bot'],
surfs['BA bot'],
surfs['grid1top'],
surfs['dashpot top'],
surfs['grid2bot'],
surfs['grid2top'],
surfs['grid3bot'],
surfs['grid3top'],
surfs['grid4bot'],
surfs['grid4top'],
surfs['top active core'],
surfs['grid5bot'],
surfs['grid5top'],
surfs['top pin plenum'],
surfs['top FR'],
surfs['bot upper nozzle'],
surfs['top upper nozzle']]
# Stack all axial pieces of control rod tubes together for each bank
univs['BA stack'] = make_stack(
'BA stack', stack_surfs_BA,
universes=[univs['water pin'],
univs['water pin'],
univs['water pin'],
univs['GTd empty'],
univs['GTd empty'],
univs['GTd empty grid (bottom)'],
univs['BA dashpot grid (bottom)'],
univs['BA dashpot'],
univs['BA'],
univs['BA grid (intermediate)'],
univs['BA'],
univs['BA grid (intermediate)'],
univs['BA'],
univs['BA grid (intermediate)'],
univs['BA'],
univs['BA blank SS'],
univs['BA blank SS'],
univs['BA blank SS'],
univs['BA blank SS'],
univs['BA blank SS'],
univs['BA blank SS bare'],
univs['water pin']])
# Fuel pin cells
univs['SS pin'] = make_pin(
'SS pin',
@ -482,16 +674,16 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
univs['end plug'] = make_pin(
'end plug',
[surfs['clad OR']],
[mats['M5'], mats['H2O']])
[mats['Zr'], mats['H2O']])
univs['pin plenum'] = make_pin(
'pin plenum',
surfaces=[surfs['plenum spring OR'],
surfs['clad IR'],
surfs['clad OR']],
materials=[mats['SS302'],
materials=[mats['In'],
mats['He'],
mats['M5'],
mats['Zr'],
mats['H2O']])
univs['pin plenum grid (intermediate)'] = make_pin(
@ -508,49 +700,47 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
#### 1.6% ENRICHED FUEL PIN CELL
if num_axial > 1:
# Determine z position between each fuel pellet, omitting the surfaces
# corresponding to the very bottom and top of the active fuel length
axial_splits = np.linspace(bottom_fuel_stack, top_active_core, num_axial + 1)[1:-1]
axial_surfs = [openmc.ZPlane(z0=z) for z in axial_splits]
if ring_radii is not None:
if num_rings > 1:
# Get z-cylinder surfaces for each ring
rings = []
for i, r in enumerate(ring_radii):
cyl = openmc.ZCylinder(r=r, name='fuel ring {}'.format(i))
for i in range(1, num_rings):
R = sqrt(i*pellet_OR**2/num_rings)
cyl = openmc.ZCylinder(R=R, name='fuel ring {}'.format(i))
rings.append(cyl)
def subdivided_fuel(fill):
# Create universe for UO2 alone with axial/radial subdivision
uo2_cells = []
if num_axial > 1:
for axial_region in subdivide(axial_surfs):
if ring_radii is not None:
for ring_region in subdivide(rings):
cell = openmc.Cell(fill=fill, region=axial_region & ring_region)
uo2_cells.append(cell)
else:
uo2_cells.append(openmc.Cell(fill=fill, region=axial_region))
"""Create universe for UO2 alone with axial/radial subdivision"""
# Create universe with radial rings
univ = openmc.Universe()
if num_rings > 1:
for ring_region in subdivide(rings):
cell = openmc.Cell(fill=fill, region=ring_region)
univ.add_cell(cell)
else:
if ring_radii is not None:
for ring_region in subdivide(rings):
cell = openmc.Cell(fill=fill, region=ring_region)
uo2_cells.append(cell)
else:
raise RuntimeError("Shouldn't call with 1 ring and 1 axial segment")
univ.add_cell(openmc.Cell(fill=fill))
return openmc.Universe(cells=uo2_cells)
# If multiple axial segments are needed, use a lattice to avoid
# performance cost of neighbor lookups
if num_axial > 1:
zlattice = openmc.RectLattice()
zlattice.lower_left = (-100., -100., bottom_fuel_stack)
zlattice.pitch = (200., 200., (top_active_core - bottom_fuel_stack)/num_axial)
zlattice.dimension = (1, 1, num_axial)
zlattice.universes = np.full((num_axial, 1, 1), univ)
zlattice.outer = univ
return zlattice
else:
return univ
# If rings/axial segments are present, create a universe for the subdivided
# fuel. Otherwise just use a plain material.
if ring_radii is not None or num_axial > 1:
if num_rings > 1 or num_axial > 1:
fuel_fill = subdivided_fuel(mats['UO2 1.6 {}'.format(fuel)])
else:
fuel_fill = mats['UO2 1.6 {}'.format(fuel)]
outside_pin_surfaces = [surfs['clad IR'], surfs['clad OR']]
outside_pin_mats = [mats['He'], mats['M5'], mats['H2O']]
outside_pin_mats = [mats['He'], mats['Zr'], mats['H2O']]
univs['Outside pin'] = make_pin(
'Outside pin',
@ -569,12 +759,6 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
materials=outside_pin_mats,
grid='intermediate')
univs['Fuel pin (1.6%) no grid'] = make_pin(
'Pin no grid',
surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
materials=[fuel_fill] + outside_pin_mats
)
# Stack all axial pieces of 1.6% enriched fuel pin cell
within_fuel_surfs = [
@ -642,17 +826,11 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
# If rings/axial segments are present, create a universe for the subdivided
# fuel. Otherwise just use a plain material.
if ring_radii is not None or num_axial > 1:
if num_rings > 1 or num_axial > 1:
fuel_fill = subdivided_fuel(mats['UO2 2.4 {}'.format(fuel)])
else:
fuel_fill = mats['UO2 2.4 {}'.format(fuel)]
univs['Fuel pin (2.4%) no grid'] = make_pin(
'Pin no grid',
surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
materials=[fuel_fill] + outside_pin_mats
)
# Stack all axial pieces of 2.4% enriched fuel pin cell
univs['Fuel pin (2.4%) stack'] = make_pin_stack(
@ -694,26 +872,11 @@ def pin_universes(ring_radii=None, num_axial=196, depleted=False):
# If rings/axial segments are present, create a universe for the subdivided
# fuel. Otherwise just use a plain material.
if ring_radii is not None or num_axial > 1:
if num_rings > 1 or num_axial > 1:
fuel_fill = subdivided_fuel(mats['UO2 3.1 {}'.format(fuel)])
else:
fuel_fill = mats['UO2 3.1 {}'.format(fuel)]
if num_axial > 1:
water_cells = []
for i, r in enumerate(subdivide(axial_surfs)):
cell = openmc.Cell(fill=mats['H2O'], region=r, name=f'Water ({i})')
water_cells.append(cell)
water_fill = openmc.Universe(cells=water_cells)
else:
water_fill = mats['H2O']
univs['Fuel pin (3.1%) no grid'] = make_pin(
'Pin no grid',
surfaces=[surfs['pellet OR']] + outside_pin_surfaces,
materials=[fuel_fill, mats['He'], mats['M5'], water_fill]
)
# Stack all axial pieces of 3.1% enriched fuel pin cell
univs['Fuel pin (3.1%) stack'] = make_pin_stack(

5
smr/smr/reflector.py
View file

@ -12,7 +12,7 @@ converted to actual dimensions by scaling according to the width of an assembly.
import openmc
from .materials import mats
from smr import surfaces
from .surfaces import lattice_pitch
def make_reflector(name, parameters):
@ -34,7 +34,7 @@ def make_reflector(name, parameters):
"""
water_holes = []
for x, y, r in parameters:
zcyl = openmc.ZCylinder(x0=x, y0=y, r=r)
zcyl = openmc.ZCylinder(x0=x, y0=y, R=r)
hole = openmc.Cell(fill=mats['H2O'], region=-zcyl)
water_holes.append(hole)
@ -91,7 +91,6 @@ def reflector_universes():
# All pixel widths are scaled according to the actual width of an assembly
# divided by the width of an assembly in pixels
lattice_pitch = surfaces.lattice_pitch
scale = lattice_pitch/width
# Physical positions

177
smr/smr/surfaces.py
View file

@ -13,6 +13,7 @@ NuScale DC application, chapter 4: https://www.nrc.gov/docs/ML1701/ML17013A274.p
"""
import copy
from math import tan, pi
import numpy as np
@ -39,7 +40,7 @@ clad_OR = 0.374*INCHES/2 # ML17013A274, Table 4.1-2
active_fuel_length = 78.74*INCHES # ML17013A274, Figure 4.2-10
plenum_length = 5.311*INCHES # ML17013A274, Figure 4.2-10
fuel_rod_length = 85.00*INCHES # ML17013A274, Table 4.1-2
lower_end_cap_length = 0.575*INCHES # ML17007A001, Table 3-2
lower_end_cap = 0.575*INCHES # ML17007A001, Table 3-2
# pin cell parameters
guide_tube_IR = 0.450*INCHES/2 # ML17013A274, Table 4.1-2
@ -78,20 +79,19 @@ top_nozzle_width = 8.406*INCHES # ML17013A274, Figure 4.2-2
core_barrel_IR = 74*INCHES/2 # ML17013A274, Table 4.1-2
core_barrel_OR = 78*INCHES/2 # ML17013A274, Table 4.1-2
neutron_shield_OR = core_barrel_OR + 2.0
rpv_IR = 96.5*INCHES/2 # ML17013A274, Table 5.3-1
rpv_OR = 105*INCHES/2 # ML17013A274, Table 5.3-1
rpv_IR = 120.0 # Estimate?
rpv_OR = 135.0 # Estimate?
# axial parameters
reference_z = -36.6205
lowest_extent = reference_z
bottom_support_plate = lowest_extent + 20.000
top_support_plate = bottom_support_plate + 5.000
bottom_lower_nozzle = bottom_support_plate + 5.000
top_lower_nozzle = bottom_lower_nozzle + 4.0*INCHES
bottom_fuel_rod = bottom_lower_nozzle + 4.0*INCHES
top_lower_thimble = bottom_fuel_rod + lower_end_cap_length
bottom_fuel_stack = bottom_fuel_rod + lower_end_cap_length
bot_burn_abs = bottom_fuel_stack + 2.0*INCHES
lowest_extent = 0.000
bottom_support_plate = 20.000
top_support_plate = 25.000
bottom_lower_nozzle = 25.000
top_lower_nozzle = 35.160
bottom_fuel_rod = 35.160
top_lower_thimble = 36.007
bottom_fuel_stack = 36.007
bot_burn_abs = 41.087
top_active_core = bottom_fuel_stack + active_fuel_length
top_plenum = top_active_core + plenum_length
top_fuel_rod = bottom_fuel_rod + fuel_rod_length
@ -127,47 +127,74 @@ neutron_shield_NEtop_SWbot = tan(-pi/6)
surfs = {}
surfs['pellet OR'] = openmc.ZCylinder(r=pellet_OR, name='Pellet OR')
surfs['plenum spring OR'] = openmc.ZCylinder(r=plenum_spring_OR, name='FR Plenum Spring OR')
surfs['clad IR'] = openmc.ZCylinder(r=clad_IR, name='Clad IR')
surfs['clad OR'] = openmc.ZCylinder(r=clad_OR, name='Clad OR')
surfs['GT IR'] = openmc.ZCylinder(r=guide_tube_IR, name='GT IR (above dashpot)')
surfs['GT OR'] = openmc.ZCylinder(r=guide_tube_OR, name='GT OR (above dashpot)')
surfs['GT dashpot IR'] = openmc.ZCylinder(r=guide_tube_dash_IR, name='GT IR (at dashpot)')
surfs['GT dashpot OR'] = openmc.ZCylinder(r=guide_tube_dash_OR, name='GT OR (at dashpot)')
surfs['CP OR'] = openmc.ZCylinder(r=boron_carbide_OR, name='Control Poison OR')
surfs['CR IR'] = openmc.ZCylinder(r=control_rod_IR, name='CR Clad IR')
surfs['CR OR'] = openmc.ZCylinder(r=control_rod_OR, name='CR Clad OR')
surfs['BA IR 1'] = openmc.ZCylinder(r=burn_abs_r1, name='BA IR 1')
surfs['BA IR 2'] = openmc.ZCylinder(r=burn_abs_r2, name='BA IR 2')
surfs['BA IR 3'] = openmc.ZCylinder(r=burn_abs_r3, name='BA IR 3')
surfs['BA IR 4'] = openmc.ZCylinder(r=burn_abs_r4, name='BA IR 4')
surfs['BA IR 5'] = openmc.ZCylinder(r=burn_abs_r5, name='BA IR 5')
surfs['BA IR 6'] = openmc.ZCylinder(r=burn_abs_r6, name='BA IR 6')
surfs['BA IR 7'] = openmc.ZCylinder(r=burn_abs_r7, name='BA IR 7')
surfs['BA IR 8'] = openmc.ZCylinder(r=burn_abs_r8, name='BA IR 8')
surfs['IT IR'] = surfs['BA IR 5']
surfs['IT OR'] = surfs['BA IR 6']
surfs['pellet OR'] = openmc.ZCylinder(
R=pellet_OR, name='Pellet OR')
surfs['plenum spring OR'] = openmc.ZCylinder(
R=plenum_spring_OR, name='FR Plenum Spring OR')
surfs['clad IR'] = openmc.ZCylinder(
R=clad_IR, name='Clad IR')
surfs['clad OR'] = openmc.ZCylinder(
R=clad_OR, name='Clad OR')
surfs['GT IR'] = openmc.ZCylinder(
R=guide_tube_IR, name='GT IR (above dashpot)')
surfs['GT OR'] = openmc.ZCylinder(
R=guide_tube_OR, name='GT OR (above dashpot)')
surfs['GT dashpot IR'] = openmc.ZCylinder(
R=guide_tube_dash_IR, name='GT IR (at dashpot)')
surfs['GT dashpot OR'] = openmc.ZCylinder(
R=guide_tube_dash_OR, name='GT OR (at dashpot)')
surfs['CP OR'] = openmc.ZCylinder(
R=boron_carbide_OR, name='Control Poison OR')
surfs['CR IR'] = openmc.ZCylinder(
R=control_rod_IR, name='CR Clad IR')
surfs['CR OR'] = openmc.ZCylinder(
R=control_rod_OR, name='CR Clad OR')
surfs['BA IR 1'] = openmc.ZCylinder(
R=burn_abs_r1, name='BA IR 1')
surfs['BA IR 2'] = openmc.ZCylinder(
R=burn_abs_r2, name='BA IR 2')
surfs['BA IR 3'] = openmc.ZCylinder(
R=burn_abs_r3, name='BA IR 3')
surfs['BA IR 4'] = openmc.ZCylinder(
R=burn_abs_r4, name='BA IR 4')
surfs['BA IR 5'] = openmc.ZCylinder(
R=burn_abs_r5, name='BA IR 5')
surfs['BA IR 6'] = openmc.ZCylinder(
R=burn_abs_r6, name='BA IR 6')
surfs['BA IR 7'] = openmc.ZCylinder(
R=burn_abs_r7, name='BA IR 7')
surfs['BA IR 8'] = openmc.ZCylinder(
R=burn_abs_r8, name='BA IR 8')
surfs['IT IR'] = copy.deepcopy(surfs['BA IR 5'])
surfs['IT OR'] = copy.deepcopy(surfs['BA IR 6'])
# Rectangular prisms for grid spacers
surfs['rod grid box'] = openmc.rectangular_prism(rod_grid_side, rod_grid_side)
surfs['rod grid box'] = \
openmc.get_rectangular_prism(rod_grid_side, rod_grid_side)
# Rectangular prisms for lattice grid sleeves
surfs['lat grid box inner'] = openmc.rectangular_prism(17.*pin_pitch, 17.*pin_pitch)
surfs['lat grid box outer'] = openmc.rectangular_prism(grid_strap_side, grid_strap_side)
surfs['lat grid box inner'] = \
openmc.get_rectangular_prism(17.*pin_pitch, 17.*pin_pitch)
surfs['lat grid box outer'] = \
openmc.get_rectangular_prism(grid_strap_side, grid_strap_side)
surfs['bot support plate'] = openmc.ZPlane(z0=bottom_support_plate, name='bot support plate')
surfs['top support plate'] = openmc.ZPlane(z0=top_support_plate, name='top support plate')
surfs['bot support plate'] = openmc.ZPlane(
z0=bottom_support_plate, name='bot support plate')
surfs['top support plate'] = openmc.ZPlane(
z0=top_support_plate, name='top support plate')
surfs['bottom FR'] = openmc.ZPlane(z0=bottom_fuel_rod, name='bottom FR')
surfs['top lower nozzle'] = surfs['bottom FR']
surfs['bot lower nozzle'] = surfs['top support plate']
surfs['top lower nozzle'] = copy.deepcopy(surfs['bottom FR'])
surfs['bot lower nozzle'] = copy.deepcopy(surfs['top support plate'])
# axial surfaces
surfs['bot active core'] = openmc.ZPlane(z0=bottom_fuel_stack, name='bot active core')
surfs['top active core'] = openmc.ZPlane(z0=top_active_core, name='top active core')
surfs['bot active core'] = openmc.ZPlane(
z0=bottom_fuel_stack, name='bot active core')
surfs['top active core'] = openmc.ZPlane(
z0=top_active_core, name='top active core')
surfs['top lower thimble'] = surfs['bot active core']
surfs['BA bot'] = openmc.ZPlane(z0=bot_burn_abs, name='bottom of BA')
surfs['BA bot'] = openmc.ZPlane(
z0=bot_burn_abs, name='bottom of BA')
for i, (bottom, top) in enumerate(zip(grid_bottom, grid_top)):
# Create plane for bottom of spacer grid
@ -180,12 +207,17 @@ for i, (bottom, top) in enumerate(zip(grid_bottom, grid_top)):
name = 'top of grid {}'.format(i + 1)
surfs[key] = openmc.ZPlane(z0=top, name=name)
surfs['dashpot top'] = openmc.ZPlane(z0=step0H, name='top dashpot')
surfs['dashpot top'] = openmc.ZPlane(
z0=step0H, name='top dashpot')
surfs['top pin plenum'] = openmc.ZPlane(z0=top_plenum, name='top pin plenum')
surfs['top FR'] = openmc.ZPlane(z0=top_fuel_rod, name='top FR')
surfs['bot upper nozzle'] = openmc.ZPlane(z0=bottom_upper_nozzle, name='bottom upper nozzle')
surfs['top upper nozzle'] = openmc.ZPlane(z0=top_upper_nozzle, name='top upper nozzle')
surfs['top pin plenum'] = openmc.ZPlane(
z0=top_plenum, name='top pin plenum')
surfs['top FR'] = openmc.ZPlane(
z0=top_fuel_rod, name='top FR')
surfs['bot upper nozzle'] = openmc.ZPlane(
z0=bottom_upper_nozzle, name='bottom upper nozzle')
surfs['top upper nozzle'] = openmc.ZPlane(
z0=top_upper_nozzle, name='top upper nozzle')
# Control rod bank surfaces for ARO configuration
for bank in ['A','B','C','D','E',]:
@ -194,37 +226,50 @@ for bank in ['A','B','C','D','E',]:
surfs['bankS{} bot'.format(bank)] = openmc.ZPlane(
z0=step248H, name='CR bankS{} bottom'.format(bank))
surfs['bankA top'] = openmc.ZPlane(z0=bank_top, name='CR bank A top')
surfs['bankA bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank A bottom')
surfs['bankB top'] = openmc.ZPlane(z0=bank_top, name='CR bank B top')
surfs['bankB bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank B bottom')
surfs['bankC top'] = openmc.ZPlane(z0=bank_top, name='CR bank C top')
surfs['bankC bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank C bottom')
surfs['bankD top'] = openmc.ZPlane(z0=bank_top, name='CR bank D top')
surfs['bankD bot'] = openmc.ZPlane(z0=bank_bot, name='CR bank D bottom')
surfs['bankA top'] = openmc.ZPlane(
z0=bank_top, name='CR bank A top')
surfs['bankA bot'] = openmc.ZPlane(
z0=bank_bot, name='CR bank A bottom')
surfs['bankB top'] = openmc.ZPlane(
z0=bank_top, name='CR bank B top')
surfs['bankB bot'] = openmc.ZPlane(
z0=bank_bot, name='CR bank B bottom')
surfs['bankC top'] = openmc.ZPlane(
z0=bank_top, name='CR bank C top')
surfs['bankC bot'] = openmc.ZPlane(
z0=bank_bot, name='CR bank C bottom')
surfs['bankD top'] = openmc.ZPlane(
z0=bank_top, name='CR bank D top')
surfs['bankD bot'] = openmc.ZPlane(
z0=bank_bot, name='CR bank D bottom')
# outer radial surfaces
surfs['core barrel IR'] = openmc.ZCylinder(r=core_barrel_IR, name='core barrel IR')
surfs['core barrel OR'] = openmc.ZCylinder(r=core_barrel_OR, name='core barrel OR')
surfs['neutron shield OR'] = openmc.ZCylinder(r=neutron_shield_OR, name='neutron shield OR')
surfs['core barrel IR'] = openmc.ZCylinder(
R=core_barrel_IR, name='core barrel IR')
surfs['core barrel OR'] = openmc.ZCylinder(
R=core_barrel_OR, name='core barrel OR')
surfs['neutron shield OR'] = openmc.ZCylinder(
R=neutron_shield_OR, name='neutron shield OR')
# neutron shield planes
surfs['neutron shield NWbot SEtop'] = openmc.Plane(
a=1., b=neutron_shield_NWbot_SEtop, c=0., d=0.,
A=1., B=neutron_shield_NWbot_SEtop, C=0., D=0.,
name='neutron shield NWbot SEtop')
surfs['neutron shield NWtop SEbot'] = openmc.Plane(
a=1., b=neutron_shield_NWtop_SEbot, c=0., d=0.,
A=1., B=neutron_shield_NWtop_SEbot, C=0., D=0.,
name='neutron shield NWtop SEbot')
surfs['neutron shield NEbot SWtop'] = openmc.Plane(
a=1., b=neutron_shield_NEbot_SWtop, c=0., d=0.,
A=1., B=neutron_shield_NEbot_SWtop, C=0., D=0.,
name='neutron shield NEbot SWtop')
surfs['neutron shield NEtop SWbot'] = openmc.Plane(
a=1., b=neutron_shield_NEtop_SWbot, c=0., d=0.,
A=1., B=neutron_shield_NEtop_SWbot, C=0., D=0.,
name='neutron shield NEtop SWbot')
# outer radial surfaces
surfs['RPV IR'] = openmc.ZCylinder(r=rpv_IR, name='RPV IR')
surfs['RPV OR'] = openmc.ZCylinder(r=rpv_OR, name='RPV OR', boundary_type='vacuum')
surfs['RPV IR'] = openmc.ZCylinder(
R=rpv_IR, name='RPV IR')
surfs['RPV OR'] = openmc.ZCylinder(
R=rpv_OR, name='RPV OR', boundary_type='vacuum')
# outer axial surfaces
surfs['upper bound'] = openmc.ZPlane(