diff --git a/smr/build-assembly-long.py b/smr/build-assembly-long.py new file mode 100644 index 0000000..623907c --- /dev/null +++ b/smr/build-assembly-long.py @@ -0,0 +1,186 @@ +#!/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 new file mode 100644 index 0000000..dc13456 --- /dev/null +++ b/smr/build-assembly-short.py @@ -0,0 +1,192 @@ +#!/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 ad1fb03..d4bc166 100644 --- a/smr/build-assembly.py +++ b/smr/build-assembly.py @@ -1,23 +1,28 @@ #!/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 -from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core +from smr.materials import materials, clone +from smr.surfaces import surfs, lattice_pitch, bottom_fuel_stack, top_active_core, pellet_OR from smr.assemblies import assembly_universes -from smr.plots import assembly_plots +from smr import inlet_temperature # Define command-line options parser = argparse.ArgumentParser() -parser.add_argument('-m', '--multipole', action='store_true', - help='Whether to use multipole cross sections') +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('-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, @@ -27,6 +32,7 @@ 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 @@ -40,26 +46,22 @@ else: directory.mkdir(exist_ok=True) # Define geometry with a single assembly -assembly = assembly_universes(args.rings, args.axial, args.depleted) -lattice_sides = openmc.model.get_rectangular_prism(lattice_pitch, lattice_pitch, - boundary_type='reflective') +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') main_cell = openmc.Cell( - fill=assembly['Assembly (3.1%) 16BA'], + fill=assembly['Assembly (3.1%)'], 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.tallies == 'mat': +if args.clone: # Count the number of instances for each cell and material geometry.determine_paths(instances_only=True) @@ -99,10 +101,15 @@ settings.inactive = 100 settings.particles = 10000 settings.output = {'tallies': False, 'summary': False} settings.source = source -settings.sourcepoint_write = False - +settings.sourcepoint = {'write': False} +settings.temperature = { + 'default': inlet_temperature, + 'method': 'interpolation', + 'range': (300.0, 1500.0), +} if args.multipole: - settings.temperature = {'multipole': True, 'tolerance': 1000} + settings.temperature['multipole'] = True + settings.temperature['tolerance'] = 1000 settings.export_to_xml(str(directory / 'settings.xml')) @@ -136,7 +143,3 @@ 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 0a48194..9747ae8 100644 --- a/smr/build-core-fresh.py +++ b/smr/build-core-fresh.py @@ -1,26 +1,30 @@ #!/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 smr.materials import materials -from smr.plots import core_plots -from smr.surfaces import lattice_pitch, bottom_fuel_stack, top_active_core +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 from smr.core import core_geometry +from smr import inlet_temperature # Define command-line options parser = argparse.ArgumentParser() -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('--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('-r', '--rings', type=int, default=10, help='Number of annular regions in fuel') parser.add_argument('-a', '--axial', type=int, default=196, @@ -28,6 +32,7 @@ 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 @@ -40,20 +45,44 @@ else: directory = args.output_dir directory.mkdir(exist_ok=True) -geometry = core_geometry(args.rings, args.axial, args.depleted) +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) -#### "Differentiate" the geometry if using distribmats -if args.tallies == 'mat': - # Count the number of instances for each cell and material +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) - # Extract all cells filled by a fuel material - fuel_mats = {m for m in materials if 'UO2 Fuel' in m.name} +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 - for cell in geometry.get_all_cells().values(): - if cell.fill in fuel_mats: - # Fill cell with list of "differentiated" materials - cell.fill = [cell.fill.clone() for i in range(cell.num_instances)] + # 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 #### Create OpenMC "materials.xml" file all_materials = geometry.get_all_materials() @@ -79,45 +108,14 @@ settings.inactive = 100 settings.particles = 10000 settings.output = {'tallies': False, 'summary': False} settings.source = source -settings.sourcepoint_write = False - +settings.sourcepoint = {'write': False} +settings.temperature = { + 'default': inlet_temperature, + 'method': 'interpolation', + 'range': (300.0, 1500.0), +} if args.multipole: - settings.temperature = {'multipole': True, 'tolerance': 1000} + settings.temperature['multipole'] = True + settings.temperature['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 new file mode 100644 index 0000000..4033d6e --- /dev/null +++ b/smr/build-core-long.py @@ -0,0 +1,118 @@ +#!/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 new file mode 100644 index 0000000..2f5339d --- /dev/null +++ b/smr/build-core-short.py @@ -0,0 +1,131 @@ +#!/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 new file mode 100644 index 0000000..bb89864 --- /dev/null +++ b/smr/milestones.md @@ -0,0 +1,20 @@ +# 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 e69de29..027517e 100644 --- a/smr/smr/__init__.py +++ b/smr/smr/__init__.py @@ -0,0 +1,6 @@ +_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 dd71d97..6f10aaa 100644 --- a/smr/smr/assemblies.py +++ b/smr/smr/assemblies.py @@ -120,13 +120,14 @@ def make_assembly(name, universes): return universe -def assembly_universes(num_rings, num_axial, depleted): +def assembly_universes(ring_radii, num_axial, depleted): """Generate universes for SMR fuel assemblies. Parameters ---------- - num_rings : int - Number of annual regions in fuel + ring_radii : iterable of float + Radii of rings in fuel (note that this doesn't need to include the + full fuel pin radius) num_axial : int Number of axial subdivisions in fuel depleted : bool @@ -138,15 +139,14 @@ def assembly_universes(num_rings, num_axial, depleted): Dictionary mapping a universe name to a openmc.Universe object """ - pins = pin_universes(num_rings, num_axial, depleted) + pins = pin_universes(ring_radii, num_axial, depleted) # Create dictionary to store assembly universes univs = {} # commonly needed universes - gtu = pins['GT empty'] + gtu = pins['GT empty stack'] 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(num_rings, num_axial, depleted): gtu, gtu, gtu, gtu, gtu, gtu, gtu, gtu, gtu, gtu ] - univs['Assembly (2.4%) no BAs' + comment] = \ + univs['Assembly (2.4%)' + comment] = \ make_assembly('Assembly (2.4%) no BAs' + comment, universes) # WITH CONTROL ROD D BANK @@ -231,33 +231,6 @@ def assembly_universes(num_rings, 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')]: @@ -288,134 +261,5 @@ def assembly_universes(num_rings, 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 4d80fbf..81763fe 100644 --- a/smr/smr/core.py +++ b/smr/smr/core.py @@ -5,18 +5,19 @@ 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(num_rings, num_axial, depleted): +def core_geometry(ring_radii, num_axial, depleted): """Generate full core SMR geometry. Parameters ---------- - num_rings : int - Number of annual regions in fuel + ring_radii : iterable of float + Radii of rings in fuel (note that this doesn't need to include the + full fuel pin radius) num_axial : int Number of axial subdivisions in fuel depleted : bool @@ -28,11 +29,12 @@ def core_geometry(num_rings, num_axial, depleted): SMR full core geometry """ - assembly = assembly_universes(num_rings, num_axial, depleted) + assembly = assembly_universes(ring_radii, 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)) @@ -45,67 +47,67 @@ def core_geometry(num_rings, num_axial, depleted): universes[1, 1] = reflector['1,1'] universes[1, 2] = reflector['NW'] - 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, 3] = assembly['Assembly (3.1%)'] + universes[1, 4] = assembly['Assembly (3.1%)'] + universes[1, 5] = assembly['Assembly (3.1%)'] 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%) 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, 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, 7] = reflector['NE'] universes[2, 8] = reflector['2,8'] universes[3, 0] = reflector['3,0'] - 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, 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, 8] = reflector['3,8'] universes[4, 0] = reflector['4,0'] - 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, 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, 8] = reflector['4,8'] universes[5, 0] = reflector['5,0'] - 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, 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, 8] = reflector['5,8'] universes[6, 0] = reflector['6,0'] universes[6, 1] = reflector['SW'] - 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, 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, 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%) instr'] - universes[7, 4] = assembly['Assembly (2.4%) CR D'] - universes[7, 5] = assembly['Assembly (3.1%) instr'] + universes[7, 3] = assembly['Assembly (3.1%)'] + universes[7, 4] = assembly['Assembly (3.1%)'] + universes[7, 5] = assembly['Assembly (3.1%)'] universes[7, 6] = reflector['SE'] universes[7, 7] = reflector['7,7'] @@ -118,6 +120,7 @@ def core_geometry(num_rings, 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') @@ -133,75 +136,10 @@ def core_geometry(num_rings, 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['neutron shield OR'] & -surfs['RPV IR'] & + cell.region = (+surfs['core barrel 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 8487b3a..ee6c754 100644 --- a/smr/smr/materials.py +++ b/smr/smr/materials.py @@ -1,7 +1,11 @@ """Instantiate the OpenMC Materials needed by the core model.""" +import copy + import openmc -from openmc.data import atomic_weight, atomic_mass +from openmc.data import atomic_weight, atomic_mass, water_density + +from . import system_pressure, core_average_temperature _DEPLETION_NUCLIDES = [ @@ -41,13 +45,11 @@ 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') @@ -56,7 +58,6 @@ 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') @@ -64,9 +65,17 @@ 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') @@ -76,7 +85,6 @@ 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') @@ -97,7 +105,6 @@ 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') @@ -105,9 +112,18 @@ 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') @@ -116,10 +132,11 @@ mats['AIC'].add_element('Cd', 0.05, 'wo') #### Borated Water -boron_ppm = 975 +# Concentration of boron at beginning of equilibrium cycle +boron_ppm = 1240 # ML17013A274, Figure 4.3-17 -# Density of clean water at 2250 psia T=560F NIST -h2o_dens = 0.73986 +# Density of water +h2o_dens = water_density(core_average_temperature, system_pressure) # Weight percent of natural boron in borated water wB_Bh2o = boron_ppm * 1.0e-6 @@ -146,7 +163,6 @@ 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') @@ -178,7 +194,6 @@ 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') @@ -191,7 +206,6 @@ 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) @@ -199,7 +213,6 @@ 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) @@ -207,7 +220,6 @@ 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) @@ -215,7 +227,6 @@ 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) @@ -224,7 +235,6 @@ 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) @@ -233,7 +243,6 @@ 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) @@ -245,3 +254,10 @@ 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 fce0bab..3966997 100644 --- a/smr/smr/pins.py +++ b/smr/smr/pins.py @@ -7,7 +7,7 @@ import openmc from openmc.model import subdivide from .materials import mats -from .surfaces import surfs, pellet_OR, bottom_fuel_rod, top_active_core +from .surfaces import surfs, pellet_OR, bottom_fuel_stack, top_active_core def make_pin(name, surfaces, materials, grid=None): @@ -142,13 +142,14 @@ def make_pin_stack(name, zsurfaces, universes, boundary, fuel_fill): return universe -def pin_universes(num_rings=10, num_axial=196, depleted=False): +def pin_universes(ring_radii=None, num_axial=196, depleted=False): """Generate universes for SMR fuel pins. Parameters ---------- - num_rings : int - Number of annual regions in fuel + ring_radii : iterable of float + Radii of rings in fuel (note that this doesn't need to include the + full fuel pin radius) num_axial : int Number of axial subdivisions in fuel depleted : bool @@ -230,8 +231,8 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False): surfs['top upper nozzle'] ] - univs['GT empty'] = make_stack( - 'GT empty', surfaces=stack_surfs, + univs['GT empty stack'] = make_stack( + 'GT empty stack', surfaces=stack_surfs, universes=[univs['water pin'], univs['water pin'], univs['water pin'], @@ -472,199 +473,6 @@ def pin_universes(num_rings=10, 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', @@ -674,16 +482,16 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False): univs['end plug'] = make_pin( 'end plug', [surfs['clad OR']], - [mats['Zr'], mats['H2O']]) + [mats['M5'], mats['H2O']]) univs['pin plenum'] = make_pin( 'pin plenum', surfaces=[surfs['plenum spring OR'], surfs['clad IR'], surfs['clad OR']], - materials=[mats['In'], + materials=[mats['SS302'], mats['He'], - mats['Zr'], + mats['M5'], mats['H2O']]) univs['pin plenum grid (intermediate)'] = make_pin( @@ -703,15 +511,14 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False): 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_rod, top_active_core, num_axial + 1)[1:-1] + 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 num_rings > 1: + if ring_radii is not None: # Get z-cylinder surfaces for each ring rings = [] - 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)) + for i, r in enumerate(ring_radii): + cyl = openmc.ZCylinder(r=r, name='fuel ring {}'.format(i)) rings.append(cyl) def subdivided_fuel(fill): @@ -719,14 +526,14 @@ def pin_universes(num_rings=10, num_axial=196, depleted=False): uo2_cells = [] if num_axial > 1: for axial_region in subdivide(axial_surfs): - if num_rings > 1: + 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)) else: - if num_rings > 1: + if ring_radii is not None: for ring_region in subdivide(rings): cell = openmc.Cell(fill=fill, region=ring_region) uo2_cells.append(cell) @@ -737,13 +544,13 @@ def pin_universes(num_rings=10, 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 num_rings > 1 or num_axial > 1: + if ring_radii is not None 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['Zr'], mats['H2O']] + outside_pin_mats = [mats['He'], mats['M5'], mats['H2O']] univs['Outside pin'] = make_pin( 'Outside pin', @@ -762,6 +569,12 @@ def pin_universes(num_rings=10, 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 = [ @@ -829,11 +642,17 @@ def pin_universes(num_rings=10, 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 num_rings > 1 or num_axial > 1: + if ring_radii is not None 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( @@ -875,11 +694,26 @@ def pin_universes(num_rings=10, 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 num_rings > 1 or num_axial > 1: + if ring_radii is not None 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 2ac4f44..74aa4bd 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 .surfaces import lattice_pitch +from smr import surfaces 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,6 +91,7 @@ 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 3aa13a8..b454c61 100644 --- a/smr/smr/surfaces.py +++ b/smr/smr/surfaces.py @@ -13,7 +13,6 @@ NuScale DC application, chapter 4: https://www.nrc.gov/docs/ML1701/ML17013A274.p """ -import copy from math import tan, pi import numpy as np @@ -40,7 +39,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 = 0.575*INCHES # ML17007A001, Table 3-2 +lower_end_cap_length = 0.575*INCHES # ML17007A001, Table 3-2 # pin cell parameters guide_tube_IR = 0.450*INCHES/2 # ML17013A274, Table 4.1-2 @@ -79,19 +78,20 @@ 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 = 120.0 # Estimate? -rpv_OR = 135.0 # Estimate? +rpv_IR = 96.5*INCHES/2 # ML17013A274, Table 5.3-1 +rpv_OR = 105*INCHES/2 # ML17013A274, Table 5.3-1 # axial parameters -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 +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 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,74 +127,47 @@ 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'] = copy.deepcopy(surfs['BA IR 5']) -surfs['IT OR'] = copy.deepcopy(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'] = surfs['BA IR 5'] +surfs['IT OR'] = surfs['BA IR 6'] # Rectangular prisms for grid spacers -surfs['rod grid box'] = \ - openmc.get_rectangular_prism(rod_grid_side, rod_grid_side) +surfs['rod grid box'] = openmc.rectangular_prism(rod_grid_side, rod_grid_side) # Rectangular prisms for lattice grid sleeves -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['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['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'] = copy.deepcopy(surfs['bottom FR']) -surfs['bot lower nozzle'] = copy.deepcopy(surfs['top support plate']) +surfs['top lower nozzle'] = surfs['bottom FR'] +surfs['bot lower nozzle'] = 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'] = copy.deepcopy(surfs['bot active core']) -surfs['BA bot'] = openmc.ZPlane( - z0=bot_burn_abs, name='bottom of BA') +surfs['top lower thimble'] = surfs['bot active core'] +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 @@ -207,17 +180,12 @@ 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',]: @@ -226,50 +194,37 @@ 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(