diff --git a/smr/build-assembly-long.py b/smr/build-assembly-long.py deleted file mode 100644 index 623907c..0000000 --- a/smr/build-assembly-long.py +++ /dev/null @@ -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')) diff --git a/smr/build-assembly-short.py b/smr/build-assembly-short.py deleted file mode 100644 index dc13456..0000000 --- a/smr/build-assembly-short.py +++ /dev/null @@ -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')) diff --git a/smr/build-assembly.py b/smr/build-assembly.py index d4bc166..ad1fb03 100644 --- a/smr/build-assembly.py +++ b/smr/build-assembly.py @@ -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')) diff --git a/smr/build-core-fresh.py b/smr/build-core-fresh.py index 9747ae8..b53fd4f 100644 --- a/smr/build-core-fresh.py +++ b/smr/build-core-fresh.py @@ -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')) diff --git a/smr/build-core-long.py b/smr/build-core-long.py deleted file mode 100644 index 4033d6e..0000000 --- a/smr/build-core-long.py +++ /dev/null @@ -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') diff --git a/smr/build-core-short.py b/smr/build-core-short.py deleted file mode 100644 index 2f5339d..0000000 --- a/smr/build-core-short.py +++ /dev/null @@ -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') diff --git a/smr/milestones.md b/smr/milestones.md deleted file mode 100644 index bb89864..0000000 --- a/smr/milestones.md +++ /dev/null @@ -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). diff --git a/smr/smr/__init__.py b/smr/smr/__init__.py index 027517e..e69de29 100644 --- a/smr/smr/__init__.py +++ b/smr/smr/__init__.py @@ -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 diff --git a/smr/smr/assemblies.py b/smr/smr/assemblies.py index 6f10aaa..dd71d97 100644 --- a/smr/smr/assemblies.py +++ b/smr/smr/assemblies.py @@ -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 diff --git a/smr/smr/core.py b/smr/smr/core.py index 81763fe..4d80fbf 100644 --- a/smr/smr/core.py +++ b/smr/smr/core.py @@ -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) diff --git a/smr/smr/materials.py b/smr/smr/materials.py index ee6c754..8487b3a 100644 --- a/smr/smr/materials.py +++ b/smr/smr/materials.py @@ -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 diff --git a/smr/smr/pins.py b/smr/smr/pins.py index 3966997..1062166 100644 --- a/smr/smr/pins.py +++ b/smr/smr/pins.py @@ -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( diff --git a/smr/smr/reflector.py b/smr/smr/reflector.py index 74aa4bd..2ac4f44 100644 --- a/smr/smr/reflector.py +++ b/smr/smr/reflector.py @@ -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 diff --git a/smr/smr/surfaces.py b/smr/smr/surfaces.py index b454c61..b8cd90b 100644 --- a/smr/smr/surfaces.py +++ b/smr/smr/surfaces.py @@ -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(