diff --git a/examples/msbr/_control_rods.py b/examples/msbr/_control_rods.py index 0f9a9fafd..6ed12c8a1 100644 --- a/examples/msbr/_control_rods.py +++ b/examples/msbr/_control_rods.py @@ -1,55 +1,41 @@ import openmc import numpy as np -def shared_cr_geometry(): - fuel_hole = openmc.ZCylinder(r=5.08, name='cr_fuel_hole') - elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary - - gr_sq_neg = openmc.rectangular_prism(7.23646*2, 7.23646*2, corner_radius=0.99) # Graphite square - # params for - r_d = 1.16 - e_d = 2 * r_d / np.sqrt(3) - r_dt = 0.8 - r_c = 0.18 - l1 = 5.8801 - l2 = 6.505 - l3 = 8.03646 - ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c) - br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c) - lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c) - ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c) - ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip') - br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip') - ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip') - lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip') - - gr_corners = elem_bound & (ul | br | lb | ru | - -ul_t | -br_t | -ru_t | -lb_t) - - inter_elem_channel = (~gr_sq_neg & elem_bound & - ~ul & ~br & ~lb & ~ru & - +ul_t & +ru_t & +ru_t & +lb_t) - - - return fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel - -def control_rod(fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel, fuel, moder): +def control_rod(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, fuel_hole, fuel, moder): s1 = openmc.ZCylinder(r=4.7625, name='control_rod') c1 = openmc.Cell(fill=moder, region=-s1, name='control_rod') c2 = openmc.Cell(fill=fuel, region=(+s1 & -fuel_hole), name='cr_fuel_inner') - moderator_block = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg | gr_corners), name='cr_moderator') - fuel_outer = openmc.Cell(fill=fuel, region=inter_elem_channel, name='cr_fuel_outer') + c3 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg & inter_elem_channel), + name='cr_moderator') + c4 = openmc.Cell(fill=fuel, region= (~gr_sq_neg & elem_bound & inter_elem_channel), + name='cr_fuel_outer') + #universe_id=3 - cr = openmc.Universe(name='control_rod') - cr.add_cells([c1, c2, moderator_block, fuel_outer]) + cr = openmc.Universe(name='control_rod', + cells=[c1, c2, c3, c4]) + + for i, (corner, name) in enumerate(gr_corners): + if i < 4: + cr.add_cell(openmc.Cell(fill=moder, region=~corner, name=f'cr_moderator_{name}')) + else: + cr.add_cell(openmc.Cell(fill=moder, region=+corner, name=f'cr_moderator_{name}')) + return cr -def control_rod_channel(fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel, fuel, moder): +def control_rod_channel(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, fuel_hole, fuel, moder): c1 = openmc.Cell(fill=fuel, region=(-fuel_hole), name='crc_fuel_inner') - moderator_block = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg | gr_corners), name='crc_moderator') - fuel_outer = openmc.Cell(fill=fuel, region=inter_elem_channel, name='crc_fuel_outer') + c2 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg & inter_elem_channel), name='crc_moderator') + c3 = openmc.Cell(fill=fuel, region=(~gr_sq_neg & elem_bound & inter_elem_channel), name='crc_fuel_outer') + # universe_id=4 - crc = openmc.Universe(name='control_rod_channel') - crc.add_cells([c1, moderator_block, fuel_outer]) + crc = openmc.Universe(name='control_rod_channel', + cells=[c1, c2, c3]) + + for i, (corner, name) in enumerate(gr_corners): + if i < 4: + crc.add_cell(openmc.Cell(fill=moder, region=~corner, name=f'crc_moderator_{name}')) + else: + crc.add_cell(openmc.Cell(fill=moder, region=+corner, name=f'crc_moderator_{name}')) + return crc diff --git a/examples/msbr/_control_rods_optimized.py b/examples/msbr/_control_rods_optimized.py new file mode 100644 index 000000000..8bb9ade7e --- /dev/null +++ b/examples/msbr/_control_rods_optimized.py @@ -0,0 +1,106 @@ +import openmc +import numpy as np + +def shared_cr_geometry(): + fuel_hole = openmc.ZCylinder(r=5.08, name='cr_fuel_hole') + elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary + eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values()) + + gr_sq_neg = openmc.rectangular_prism(7.23646*2, 7.23646*2, corner_radius=0.99) # Graphite square + # params for + (gr_minx, gr_maxx, gr_miny, gr_maxy, + gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, + gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, + gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values()) + + r_d = 1.16 + e_d = 2 * r_d / np.sqrt(3) + r_dt = 0.8 + r_c = 0.18 + l1 = 5.8801 + l2 = 6.505 + l3 = 8.03646 + ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c) + br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c) + lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c) + ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c) + ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip') + br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip') + ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip') + lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip') + + gr_ul = ul & -eb_maxy & +gr_maxy + gr_ul_fill = ul & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy + gr_br = br & +eb_miny & -gr_miny + gr_br_fill = br & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny + gr_lb = lb & +eb_minx & -gr_minx + gr_lb_fill = lb & +gr_cyl_lb & +gr_minx & -gr_cyl_minx & -gr_cyl_miny + gr_ru = ru & -eb_maxx & +gr_maxx + gr_ru_fill = ru & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxx & +gr_cyl_maxy + + gr_ul_t = -ul_t & +eb_miny & -gr_miny + gr_br_t = -br_t & -eb_maxy & +gr_maxy + gr_ru_t = -ru_t & +eb_minx & -gr_minx + gr_lb_t = -lb_t & -eb_maxx & +gr_maxx + + gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill) + gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t) + + inter_elem_channel = (~gr_sq_neg & elem_bound & + ~ul & ~br & ~lb & ~ru & + +ul_t & +br_t & +ru_t & +lb_t) + + + return fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel + +def control_rod(fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, fuel, moder): + s1 = openmc.ZCylinder(r=4.7625, name='control_rod') + + c1 = openmc.Cell(fill=moder, region=-s1, name='control_rod') + c2 = openmc.Cell(fill=fuel, region=(+s1 & -fuel_hole), name='cr_fuel_inner') + c3 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg), name='cr_moderator') + c4 = openmc.Cell(fill=fuel, region=inter_elem_channel, name='cr_fuel_outer') + + c3_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='cr_moderator_ul') + c3_br = openmc.Cell(fill=moder, region=gr_corners[1], name='cr_moderator_br') + c3_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='cr_moderator_ru') + c3_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='cr_moderator_lb') + c3_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='cr_moderator_ul_fill') + c3_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='cr_moderator_br_fill') + c3_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='cr_moderator_ru_fill') + c3_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='cr_moderator_lb_fill') + c3_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='cr_moderator_ul_t') + c3_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='cr_moderator_br_t') + c3_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='cr_moderator_ru_t') + c3_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='cr_moderator_lb_t') + + + #universe_id=3 + cr = openmc.Universe(name='control_rod') + cr.add_cells([c1, c2, c3, c4, c3_ul, c3_br, c3_ru, c3_lb, c3_ulf, c3_brf, c3_ruf, c3_lbf, c3_ul_t, c3_br_t, c3_ru_t, c3_lb_t]) + + return cr + +def control_rod_channel(fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, fuel, moder): + c1 = openmc.Cell(fill=fuel, region=(-fuel_hole), name='crc_fuel_inner') + c2 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq_neg), name='crc_moderator') + c3 = openmc.Cell(fill=fuel, region=inter_elem_channel, name='crc_fuel_outer') + + c2_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='crc_moderator_ul') + c2_br = openmc.Cell(fill=moder, region=gr_corners[1], name='crc_moderator_br') + c2_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='crc_moderator_ru') + c2_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='crc_moderator_lb') + c2_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='crc_moderator_ul_fill') + c2_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='crc_moderator_br_fill') + c2_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='crc_moderator_ru_fill') + c2_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='crc_moderator_lb_fill') + c2_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='crc_moderator_ul_t') + c2_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='crc_moderator_br_t') + c2_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='crc_moderator_ru_t') + c2_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='crc_moderator_lb_t') + + # universe_id=4 + crc = openmc.Universe(name='control_rod_channel') + crc.add_cells([c1, c2, c3, c2_ul, c2_br, c2_ru, c2_lb, c2_ulf, c2_brf, c2_ruf, c2_lbf, c2_ul_t, c2_br_t, c2_ru_t, c2_lb_t]) + + return crc diff --git a/examples/msbr/_core_elements.py b/examples/msbr/_core_elements.py index e2b276814..e254cbe01 100644 --- a/examples/msbr/_core_elements.py +++ b/examples/msbr/_core_elements.py @@ -1,53 +1,33 @@ import openmc import numpy as np - def _bound_zone_cells(cells_tuples, levels): """Helper function that moves Zone IA and Zone IIA cells to the appropriate height.""" cell_list = [] + n_levels = len(cells_tuples) for i, cells in enumerate(cells_tuples): - lower_bound = levels[i] - upper_bound = levels[i+1] + if i == 0: + lower_bound = None + upper_bound = levels[i] + elif i == n_levels - 1: + lower_bound = levels[i-1] + upper_bound = None + else: + lower_bound = levels[i-1] + upper_bound = levels[i] for j, cell in enumerate(cells): - cell.region = cell.region & +lower_bound & -upper_bound + if lower_bound is None: + cell.region = cell.region & -upper_bound + elif upper_bound is None: + cell.region = cell.region & +lower_bound + else: + cell.region = cell.region & +lower_bound & -upper_bound cell_list.append(cell) return cell_list -def shared_elem_geometry(): - """Surfaces and regions shared by Zone IA and Zone IIA elements. - Specs found in Robertson, 1971 Fig 3.4 (p. 17) and Fig 3.5 (p.18)""" - - elem_bound = openmc.rectangular_prism(5.08*2, 5.08*2) # Pin outer boundary - - gr_sq_neg = openmc.rectangular_prism(4.953*2, 4.953*2, corner_radius=0.46) # Graphite square - - # params for main pin section for both I-A and II-A - r_d = 0.66802 - l1 = 4.28498 - l2 = 4.53898 - l3 = 5.62102 - ul = openmc.ZCylinder(-l1, l2, r_d, name='corner_ul') - br = openmc.ZCylinder(l1, -l2, r_d, name='corner_br') - lb = openmc.ZCylinder(-l2, -l1, r_d, name='corner_lb') - ru = openmc.ZCylinder(l2, l1, r_d, name='corner_ru') - ul_t = openmc.ZCylinder(-l1, -l3, r_d, name='corner_ul_tip') - br_t = openmc.ZCylinder(l1, l3, r_d, name='corner_br_tip') - ru_t = openmc.ZCylinder(-l3, l1, r_d, name='corner_ru_tip') - lb_t = openmc.ZCylinder(l3, -l1, r_d, name='corner_lb_tip') - - gr_corners = elem_bound & (-ul | -br | -lb | -ru | - -ul_t | -br_t | -ru_t | -lb_t) - - inter_elem_channel = (~gr_sq_neg & elem_bound & - +ul & +br & +lb & +ru & - +ul_t & +ru_t & +ru_t & +lb_t) - gr_round_4 = openmc.ZCylinder(r=2.2225, name='gr_round_4') - - return elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4 - def zoneIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, moder, fuel, hast): """Zone IA element. Specs found in Robertson, 1971 Fig 3.4 (p. 17)""" - elem_levels = [0.0, 22.86, 419.10, 438.15, 445.135, 449.58] + elem_levels = [22.86, 419.10, 438.15, 445.135] level_bounds = [] for level in elem_levels: level_bounds.append(openmc.ZPlane(z0=level)) @@ -57,7 +37,7 @@ def zoneIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, mo h = 12.66 theta = np.arctan(4.953 / h) r2 = (1 / np.cos(theta))**2 - 1 - s3 = openmc.ZCone(z0=h + elem_levels[3], r2=r2, name='cone_i') + s3 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_i') c1 = openmc.Cell(fill=fuel, region=(-s2), name='ia_fuel_inner_1') c2 = openmc.Cell(fill=moder, region=(+s2 & -s1), name='ia_moderator_1') @@ -67,11 +47,13 @@ def zoneIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, mo # I-A main (lower 2) c4 = c1.clone(clone_materials=False) c4.name = 'ia_fuel_inner_main' - c5 = openmc.Cell(fill=moder, region=(+s2 & - gr_sq_neg | - gr_corners), name='ia_moderator_main') - c6 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='ia_fuel_outer_main') - iam = (c4, c5, c6) + c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq_neg & inter_elem_channel), + name='ia_moderator_main') + c6 = openmc.Cell(fill=fuel, region=(~gr_sq_neg & elem_bound & inter_elem_channel), + name='ia_fuel_outer_main') + iam = [c4, c5, c6] + for (corner, name) in gr_corners: + iam.append(openmc.Cell(fill=moder, region=+corner, name=f'ia_moderator_main_{name}')) # I-A 2 (upper 1) c7 = c1.clone(clone_materials=False) @@ -103,7 +85,7 @@ def zoneIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, mo def zoneIIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, moder, fuel): """Zone IIA element. Specs found in Robertson, 1971 Fig 3.5 (p. 18)""" - elem_levels = [0.0, 434.34, 436.88, 439.42, 441.96, 449.58] + elem_levels = [434.34, 436.88, 439.42, 441.96] level_bounds = [] for level in elem_levels: level_bounds.append(openmc.ZPlane(z0=level)) @@ -113,21 +95,24 @@ def zoneIIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, m h = 6.5 theta = np.arctan(3.65125 / h) r2 = (1 / np.cos(theta))**2 - 1 - s4 = openmc.ZCone(z0=h + elem_levels[3], r2=r2, name='cone_ii') + s4 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_ii') # II-A main (lower 1) c1 = openmc.Cell(fill=fuel, region=(-s1), name='iia_fuel_inner_main') - c2 = openmc.Cell(fill=moder, region=(+s1 & gr_sq_neg | gr_corners), name='iia_moderator_main') - c3 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='iia_fuel_outer_main') - c3.name = 'iia_fuel_outer_main' - iiam = (c1, c2, c3) + c2 = openmc.Cell(fill=moder, region=(+s1 & gr_sq_neg & inter_elem_channel), name='iia_moderator_main') + c3 = openmc.Cell(fill=fuel, region=(~gr_sq_neg & elem_bound & inter_elem_channel), name='iia_fuel_outer_main') + iiam = [c1, c2, c3] # II-A 2 (upper 1) c4 = openmc.Cell(fill=fuel, region=(-s2), name='iia_fuel_inner_2') - c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq_neg | gr_corners), name='iia_moderator_2') + c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq_neg & inter_elem_channel), name='iia_moderator_2') c6 = c3.clone(clone_materials=False) c6.name = 'iia_fuel_outer_2' - iia2 = (c4, c5, c6) + iia2 = [c4, c5, c6] + + for (corner, name) in gr_corners: + iiam.append(openmc.Cell(fill=moder, region=+corner, name=f'iia_moderator_main_{name}')) + iia2.append(openmc.Cell(fill=moder, region=+corner, name=f'iia_moderator_2_{name}')) # II-A 3 (upper 2) c7 = c4.clone(clone_materials=False) diff --git a/examples/msbr/_core_elements_optimized.py b/examples/msbr/_core_elements_optimized.py new file mode 100644 index 000000000..0db560c97 --- /dev/null +++ b/examples/msbr/_core_elements_optimized.py @@ -0,0 +1,266 @@ +import openmc +import numpy as np + +def _bound_zone_cells(cells_tuples, levels): + """Helper function that moves Zone IA and Zone IIA cells to the + appropriate height.""" + cell_list = [] + n_levels = len(cells_tuples) + for i, cells in enumerate(cells_tuples): + if i == 0: + lower_bound = None + upper_bound = levels[i] + elif i == n_levels - 1: + lower_bound = levels[i-1] + upper_bound = None + else: + lower_bound = levels[i-1] + upper_bound = levels[i] + + for j, cell in enumerate(cells): + if lower_bound is None: + cell.region = cell.region & -upper_bound + elif upper_bound is None: + cell.region = cell.region & +lower_bound + else: + cell.region = cell.region & +lower_bound & -upper_bound + cell_list.append(cell) + return cell_list + +def shared_elem_geometry(): + """Surfaces and regions shared by Zone IA and Zone IIA elements. + Specs found in Robertson, 1971 Fig 3.4 (p. 17) and Fig 3.5 (p.18)""" + + elem_bound = openmc.rectangular_prism(5.08*2, 5.08*2) # Pin outer boundary + eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values()) + + gr_sq_neg = openmc.rectangular_prism(4.953*2, 4.953*2, corner_radius=0.46) # Graphite square + (gr_minx, gr_maxx, gr_miny, gr_maxy, + gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, + gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, + gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values()) + + # params for main pin section for both I-A and II-A + r_d = 0.66802 + l1 = 4.28498 + l2 = 4.53898 + l3 = 5.62102 + ul = openmc.ZCylinder(-l1, l2, r_d, name='corner_ul') + br = openmc.ZCylinder(l1, -l2, r_d, name='corner_br') + lb = openmc.ZCylinder(-l2, -l1, r_d, name='corner_lb') + ru = openmc.ZCylinder(l2, l1, r_d, name='corner_ru') + ul_t = openmc.ZCylinder(-l1, -l3, r_d, name='corner_ul_tip') + br_t = openmc.ZCylinder(l1, l3, r_d, name='corner_br_tip') + ru_t = openmc.ZCylinder(-l3, l1, r_d, name='corner_ru_tip') + lb_t = openmc.ZCylinder(l3, -l1, r_d, name='corner_lb_tip') + + gr_ul = -ul & -eb_maxy & +gr_maxy# | + gr_ul_fill = -ul & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy + gr_br = -br & +eb_miny & -gr_miny# | + gr_br_fill = -br & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny + gr_lb = -lb & +eb_minx & -gr_minx# | + gr_lb_fill = -lb & +gr_cyl_lb & +gr_minx & -gr_cyl_minx & -gr_cyl_miny + gr_ru = -ru & -eb_maxx & +gr_maxx# | + gr_ru_fill = -ru & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxx & +gr_cyl_maxy + + gr_ul_t = -ul_t & +eb_miny & -gr_miny + gr_br_t = -br_t & -eb_maxy & +gr_maxy + gr_ru_t = -ru_t & +eb_minx & -gr_minx + gr_lb_t = -lb_t & -eb_maxx & +gr_maxx + + gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill) + gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t) + + inter_elem_channel = (~gr_sq_neg & elem_bound & + +ul & +br & +lb & +ru & + +ul_t & +br_t & +ru_t & +lb_t) + + gr_round_4 = openmc.ZCylinder(r=2.2225, name='gr_round_4') + + return elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4 + +def zoneIA(elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4, moder, fuel, hast): + """Zone IA element. Specs found in Robertson, 1971 Fig 3.4 (p. 17)""" + elem_levels = [22.86, 419.10, 438.15, 445.135] + level_bounds = [] + for level in elem_levels: + level_bounds.append(openmc.ZPlane(z0=level)) + s1 = openmc.ZCylinder(r=4.953, name='ia_gr_round_1') + s2 = openmc.ZCylinder(r=1.71069, name='ia_fuel_hole') + + h = 12.66 + theta = np.arctan(4.953 / h) + r2 = (1 / np.cos(theta))**2 - 1 + s3 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_i') + + c1 = openmc.Cell(fill=fuel, region=(-s2), name='ia_fuel_inner_1') + c2 = openmc.Cell(fill=moder, region=(+s2 & -s1), name='ia_moderator_1') + c3 = openmc.Cell(fill=fuel, region=(+s1 & elem_bound), name='ia_fuel_outer_1') + ia1 = (c1, c2, c3) + + # I-A main (lower 2) + s2 = s2.clone() + gr_sq_neg = gr_sq_neg.clone() + c4 = c1.clone(clone_materials=False) + c4.name = 'ia_fuel_inner_main' + c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq_neg), name='ia_moderator_main') + c6 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='ia_fuel_outer_main') + c5_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='ia_moderator_main_ul') + c5_br = openmc.Cell(fill=moder, region=gr_corners[1], name='ia_moderator_main_br') + c5_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='ia_moderator_main_ru') + c5_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='ia_moderator_main_lb') + c5_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='ia_moderator_main_ul_fill') + c5_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='ia_moderator_main_br_fill') + c5_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='ia_moderator_main_ru_fill') + c5_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='ia_moderator_main_lb_fill') + c5_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='ia_moderator_main_ul_t') + c5_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='ia_moderator_main_br_t') + c5_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='ia_moderator_main_ru_t') + c5_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='ia_moderator_main_lb_t') + + iam = (c4, c5, c6, + c5_ul, c5_br, c5_ru, c5_lb, + c5_ulf, c5_brf, c5_ruf, c5_lbf, + c5_ul_t, c5_br_t, c5_ru_t, c5_lb_t) + + # I-A 2 (upper 1) + c7 = c1.clone(clone_materials=False) + c7.name = 'ia_fuel_inner_2' + c8 = c2.clone(clone_materials=False) + c8.name = 'ia_moderator_2' + c9 = c3.clone(clone_materials=False) + c9.name = 'ia_fuel_outer_2' + ia2 = (c7, c8, c9) + + # I-A 3 (upper 2)' + s2 = s2.clone() + s3 = s3.clone() + elem_bound = elem_bound.clone() + c10 = c1.clone(clone_materials=False) + c10.name = 'ia_fuel_inner_3' + c11 = openmc.Cell(fill=moder, region=(+s2 & -s3), name='ia_moderator_3') + c12 = openmc.Cell(fill=fuel, region=(+s3 & elem_bound), name='ia_fuel_outer_3') + ia3 = (c10, c11, c12) + + # I-A 4 (upper 3) + s2 = s2.clone() + elem_bound = elem_bound.clone() + c13 = openmc.Cell(fill=hast, region=(-s2), name='ia_hast') + c14 = openmc.Cell(fill=moder, region=(+s2 & -gr_round_4), name='ia_moderator_4') + c15 = openmc.Cell(fill=fuel, region=(+gr_round_4 & elem_bound), name='ia_fuel_outer_4') + ia4 = (c13, c14, c15) + + elem_cells = [ia1, iam, ia2, ia3, ia4] + # universe_id=1 + ia = openmc.Universe(name='zone_ia') + ia.add_cells(_bound_zone_cells(elem_cells, level_bounds)) + return ia + +def zoneIIA(elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4, moder, fuel): + """Zone IIA element. Specs found in Robertson, 1971 Fig 3.5 (p. 18)""" + elem_levels = [434.34, 436.88, 439.42, 441.96] + level_bounds = [] + for level in elem_levels: + level_bounds.append(openmc.ZPlane(z0=level)) + s1 = openmc.ZCylinder(r=3.302, name='iia_fuel_hole_main') # Hole with fuel salt - Fig 3.5, Roberton 1971 (3.27787 - p.47) + s2 = openmc.ZCylinder(r=0.635, name='iia_fuel_hole_2') + s3 = openmc.ZCylinder(r=3.65125, name='iia_gr_round_3') + h = 6.5 + theta = np.arctan(3.65125 / h) + r2 = (1 / np.cos(theta))**2 - 1 + s4 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_ii') + + # II-A main (lower 1) + c1 = openmc.Cell(fill=fuel, region=(-s1), name='iia_fuel_inner_main') + c2 = openmc.Cell(fill=moder, region=(+s1 & gr_sq_neg), name='iia_moderator_main') + c3 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='iia_fuel_outer_main') + c3.name = 'iia_fuel_outer_main' + c2_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='iia_moderator_main_ul') + c2_br = openmc.Cell(fill=moder, region=gr_corners[1], name='iia_moderator_main_br') + c2_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='iia_moderator_main_ru') + c2_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='iia_moderator_main_lb') + c2_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='iia_moderator_main_ul_fill') + c2_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='iia_moderator_main_br_fill') + c2_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='iia_moderator_main_ru_fill') + c2_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='iia_moderator_main_lb_fill') + c2_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='iia_moderator_main_ul_t') + c2_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='iia_moderator_main_br_t') + c2_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='iia_moderator_main_ru_t') + c2_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='iia_moderator_main_lb_t') + iiam = (c1, c2, c3, c2_ul, c2_br, c2_ru, c2_lb, c2_ulf, c2_brf, c2_ruf, c2_lbf, c2_ul_t, c2_br_t, c2_ru_t, c2_lb_t) + + # II-A 2 (upper 1) + gr_sq_neg = gr_sq_neg.clone() + c4 = openmc.Cell(fill=fuel, region=(-s2), name='iia_fuel_inner_2') + c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq_neg), name='iia_moderator_2') + c6 = c3.clone(clone_materials=False) + c6.name = 'iia_fuel_outer_2' + c5_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='iia_moderator_2_ul') + c5_br = openmc.Cell(fill=moder, region=gr_corners[1], name='iia_moderator_2_br') + c5_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='iia_moderator_2_ru') + c5_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='iia_moderator_2_lb') + c5_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='iia_moderator_2_ul_fill') + c5_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='iia_moderator_2_br_fill') + c5_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='iia_moderator_2_ru_fill') + c5_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='iia_moderator_2_lb_fill') + c5_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='iia_moderator_2_ul_t') + c5_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='iia_moderator_2_br_t') + c5_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='iia_moderator_2_ru_t') + c5_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='iia_moderator_2_lb_t') + iia2 = (c4, c5, c6, c5_ul, c5_br, c5_ru, c5_lb, c5_ulf, c5_brf, c5_ruf, c5_lbf, c5_ul_t, c5_br_t, c5_ru_t, c5_lb_t) + + # II-A 3 (upper 2) + s2 = s2.clone() + elem_bound = elem_bound.clone() + c7 = c4.clone(clone_materials=False) + c7.name = 'iia_fuel_inner_3' + c8 = openmc.Cell(fill=moder, region=(+s2 & -s3), name='iia_moderator_3') + c9 = openmc.Cell(fill=fuel, region=(+s3 & elem_bound), name='iia_fuel_outer_3') + iia3 = (c7, c8, c9) + + # II-A 4 (upper 3) + elem_bound = elem_bound.clone() + c10 = openmc.Cell(fill=moder, region=(-s4), name='iia_moderator_4') + c11 = openmc.Cell(fill=fuel, region=(+s4 & elem_bound), name='iia_fuel_outer_4') + iia4 = (c10, c11) + + # II-A 5 (upper 4) + elem_bound = elem_bound.clone() + c12 = openmc.Cell(fill=moder, region=(-gr_round_4), name='iia_moderator_5') + c13 = openmc.Cell(fill=fuel, region=(+gr_round_4 & elem_bound), name='iia_fuel_outer_5') + iia5 = (c12, c13) + + elem_cells = [iiam, iia2, iia3, iia4, iia5] + # universe_id=2 + iia = openmc.Universe(name='zone_iia') + iia.add_cells(_bound_zone_cells(elem_cells, level_bounds)) + return iia + +def void_cell(elem_bound): + elem_bound = elem_bound.clone() + c1 = openmc.Cell(region=elem_bound, name='lattice_void') + #universe_id=5 + v = openmc.Universe(name='lattice_void') + v.add_cell(c1) + return v + +def graphite_triangles(elem_bound, moder, fuel): + s1 = openmc.Plane(1.0, 1.0, 0.0, 0.0) + s2 = openmc.Plane(-1.0, 1.0, 0.0, 0.0) + s3 = openmc.Plane(1.0, -1.0, 0.0, 0.0) + s4 = openmc.Plane(-1.0, -1.0, 0.0, 0.0) + surfs = [(s4, 'bottom_left'), + (s1, 'top_right'), + (s2, 'top_left'), + (s3, 'bottom_right')] + univs = [] + for i, (s, name) in enumerate(surfs): + elem_bound = elem_bound.clone() + c1 = openmc.Cell(fill=moder, region=(-s & elem_bound)) + elem_bound = elem_bound.clone() + c2 = openmc.Cell(fill=fuel, region=(+s & elem_bound)) + # universe_id = 6+i + gr_tri = openmc.Universe(name=f'{name}_triangle') + gr_tri.add_cells([c1, c2]) + univs.append(gr_tri) + return univs diff --git a/examples/msbr/_root_geometry.py b/examples/msbr/_root_geometry.py index 2b14dcc88..c5e1bec24 100644 --- a/examples/msbr/_root_geometry.py +++ b/examples/msbr/_root_geometry.py @@ -9,7 +9,7 @@ def shared_root_geometry(): s1 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=208.28, r2=222.71, name='base_octader') s2 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=218.44, r2=215.53, name='smaller_octader') s3 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=228.60, r2=193.97, name='smallest_octader') - zone_i_boundary = -s1 | -s2 | -s3 + zone_i_boundary = (s1, s2, s3) zone_ii_boundary = openmc.ZCylinder(r=256.032, name='iib_boundary') annulus_boundary = openmc.ZCylinder(r=261.112, name='annulus_boundary') @@ -63,36 +63,47 @@ def zoneIIB(zone_i_boundary, zone_ii_boundary, annulus_boundary, core_base, core small_elems_per_octant = 25 elem_cells = [] + zone_iib_reg = None for i, pos in enumerate(large_positions): - r1, r2 = big_radii[i] - t1 = pos - large_half_w - t2 = pos + large_half_w - large_elem = openmc.model.CylinderSector(r1, r2, t1, t2) - elem_hole = hole_region.rotate((0.0, 0.0, pos)) - elem_reg = -large_elem & elem_hole + pos = np.round(pos, 3) + r1_big, r2_big = big_radii[i] + t1_big = pos - large_half_w + t2_big = pos + large_half_w + s1 = openmc.model.CylinderSector(r1_big, r2_big, t1_big, t2_big, name=f'iib_large_element_{pos}') + s2 = openmc.ZCylinder(**hole_args[i]) if isinstance(zone_iib_reg, openmc.Region): - zone_iib_reg = zone_iib_reg | elem_reg + zone_iib_reg = zone_iib_reg & +s1 else: - zone_iib_reg = elem_reg - small_start = t2 + adjacent_angular_offset - r1, r2 = small_radii - t1 = small_start + zone_iib_reg = +s1 + + elem_cells.append(openmc.Cell(fill=moder, region=(-s1 & +s2), name=f'iib_large_element_{pos}')) + elem_cells.append(openmc.Cell(fill=fuel, region=(-s2), + name=f'iib_large_element_fuel_hole_{pos}')) + t1_small = t2_big + adjacent_angular_offset + r1_small, r2_small = small_radii + for i in range(0, small_elems_per_octant): - t2 = t1 + small_angular_width - elem_reg = -openmc.model.CylinderSector(r1, r2, t1, t2) - pos = t2 - (small_angular_width / 2) - zone_iib_reg = zone_iib_reg | elem_reg - t1 = t2 + adjacent_angular_offset - c1 = openmc.Cell(fill=moder, region=(zone_iib_reg & - ~zone_i_boundary & - +core_base & - -core_top), name='iib_moderator') - c2 = openmc.Cell(fill=fuel, region=(~zone_iib_reg & - ~zone_i_boundary & - -zone_ii_boundary & - +core_base & - -core_top), name='iib_fuel') - return c1, c2 + t2_small = t1_small + small_angular_width + + # reflector element + s5 = openmc.model.CylinderSector(r1_small, r2_small, t1_small, t2_small) + pos = t2_small - (small_angular_width / 2) + pos = np.round(pos, 3) + elem_cells.append(openmc.Cell(fill=moder, region=-s5, name=f'iib_small_element_{pos}')) + zone_iib_reg = zone_iib_reg & +s5 + t1_small = t2_small + adjacent_angular_offset + + #universe_id=10 + iib = openmc.Universe(name='zone_iib', cells=elem_cells) + iib.add_cell(openmc.Cell(fill=fuel, region=zone_iib_reg, name='zone_iib_fuel')) + + c1 = openmc.Cell(fill=iib, region=(+zone_i_boundary[0] & + +zone_i_boundary[1] & + +zone_i_boundary[2] & + -zone_ii_boundary & + +core_base & + -core_top), name='zone_iib') + return c1 def annulus(zone_ii_boundary, annulus_boundary, core_base, core_top, fuel): annulus_reg = +zone_ii_boundary & -annulus_boundary & +core_base & -core_top diff --git a/examples/msbr/_root_geometry_optimized.py b/examples/msbr/_root_geometry_optimized.py new file mode 100644 index 000000000..e2f81e2d2 --- /dev/null +++ b/examples/msbr/_root_geometry_optimized.py @@ -0,0 +1,158 @@ +import openmc +import numpy as np + +def shared_root_geometry(): + cr_boundary = openmc.model.rectangular_prism(15.24*2, 15.24*2) + core_base = openmc.ZPlane(z0=0.0, name='core_base') + core_top = openmc.ZPlane(z0=449.58, name='core_top') + + s1 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=208.28, r2=222.71, name='base_octader') + s2 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=218.44, r2=215.53, name='smaller_octader') + s3 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=228.60, r2=193.97, name='smallest_octader') + zone_i_boundary = -s1 | -s2 | -s3 + + zone_ii_boundary = openmc.ZCylinder(r=256.032, name='iib_boundary') + annulus_boundary = openmc.ZCylinder(r=261.112, name='annulus_boundary') + lower_plenum_boundary = openmc.ZPlane(z0=-7.62, name='lower_plenum_boundary') + + zone_bounds = (cr_boundary, zone_i_boundary, zone_ii_boundary) + core_bounds = (annulus_boundary, lower_plenum_boundary, core_base, core_top) + radial_reflector_boundary = openmc.ZCylinder(r=338.328, name='radial_reflector_boundary') + bottom_reflector_boundary = openmc.ZPlane(z0=-76.2, name='bottom_axial_reflector_boundary') + top_reflector_boundary = openmc.ZPlane(z0=525.78, name='top_axial_reflector_boundary') + reflector_bounds = (radial_reflector_boundary, + bottom_reflector_boundary, + top_reflector_boundary) + + radial_vessel_boundary = openmc.ZCylinder(r=343.408, name='radial_vessel_wall', boundary_type='vacuum') + bottom_vessel_boundary = openmc.ZPlane(z0=-81.28, name='bottom_vessel_wall', boundary_type='vacuum') + top_vessel_boundary = openmc.ZPlane(z0=530.86, name='top_vessel_wall', boundary_type='vacuum') + + vessel_bounds = (radial_vessel_boundary, + bottom_vessel_boundary, + top_vessel_boundary) + + return zone_bounds, core_bounds, reflector_bounds, vessel_bounds + +def zoneIIB(zone_i_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top, moder, fuel): + # Large elements + large_angular_width = 3.538 + large_half_w = large_angular_width / 2 + large_positions = np.linspace(0, 315, 8) + r_outer = 256.032 + r_big1 = 229.6 + r_big2 = 223.6 + rb_1 = (r_big1, r_outer) + rb_2 = (r_big2, r_outer) + big_radii = [rb_1, rb_2] * 4 + small_radii = (207.28, r_outer) + + r_hole = 3.0875 + hole_args = ({'x0': 242.679, 'y0': 0.0, 'r': r_hole}, + {'x0': 171.60, 'y0': 171.60, 'r': r_hole}, + {'x0': 0.0, 'y0': 242.679, 'r': r_hole}, + {'x0': -171.60, 'y0': 171.60, 'r': r_hole}, + {'x0': -242.679, 'y0': 0.0, 'r': r_hole}, + {'x0': -171.60, 'y0': -171.60, 'r': r_hole}, + {'x0': 0.0, 'y0': -242.697, 'r': r_hole}, + {'x0': 171.60, 'y0': -171.60, 'r': r_hole}) + + # Small elements + small_angular_width = 0.96 + adjacent_angular_offset = 0.675 #27/40 + small_elems_per_octant = 25 + + elem_cells = [] + for i, pos in enumerate(large_positions): + pos = np.round(pos, 3) + r1_big, r2_big = big_radii[i] + t1_big = pos - large_half_w + t2_big = pos + large_half_w + s1 = openmc.model.CylinderSector(r1_big, r2_big, t1_big, t2_big) + s2 = openmc.ZCylinder(**hole_args[i]) + elem_cells.append(openmc.Cell(fill=moder, region=(-s1 & +s2), name=f'iib_large_element_{pos}')) + elem_cells.append(openmc.Cell(fill=fuel, region=(-s2), + name=f'iib_large_element_fuel_hole_{pos}')) + t1_small = t2_big + adjacent_angular_offset + r1_small, r2_small = small_radii + + # Inter element fuel channel + s3 = openmc.model.CylinderSector(r1_small, r2_small, t2_big, t1_small) + cpos = t2_big + (adjacent_angular_offset / 2) + cpos = np.round(cpos, 3) + elem_cells.append(openmc.Cell(fill=fuel, region=-s3, + name=f'inter_element_fuel_channel_{cpos}')) + + t4a = t1_big - adjacent_angular_offset + s4 = openmc.model.CylinderSector(r1_small, r1_big, t4a, t1_small) + elem_cells.append(openmc.Cell(fill=fuel, region=-s4, + name=f'inter_element_fuel_channel_{pos}')) + + for i in range(0, small_elems_per_octant): + t2_small = t1_small + small_angular_width + + # reflector element + s5 = openmc.model.CylinderSector(r1_small, r2_small, t1_small, t2_small) + pos = t2_small - (small_angular_width / 2) + pos = np.round(pos, 3) + elem_cells.append(openmc.Cell(fill=moder, region=-s5, name=f'iib_small_element_{pos}')) + t1_small = t2_small + adjacent_angular_offset + + # inter-element fuel channel + s6 = openmc.model.CylinderSector(r1_small, r2_small, t2_small, t1_small) + cpos = t2_small + (adjacent_angular_offset/2) + cpos = np.round(cpos, 3) + elem_cells.append(openmc.Cell(fill=fuel, region=-s6, + name=f'inter_element_fuel_channel_{cpos}')) + + #universe_id=10 + iib = openmc.Universe(name='zone_iib') + iib.add_cells(elem_cells) + + c1 = openmc.Cell(fill=iib, region=(~zone_i_boundary & + -zone_ii_boundary & + +core_base & + -core_top), name='zone_iib') + return c1 + +def annulus(zone_ii_boundary, annulus_boundary, core_base, core_top, fuel): + annulus_reg = +zone_ii_boundary & -annulus_boundary & +core_base & -core_top + c1 = openmc.Cell(fill=fuel, region=annulus_reg, name='annulus') + return c1 + +def lower_plenum(core_base, lower_plenum_boundary, annulus_boundary, fuel): + lower_plenum_reg = -core_base & +lower_plenum_boundary & -annulus_boundary + c1 = openmc.Cell(fill=fuel, region=lower_plenum_reg, name='lower_plenum') + return c1 + +def reflectors(annulus_boundary, + radial_reflector_boundary, + lower_plenum_boundary, + bottom_reflector_boundary, + core_top, + top_reflector_boundary, + moder): + radial_reflector_reg = +annulus_boundary & -radial_reflector_boundary & +bottom_reflector_boundary & -top_reflector_boundary + bottom_reflector_reg = -annulus_boundary & -lower_plenum_boundary & +bottom_reflector_boundary + top_reflector_reg = -annulus_boundary & +core_top & -top_reflector_boundary + + c1 = openmc.Cell(fill=moder, region=radial_reflector_reg, name='radial_reflector') + c2 = openmc.Cell(fill=moder, region=bottom_reflector_reg, name='bottom_axial_reflector') + c3 = openmc.Cell(fill=moder, region=top_reflector_reg, name='top_axial_reflector') + return c1, c2, c3 + +def vessel(radial_reflector_boundary, + radial_vessel_boundary, + bottom_vessel_boundary, + top_vessel_boundary, + top_reflector_boundary, + bottom_reflector_boundary, + hast): + radial_vessel_reg = +radial_reflector_boundary & -radial_vessel_boundary & -top_vessel_boundary & +bottom_vessel_boundary + bottom_vessel_reg = -radial_reflector_boundary & -bottom_reflector_boundary & +bottom_vessel_boundary + top_vessel_reg = -radial_reflector_boundary & -top_vessel_boundary & +top_reflector_boundary + + c1 = openmc.Cell(fill=hast, region=radial_vessel_reg, name='radial_vessel_wall') + c2 = openmc.Cell(fill=hast, region=bottom_vessel_reg, name='bottom_vessel_wall') + c3 = openmc.Cell(fill=hast, region=top_vessel_reg, name='top_vessel_wall') + return c1, c2, c3 diff --git a/examples/msbr/model-plotting-optimized.ipynb b/examples/msbr/model-plotting-optimized.ipynb new file mode 100644 index 000000000..018e1b3e8 --- /dev/null +++ b/examples/msbr/model-plotting-optimized.ipynb @@ -0,0 +1,2258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "bdbd4aaa", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56642afc", + "metadata": {}, + "outputs": [], + "source": [ + "# Create materials\n", + "fuel = openmc.Material(name='fuel')\n", + "fuel.set_density('g/cm3', density=3.35)\n", + "fuel.add_components({'Li7': 0.0787474673879085,\n", + " 'Be9': 0.0225566879138321,\n", + " 'F19': 0.454003012179284,\n", + " 'Th232': 0.435579130482336,\n", + " 'U233': 0.00911370203663893},\n", + " percent_type='wo')\n", + "fuel.depletable = True\n", + "fuel.volume = 48710000.0\n", + "\n", + "moder = openmc.Material(name='graphite', temperature=900.0)\n", + "moder.set_density('g/cm3', density=1.84)\n", + "moder.add_nuclide('C0', 1.000, percent_type='wo')\n", + "moder.add_s_alpha_beta('c_Graphite')\n", + "\n", + "hast = openmc.Material(name='hastelloyN', temperature=900.0)\n", + "hast.set_density('g/cm3', density=8.671)\n", + "hast.add_components({'Al27': 0.003,\n", + " 'Ni': 0.677,\n", + " 'W': 0.250,\n", + " 'Cr': 0.070},\n", + " percent_type='wo')\n", + "\n", + "mat = openmc.Materials(materials=[fuel, moder, hast])\n", + "mat.export_to_xml()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bf1ab77", + "metadata": {}, + "outputs": [], + "source": [ + "colormap = {moder: 'purple',\n", + " hast: 'blue',\n", + " fuel: 'yellow'}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e341197", + "metadata": {}, + "outputs": [], + "source": [ + " def _bound_zone_cells(cells_tuples, levels):\n", + " \"\"\"Helper function that moves Zone IA and Zone IIA cells to the\n", + " appropriate height.\"\"\"\n", + " cell_list = []\n", + " n_levels = len(cells_tuples)\n", + " for i, cells in enumerate(cells_tuples):\n", + " if i == 0:\n", + " lower_bound = None\n", + " upper_bound = levels[i]\n", + " elif i == n_levels - 1:\n", + " lower_bound = levels[i-1]\n", + " upper_bound = None\n", + " else:\n", + " lower_bound = levels[i-1]\n", + " upper_bound = levels[i]\n", + " for j, cell in enumerate(cells):\n", + " if lower_bound is None:\n", + " cell.region = cell.region & -upper_bound\n", + " elif upper_bound is None:\n", + " cell.region = cell.region & +lower_bound\n", + " else:\n", + " cell.region = cell.region & +lower_bound & -upper_bound\n", + " cell_list.append(cell)\n", + " return cell_list\n", + "\n", + "def shared_elem_geometry():\n", + " \"\"\"Surfaces and regions shared by Zone IA and Zone IIA elements. \n", + " Specs found in Robertson, 1971 Fig 3.4 (p. 17) and Fig 3.5 (p.18)\"\"\"\n", + "\n", + " elem_bound = openmc.rectangular_prism(5.08*2, 5.08*2) # Pin outer boundary\n", + " eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values())\n", + "\n", + " gr_sq_neg = openmc.rectangular_prism(4.953*2, 4.953*2, corner_radius=0.46) # Graphite square\n", + " (gr_minx, gr_maxx, gr_miny, gr_maxy, \n", + " gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, \n", + " gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, \n", + " gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values())\n", + " \n", + " # slabs that line up with rounded edges\n", + " slab_u = -gr_maxy & +gr_cyl_maxy & -gr_cyl_maxx & +gr_cyl_minx\n", + " slab_l = +gr_minx & -gr_cyl_minx & -gr_cyl_maxy & +gr_cyl_miny\n", + " slab_b = +gr_miny & -gr_cyl_miny & -gr_cyl_maxx & +gr_cyl_minx\n", + " slab_r = -gr_maxx & +gr_cyl_maxx & -gr_cyl_maxy & +gr_cyl_miny\n", + " slabs = (slab_u, slab_l, slab_b, slab_r)\n", + " \n", + " # the rounded edges themselves\n", + " quarter_ul = -gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy\n", + " quarter_br = -gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny\n", + " quarter_lb = -gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny\n", + " quarter_ru = -gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy\n", + " quarters = (quarter_ul, quarter_br, quarter_lb, quarter_ru)\n", + " \n", + " # remaining square\n", + " gr_sq = +gr_cyl_minx & -gr_cyl_maxx & +gr_cyl_miny & -gr_cyl_maxy\n", + "\n", + " # params for main pin section for both I-A and II-A\n", + " r_d = 0.66802\n", + " l1 = 4.28498\n", + " l2 = 4.53898\n", + " l3 = 5.62102\n", + " ul = openmc.ZCylinder(-l1, l2, r_d, name='corner_ul')\n", + " br = openmc.ZCylinder(l1, -l2, r_d, name='corner_br')\n", + " lb = openmc.ZCylinder(-l2, -l1, r_d, name='corner_lb')\n", + " ru = openmc.ZCylinder(l2, l1, r_d, name='corner_ru')\n", + " ul_t = openmc.ZCylinder(-l1, -l3, r_d, name='corner_ul_tip')\n", + " br_t = openmc.ZCylinder(l1, l3, r_d, name='corner_br_tip')\n", + " ru_t = openmc.ZCylinder(-l3, l1, r_d, name='corner_ru_tip')\n", + " lb_t = openmc.ZCylinder(l3, -l1, r_d, name='corner_lb_tip')\n", + " \n", + " gr_ul = -ul & -eb_maxy & +gr_maxy# | \n", + " gr_ul_fill = -ul & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy\n", + " gr_br = -br & +eb_miny & -gr_miny# | \n", + " gr_br_fill = -br & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny \n", + " gr_lb = -lb & +eb_minx & -gr_minx# | \n", + " gr_lb_fill = -lb & +gr_cyl_lb & +gr_minx & -gr_cyl_minx & -gr_cyl_miny \n", + " gr_ru = -ru & -eb_maxx & +gr_maxx# | \n", + " gr_ru_fill = -ru & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxx & +gr_cyl_maxy\n", + "\n", + " gr_ul_t = -ul_t & +eb_miny & -gr_miny\n", + " gr_br_t = -br_t & -eb_maxy & +gr_maxy\n", + " gr_ru_t = -ru_t & +eb_minx & -gr_minx\n", + " gr_lb_t = -lb_t & -eb_maxx & +gr_maxx \n", + " \n", + " gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill)\n", + " gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t)\n", + " \n", + " iec_1 = +gr_cyl_miny & -gr_cyl_maxy & +gr_maxx & -eb_maxx & +lb_t & +ru & +gr_cyl_ru\n", + " iec_2 = (+gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy) & +ru & +br_t & -eb_maxx & -eb_maxy\n", + " iec_3 = +gr_cyl_minx & -gr_cyl_maxx & +gr_maxy & -eb_maxy & +br_t & +ul & +gr_cyl_ul\n", + " iec_4 = (+gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy) & +ul & +ru_t & +eb_minx & -eb_maxy\n", + " iec_5 = +gr_cyl_miny & -gr_cyl_maxy & -gr_minx & +eb_minx & +ru_t & +lb & +gr_cyl_lb\n", + " iec_6 = (+gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny) & +lb & +ul_t & +eb_minx & +eb_miny\n", + " iec_7 = +gr_cyl_minx & -gr_cyl_maxx & -gr_miny & +eb_miny & +ul_t & +br & +gr_cyl_br\n", + " iec_8 = (+gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny) & +br & +lb_t & -eb_maxx & +eb_miny\n", + " \n", + " inter_elem_channel = (iec_1, iec_2, iec_3, iec_4, iec_5, iec_6, iec_7, iec_8)\n", + " #inter_elem_channel = (~gr_sq_neg & elem_bound &\n", + " # +ul & +br & +lb & +ru &\n", + " # +ul_t & +br_t & +ru_t & +lb_t)\n", + " \n", + " \n", + " gr_round_4 = openmc.ZCylinder(r=2.2225, name='gr_round_4')\n", + "\n", + " \n", + " return elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71b4a322", + "metadata": {}, + "outputs": [], + "source": [ + "def zoneIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4):\n", + " \"\"\"Zone IA element. Specs found in Robertson, 1971 Fig 3.4 (p. 17)\"\"\"\n", + " elem_levels = [22.86, 419.10, 438.15, 445.135]\n", + " level_bounds = []\n", + " for level in elem_levels:\n", + " level_bounds.append(openmc.ZPlane(z0=level))\n", + " s1 = openmc.ZCylinder(r=4.953, name='ia_gr_round_1')\n", + " s2 = openmc.ZCylinder(r=1.71069, name='ia_fuel_hole')\n", + " \n", + " h = 12.66\n", + " theta = np.arctan(4.953 / h)\n", + " r2 = (1 / np.cos(theta))**2 - 1\n", + " s3 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_i')\n", + "\n", + " c1 = openmc.Cell(fill=fuel, region=(-s2), name='ia_fuel_inner_1')\n", + " c2 = openmc.Cell(fill=moder, region=(+s2 & -s1), name='ia_moderator_1')\n", + " c3 = openmc.Cell(fill=fuel, region=(+s1 & elem_bound), name='ia_fuel_outer_1')\n", + " ia1 = (c1, c2, c3)\n", + " \n", + " # I-A main (lower 2)\n", + " s2 = s2.clone()\n", + " c4 = c1.clone(clone_materials=False)\n", + " c4.name = 'ia_fuel_inner_main'\n", + " c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq), name='ia_moderator_main')\n", + " #c6 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='ia_fuel_outer_main')\n", + "\n", + " c5_slab_u = openmc.Cell(fill=moder, region=slabs[0], name='ia_moderator_main_slab_u')\n", + " c5_slab_l = openmc.Cell(fill=moder, region=slabs[1], name='ia_moderator_main_slab_l')\n", + " c5_slab_b = openmc.Cell(fill=moder, region=slabs[2], name='ia_moderator_main_slab_b')\n", + " c5_slab_r = openmc.Cell(fill=moder, region=slabs[3], name='ia_moderator_main_slab_r')\n", + " \n", + " c5_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='ia_moderator_main_quarter_u')\n", + " c5_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='ia_moderator_main_quarter_l')\n", + " c5_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='ia_moderator_main_quarter_b')\n", + " c5_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='ia_moderator_main_quarter_r')\n", + " \n", + " c5_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='ia_moderator_main_ul')\n", + " c5_br = openmc.Cell(fill=moder, region=gr_corners[1], name='ia_moderator_main_br')\n", + " c5_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='ia_moderator_main_ru')\n", + " c5_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='ia_moderator_main_lb')\n", + " \n", + " c5_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='ia_moderator_main_ul_fill')\n", + " c5_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='ia_moderator_main_br_fill')\n", + " c5_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='ia_moderator_main_ru_fill')\n", + " c5_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='ia_moderator_main_lb_fill')\n", + " \n", + " c5_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='ia_moderator_main_ul_t')\n", + " c5_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='ia_moderator_main_br_t')\n", + " c5_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='ia_moderator_main_ru_t')\n", + " c5_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='ia_moderator_main_lb_t')\n", + " \n", + " c61 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='ia_fuel_outer_main_1')\n", + " c62 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='ia_fuel_outer_main_2')\n", + " c63 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='ia_fuel_outer_main_3')\n", + " c64 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='ia_fuel_outer_main_4')\n", + " c65 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='ia_fuel_outer_main_5')\n", + " c66 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='ia_fuel_outer_main_6')\n", + " c67 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='ia_fuel_outer_main_7')\n", + " c68 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='ia_fuel_outer_main_8')\n", + "\n", + " iam = (c4, c5,# c6, \n", + " c5_slab_u, c5_slab_l, c5_slab_b, c5_slab_r,\n", + " c5_quarter_u, c5_quarter_l, c5_quarter_b, c5_quarter_r,\n", + " c5_ul, c5_br, c5_ru, c5_lb, \n", + " c5_ulf, c5_brf, c5_ruf, c5_lbf, \n", + " c5_ul_t, c5_br_t, c5_ru_t, c5_lb_t,\n", + " c61, c62, c63, c64, c65, c66, c67, c68)\n", + " \n", + " # I-A 2 (upper 1)\n", + " c7 = c1.clone(clone_materials=False)\n", + " c7.name = 'ia_fuel_inner_2'\n", + " c8 = c2.clone(clone_materials=False)\n", + " c8.name = 'ia_moderator_2'\n", + " c9 = c3.clone(clone_materials=False)\n", + " c9.name = 'ia_fuel_outer_2'\n", + " ia2 = (c7, c8, c9)\n", + "\n", + " # I-A 3 (upper 2)'\n", + " s2 = s2.clone()\n", + " s3 = s3.clone()\n", + " elem_bound = elem_bound.clone()\n", + " c10 = c1.clone(clone_materials=False)\n", + " c10.name = 'ia_fuel_inner_3'\n", + " c11 = openmc.Cell(fill=moder, region=(+s2 & -s3), name='ia_moderator_3')\n", + " c12 = openmc.Cell(fill=fuel, region=(+s3 & elem_bound), name='ia_fuel_outer_3')\n", + " ia3 = (c10, c11, c12) \n", + " \n", + " # I-A 4 (upper 3)\n", + " s2 = s2.clone()\n", + " elem_bound = elem_bound.clone()\n", + " c13 = openmc.Cell(fill=hast, region=(-s2), name='ia_hast')\n", + " c14 = openmc.Cell(fill=moder, region=(+s2 & -gr_round_4), name='ia_moderator_4')\n", + " c15 = openmc.Cell(fill=fuel, region=(+gr_round_4 & elem_bound), name='ia_fuel_outer_4')\n", + " ia4 = (c13, c14, c15)\n", + " \n", + " elem_cells = [ia1, iam, ia2, ia3, ia4]\n", + " # universe_id=1\n", + " ia = openmc.Universe(name='zone_ia')\n", + " ia.add_cells(_bound_zone_cells(elem_cells, level_bounds))\n", + " \n", + " return ia\n", + "\n", + "def zoneIIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4):\n", + " \"\"\"Zone IIA element. Specs found in Robertson, 1971 Fig 3.5 (p. 18)\"\"\"\n", + " elem_levels = [434.34, 436.88, 439.42, 441.96]\n", + " level_bounds = []\n", + " for level in elem_levels:\n", + " level_bounds.append(openmc.ZPlane(z0=level))\n", + " s1 = openmc.ZCylinder(r=3.302, name='iia_fuel_hole_main') # Hole with fuel salt - Fig 3.5, Roberton 1971 (3.27787 - p.47)\n", + " s2 = openmc.ZCylinder(r=0.635, name='iia_fuel_hole_2')\n", + " s3 = openmc.ZCylinder(r=3.65125, name='iia_gr_round_3')\n", + " h = 6.5\n", + " theta = np.arctan(3.65125 / h)\n", + " r2 = (1 / np.cos(theta))**2 - 1\n", + " s4 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_ii')\n", + "\n", + " # II-A main (lower 1)\n", + " c1 = openmc.Cell(fill=fuel, region=(-s1), name='iia_fuel_inner_main')\n", + " c2 = openmc.Cell(fill=moder, region=(+s1 & gr_sq), name='iia_moderator_main')\n", + " \n", + " c2_slab_u = openmc.Cell(fill=moder, region=slabs[0], name='iia_moderator_main_slab_u')\n", + " c2_slab_l = openmc.Cell(fill=moder, region=slabs[1], name='iia_moderator_main_slab_l')\n", + " c2_slab_b = openmc.Cell(fill=moder, region=slabs[2], name='iia_moderator_main_slab_b')\n", + " c2_slab_r = openmc.Cell(fill=moder, region=slabs[3], name='iia_moderator_main_slab_r')\n", + " \n", + " c2_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='iia_moderator_main_quarter_u')\n", + " c2_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='iia_moderator_main_quarter_l')\n", + " c2_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='iia_moderator_main_quarter_b')\n", + " c2_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='iia_moderator_main_quarter_r')\n", + " \n", + " c2_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='iia_moderator_main_ul')\n", + " c2_br = openmc.Cell(fill=moder, region=gr_corners[1], name='iia_moderator_main_br')\n", + " c2_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='iia_moderator_main_ru')\n", + " c2_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='iia_moderator_main_lb')\n", + " \n", + " c2_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='iia_moderator_main_ul_fill')\n", + " c2_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='iia_moderator_main_br_fill')\n", + " c2_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='iia_moderator_main_ru_fill')\n", + " c2_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='iia_moderator_main_lb_fill')\n", + " \n", + " # repeat section tips\n", + " c2_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='iia_moderator_main_ul_t')\n", + " c2_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='iia_moderator_main_br_t')\n", + " c2_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='iia_moderator_main_ru_t')\n", + " c2_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='iia_moderator_main_lb_t')\n", + " \n", + " #c3 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='iia_fuel_outer_main')\n", + " # interelement fuel channel\n", + " c31 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='iia_fuel_outer_main_1')\n", + " c32 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='iia_fuel_outer_main_2')\n", + " c33 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='iia_fuel_outer_main_3')\n", + " c34 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='iia_fuel_outer_main_4')\n", + " c35 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='iia_fuel_outer_main_5')\n", + " c36 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='iia_fuel_outer_main_3')\n", + " c37 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='iia_fuel_outer_main_7')\n", + " c38 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='iia_fuel_outer_main_8')\n", + "\n", + " iiam = (c1, c2, #c3,\n", + " c2_slab_u, c2_slab_l, c2_slab_b, c2_slab_r,\n", + " c2_quarter_u, c2_quarter_l, c2_quarter_b, c2_quarter_r,\n", + " c2_ul, c2_br, c2_ru, c2_lb, \n", + " c2_ulf, c2_brf, c2_ruf, c2_lbf, \n", + " c2_ul_t, c2_br_t, c2_ru_t, c2_lb_t,\n", + " c31, c32, c33, c34, c35, c36, c37, c38)\n", + "\n", + " # II-A 2 (upper 1)\n", + " #gr_sq_neg = gr_sq_neg.clone()\n", + " gr_sq = gr_sq.clone()\n", + " c4 = openmc.Cell(fill=fuel, region=(-s2), name='iia_fuel_inner_2')\n", + " c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq), name='iia_moderator_2')\n", + " c5_slab_u = c2_slab_u.clone(clone_materials=False)\n", + " c5_slab_u.name = 'iia_moderator_2_slab_u'\n", + " c5_slab_l = c2_slab_l.clone(clone_materials=False)\n", + " c5_slab_l.name = 'iia_moderator_2_slab_l'\n", + " c5_slab_b = c2_slab_b.clone(clone_materials=False)\n", + " c5_slab_b.name = 'iia_moderator_2_slab_b'\n", + " c5_slab_r = c2_slab_r.clone(clone_materials=False)\n", + " c5_slab_r.name = 'iia_moderator_2_slab_r'\n", + " \n", + " c5_quarter_u = c2_quarter_u.clone(clone_materials=False)\n", + " c5_quarter_u.name = 'iia_moderator_2_quarter_u'\n", + " c5_quarter_l = c2_quarter_l.clone(clone_materials=False)\n", + " c5_quarter_l.name = 'iia_moderator_2_quarter_l'\n", + " c5_quarter_b = c2_quarter_b.clone(clone_materials=False)\n", + " c5_quarter_b.name = 'iia_moderator_2_quarter_b'\n", + " c5_quarter_r = c2_quarter_r.clone(clone_materials=False)\n", + " c5_quarter_r.name = 'iia_moderator_2_quarter_r'\n", + " \n", + " c5_ul = c2_ul.clone(clone_materials=False)\n", + " c5_ul.name = 'iia_moderator_2_ul'\n", + " c5_br = c2_br.clone(clone_materials=False)\n", + " c5_br.name = 'iia_moderator_2_br'\n", + " c5_ru = c2_ru.clone(clone_materials=False)\n", + " c5_ru.name = 'iia_moderator_2_ru'\n", + " c5_lb = c2_lb.clone(clone_materials=False)\n", + " c5_lb.name = 'iia_moderator_2_lb'\n", + " c5_ulf = c2_ulf.clone(clone_materials=False)\n", + " c5_ulf.name = 'iia_moderator_2_ulf'\n", + " c5_brf = c2_brf.clone(clone_materials=False)\n", + " c5_brf.name = 'iia_moderator_2_brf'\n", + " c5_ruf = c2_ruf.clone(clone_materials=False)\n", + " c5_ruf.name = 'iia_moderator_2_ruf'\n", + " c5_lbf = c2_lbf.clone(clone_materials=False)\n", + " c5_lbf.name = 'iia_moderator_2_lbf'\n", + " c5_ul_t = c2_ul_t.clone(clone_materials=False)\n", + " c5_ul_t.name = 'iia_moderator_2_ul_t'\n", + " c5_br_t = c2_br_t.clone(clone_materials=False)\n", + " c5_br_t.name = 'iia_moderator_2_br_t'\n", + " c5_ru_t = c2_ru_t.clone(clone_materials=False)\n", + " c5_ru_t.name = 'iia_moderator_2_ru_t'\n", + " c5_lb_t = c2_lb_t.clone(clone_materials=False)\n", + " c5_lb_t.name = 'iia_moderator_2_lb_t' \n", + " \n", + " c61 = c31.clone(clone_materials=False)\n", + " c61.name = 'iia_fuel_outer_2_1'\n", + " c62 = c32.clone(clone_materials=False)\n", + " c62.name = 'iia_fuel_outer_2_2'\n", + " c63 = c33.clone(clone_materials=False)\n", + " c63.name = 'iia_fuel_outer_2_3'\n", + " c64 = c34.clone(clone_materials=False)\n", + " c64.name = 'iia_fuel_outer_2_4'\n", + " c65 = c35.clone(clone_materials=False)\n", + " c65.name = 'iia_fuel_outer_2_5'\n", + " c66 = c36.clone(clone_materials=False)\n", + " c66.name = 'iia_fuel_outer_2_6'\n", + " c67 = c37.clone(clone_materials=False)\n", + " c67.name = 'iia_fuel_outer_2_7'\n", + " c68 = c38.clone(clone_materials=False)\n", + " c68.name = 'iia_fuel_outer_2_8'\n", + " iia2 = (c4, c5, #c6,\n", + " c5_slab_u, c5_slab_l, c5_slab_b, c5_slab_r,\n", + " c5_quarter_u, c5_quarter_l, c5_quarter_b, c5_quarter_r,\n", + " c5_ul, c5_br, c5_ru, c5_lb,\n", + " c5_ulf, c5_brf, c5_ruf, c5_lbf, \n", + " c5_ul_t, c5_br_t, c5_ru_t, c5_lb_t,\n", + " c61, c62, c63, c64, c65, c66, c67, c68)\n", + "\n", + " # II-A 3 (upper 2)\n", + " s2 = s2.clone()\n", + " elem_bound = elem_bound.clone()\n", + " c7 = c4.clone(clone_materials=False)\n", + " c7.name = 'iia_fuel_inner_3'\n", + " c8 = openmc.Cell(fill=moder, region=(+s2 & -s3), name='iia_moderator_3')\n", + " c9 = openmc.Cell(fill=fuel, region=(+s3 & elem_bound), name='iia_fuel_outer_3')\n", + " iia3 = (c7, c8, c9)\n", + "\n", + " # II-A 4 (upper 3) \n", + " elem_bound = elem_bound.clone()\n", + " c10 = openmc.Cell(fill=moder, region=(-s4), name='iia_moderator_4')\n", + " c11 = openmc.Cell(fill=fuel, region=(+s4 & elem_bound), name='iia_fuel_outer_4')\n", + " iia4 = (c10, c11)\n", + "\n", + " # II-A 5 (upper 4)\n", + " elem_bound = elem_bound.clone()\n", + " c12 = openmc.Cell(fill=moder, region=(-gr_round_4), name='iia_moderator_5')\n", + " c13 = openmc.Cell(fill=fuel, region=(+gr_round_4 & elem_bound), name='iia_fuel_outer_5')\n", + " iia5 = (c12, c13)\n", + " \n", + " elem_cells = [iiam, iia2, iia3, iia4, iia5]\n", + " # universe_id=2\n", + " iia = openmc.Universe(name='zone_iia')\n", + " iia.add_cells(_bound_zone_cells(elem_cells, level_bounds))\n", + " return iia\n", + "\n", + "def void_cell(elem_bound):\n", + " elem_bound = elem_bound.clone()\n", + " c1 = openmc.Cell(region=elem_bound, name='lattice_void')\n", + " #universe_id=5\n", + " v = openmc.Universe(name='lattice_void')\n", + " v.add_cell(c1)\n", + " return v\n", + "\n", + "def graphite_triangles(elem_bound):\n", + " s1 = openmc.Plane(1.0, 1.0, 0.0, 0.0)\n", + " s2 = openmc.Plane(-1.0, 1.0, 0.0, 0.0)\n", + " s3 = openmc.Plane(1.0, -1.0, 0.0, 0.0)\n", + " s4 = openmc.Plane(-1.0, -1.0, 0.0, 0.0)\n", + " \n", + " surfs = [(s4, 'bottom_left'),\n", + " (s1, 'top_right'),\n", + " (s2, 'top_left'),\n", + " (s3, 'bottom_right')]\n", + " univs = []\n", + " for i, (s, name) in enumerate(surfs):\n", + " elem_bound = elem_bound.clone()\n", + " c1 = openmc.Cell(fill=moder, region=(-s & elem_bound))\n", + " elem_bound = elem_bound.clone()\n", + " c2 = openmc.Cell(fill=fuel, region=(+s & elem_bound))\n", + " \n", + " # universe_id = 6+i\n", + " gr_tri = openmc.Universe(name=f'{name}_triangle')\n", + " gr_tri.add_cells([c1, c2])\n", + " univs.append(gr_tri)\n", + " \n", + " return univs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1dc06889", + "metadata": {}, + "outputs": [], + "source": [ + "elem_bound = openmc.rectangular_prism(5.08*2, 5.08*2) # Pin outer boundary\n", + "eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values())\n", + "\n", + "gr_sq_neg = openmc.rectangular_prism(4.953*2, 4.953*2, corner_radius=0.46) # Graphite square\n", + "(gr_minx, gr_maxx, gr_miny, gr_maxy, \n", + "gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, \n", + "gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, \n", + "gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values())\n", + "\n", + "# slabs that line up with rounded edges\n", + "slab_u = -gr_maxy & +gr_cyl_maxy & -gr_cyl_maxx & +gr_cyl_minx\n", + "slab_l = +gr_minx & -gr_cyl_minx & -gr_cyl_maxy & +gr_cyl_miny\n", + "slab_b = +gr_miny & -gr_cyl_miny & -gr_cyl_maxx & +gr_cyl_minx\n", + "slab_r = -gr_maxx & +gr_cyl_maxx & -gr_cyl_maxy & +gr_cyl_miny\n", + "slabs = (slab_u, slab_l, slab_b, slab_r)\n", + "\n", + "# the rounded edges themselves\n", + "quarter_ul = -gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy\n", + "quarter_br = -gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny\n", + "quarter_lb = -gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny\n", + "quarter_ru = -gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy\n", + "quarters = (quarter_ul, quarter_br, quarter_lb, quarter_ru)\n", + "\n", + "# remaining square\n", + "gr_sq = +gr_cyl_minx & -gr_cyl_maxx & +gr_cyl_miny & -gr_cyl_maxy\n", + "\n", + "# params for main pin section for both I-A and II-A\n", + "r_d = 0.66802\n", + "l1 = 4.28498\n", + "l2 = 4.53898\n", + "l3 = 5.62102\n", + "ul = openmc.ZCylinder(-l1, l2, r_d, name='corner_ul')\n", + "br = openmc.ZCylinder(l1, -l2, r_d, name='corner_br')\n", + "lb = openmc.ZCylinder(-l2, -l1, r_d, name='corner_lb')\n", + "ru = openmc.ZCylinder(l2, l1, r_d, name='corner_ru')\n", + "ul_t = openmc.ZCylinder(-l1, -l3, r_d, name='corner_ul_tip')\n", + "br_t = openmc.ZCylinder(l1, l3, r_d, name='corner_br_tip')\n", + "ru_t = openmc.ZCylinder(-l3, l1, r_d, name='corner_ru_tip')\n", + "lb_t = openmc.ZCylinder(l3, -l1, r_d, name='corner_lb_tip')\n", + "\n", + "gr_ul = -ul & -eb_maxy & +gr_maxy# | \n", + "gr_ul_fill = -ul & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy\n", + "gr_br = -br & +eb_miny & -gr_miny# | \n", + "gr_br_fill = -br & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny \n", + "gr_lb = -lb & +eb_minx & -gr_minx# | \n", + "gr_lb_fill = -lb & +gr_cyl_lb & +gr_minx & -gr_cyl_minx & -gr_cyl_miny \n", + "gr_ru = -ru & -eb_maxx & +gr_maxx# | \n", + "gr_ru_fill = -ru & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxx & +gr_cyl_maxy\n", + "\n", + "gr_ul_t = -ul_t & +eb_miny & -gr_miny\n", + "gr_br_t = -br_t & -eb_maxy & +gr_maxy\n", + "gr_ru_t = -ru_t & +eb_minx & -gr_minx\n", + "gr_lb_t = -lb_t & -eb_maxx & +gr_maxx \n", + "\n", + "gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill)\n", + "gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t)\n", + "\n", + "iec_1 = +gr_cyl_miny & -gr_cyl_maxy & +gr_maxx & -eb_maxx & +lb_t & +ru & +gr_cyl_ru\n", + "iec_2 = (+gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy) & +ru & +br_t & -eb_maxx & -eb_maxy\n", + "iec_3 = +gr_cyl_minx & -gr_cyl_maxx & +gr_maxy & -eb_maxy & +br_t & +ul & +gr_cyl_ul\n", + "iec_4 = (+gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy) & +ul & +ru_t & +eb_minx & -eb_maxy\n", + "iec_5 = +gr_cyl_miny & -gr_cyl_maxy & -gr_minx & +eb_minx & +ru_t & +lb & +gr_cyl_lb\n", + "iec_6 = (+gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny) & +lb & +ul_t & +eb_minx & +eb_miny\n", + "iec_7 = +gr_cyl_minx & -gr_cyl_maxx & -gr_miny & +eb_miny & +ul_t & +br & +gr_cyl_br\n", + "iec_8 = (+gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny) & +br & +lb_t & -eb_maxx & +eb_miny\n", + "\n", + "inter_elem_channel = (iec_1, iec_2, iec_3, iec_4, iec_5, iec_6, iec_7, iec_8)\n", + "#inter_elem_channel = (~gr_sq_neg & elem_bound &\n", + "# +ul & +br & +lb & +ru &\n", + "# +ul_t & +br_t & +ru_t & +lb_t)\n", + "\n", + "\n", + "gr_round_4 = openmc.ZCylinder(r=2.2225, name='gr_round_4')\n", + "\n", + "\"\"\"Zone IA element. Specs found in Robertson, 1971 Fig 3.4 (p. 17)\"\"\"\n", + "elem_levels = [22.86, 419.10, 438.15, 445.135]\n", + "level_bounds = []\n", + "for level in elem_levels:\n", + " level_bounds.append(openmc.ZPlane(z0=level))\n", + "s1 = openmc.ZCylinder(r=4.953, name='ia_gr_round_1')\n", + "s2 = openmc.ZCylinder(r=1.71069, name='ia_fuel_hole')\n", + "\n", + "h = 12.66\n", + "theta = np.arctan(4.953 / h)\n", + "r2 = (1 / np.cos(theta))**2 - 1\n", + "s3 = openmc.ZCone(z0=h + elem_levels[2], r2=r2, name='cone_i')\n", + "\n", + "c1 = openmc.Cell(fill=fuel, region=(-s2), name='ia_fuel_inner_1')\n", + "c2 = openmc.Cell(fill=moder, region=(+s2 & -s1), name='ia_moderator_1')\n", + "c3 = openmc.Cell(fill=fuel, region=(+s1 & elem_bound), name='ia_fuel_outer_1')\n", + "ia1 = (c1, c2, c3)\n", + "\n", + "# I-A main (lower 2)\n", + "s2 = s2.clone()\n", + "c4 = c1.clone(clone_materials=False)\n", + "c4.name = 'ia_fuel_inner_main'\n", + "c5 = openmc.Cell(fill=moder, region=(+s2 & gr_sq), name='ia_moderator_main')\n", + "#c6 = openmc.Cell(fill=fuel, region=(inter_elem_channel), name='ia_fuel_outer_main')\n", + "\n", + "c5_slab_u = openmc.Cell(fill=moder, region=slabs[0], name='ia_moderator_main_slab_u')\n", + "c5_slab_l = openmc.Cell(fill=moder, region=slabs[1], name='ia_moderator_main_slab_l')\n", + "c5_slab_b = openmc.Cell(fill=moder, region=slabs[2], name='ia_moderator_main_slab_b')\n", + "c5_slab_r = openmc.Cell(fill=moder, region=slabs[3], name='ia_moderator_main_slab_r')\n", + "\n", + "c5_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='ia_moderator_main_quarter_u')\n", + "c5_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='ia_moderator_main_quarter_l')\n", + "c5_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='ia_moderator_main_quarter_b')\n", + "c5_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='ia_moderator_main_quarter_r')\n", + "\n", + "c5_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='ia_moderator_main_ul')\n", + "c5_br = openmc.Cell(fill=moder, region=gr_corners[1], name='ia_moderator_main_br')\n", + "c5_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='ia_moderator_main_ru')\n", + "c5_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='ia_moderator_main_lb')\n", + "\n", + "c5_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='ia_moderator_main_ul_fill')\n", + "c5_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='ia_moderator_main_br_fill')\n", + "c5_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='ia_moderator_main_ru_fill')\n", + "c5_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='ia_moderator_main_lb_fill')\n", + "\n", + "c5_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='ia_moderator_main_ul_t')\n", + "c5_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='ia_moderator_main_br_t')\n", + "c5_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='ia_moderator_main_ru_t')\n", + "c5_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='ia_moderator_main_lb_t')\n", + "\n", + "c61 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='ia_fuel_outer_main_1')\n", + "c62 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='ia_fuel_outer_main_2')\n", + "c63 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='ia_fuel_outer_main_3')\n", + "c64 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='ia_fuel_outer_main_4')\n", + "c65 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='ia_fuel_outer_main_5')\n", + "c66 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='ia_fuel_outer_main_6')\n", + "c67 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='ia_fuel_outer_main_7')\n", + "c68 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='ia_fuel_outer_main_8')\n", + "\n", + "iam = (#c4, c5,# c6, \n", + " #c5_slab_u, c5_slab_l, c5_slab_b, c5_slab_r,\n", + " #c5_quarter_u, c5_quarter_l, c5_quarter_b, c5_quarter_r,\n", + " c5_ul, c5_br, c5_ru, c5_lb, \n", + " c5_ulf, c5_brf, c5_ruf, c5_lbf, \n", + " #c5_ul_t, c5_br_t, c5_ru_t, c5_lb_t,\n", + " c61, c62, c63, c64, c65, c66, c67, c68)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62c7dfbb", + "metadata": {}, + "outputs": [], + "source": [ + "u = openmc.Universe()\n", + "u.add_cells(iam)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9681e5c6", + "metadata": {}, + "outputs": [], + "source": [ + "u.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0., 446),\n", + " #color_by='material',\n", + " pixels=(2000,2000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07136110", + "metadata": {}, + "outputs": [], + "source": [ + "elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry()\n", + "ia = zoneIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4)\n", + "elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry()\n", + "iia = zoneIIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4)\n", + "# tres, uno, dos, quatro\n", + "bl, ur, ul, br = graphite_triangles(elem_bound)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "752c5160", + "metadata": {}, + "outputs": [], + "source": [ + "ia.plot(width=(20,20),\n", + " basis='xz',\n", + " colors=colormap,\n", + " origin=(0.,0., 446),\n", + " color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee56bc53", + "metadata": {}, + "outputs": [], + "source": [ + "iia.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,435),\n", + " color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f71ac1b6", + "metadata": {}, + "outputs": [], + "source": [ + "def shared_cr_geometry():\n", + " fuel_hole = openmc.ZCylinder(r=5.08, name='cr_fuel_hole')\n", + " \n", + " elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary\n", + " eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values())\n", + " \n", + " gr_sq_neg = openmc.rectangular_prism(7.23646*2, 7.23646*2, corner_radius=0.99) # Graphite square\n", + " # params for \n", + " (gr_minx, gr_maxx, gr_miny, gr_maxy, \n", + " gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, \n", + " gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, \n", + " gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values())\n", + " \n", + " # slabs that line up with rounded edges\n", + " slab_u = -gr_maxy & +gr_cyl_maxy & -gr_cyl_maxx & +gr_cyl_minx\n", + " slab_l = +gr_minx & -gr_cyl_minx & -gr_cyl_maxy & +gr_cyl_miny\n", + " slab_b = +gr_miny & -gr_cyl_miny & -gr_cyl_maxx & +gr_cyl_minx\n", + " slab_r = -gr_maxx & +gr_cyl_maxx & -gr_cyl_maxy & +gr_cyl_miny\n", + " slabs = (slab_u, slab_l, slab_b, slab_r)\n", + " \n", + " # the rounded edges themselves\n", + " quarter_ul = -gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy\n", + " quarter_br = -gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny\n", + " quarter_lb = -gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny\n", + " quarter_ru = -gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy\n", + " quarters = (quarter_ul, quarter_br, quarter_lb, quarter_ru)\n", + " \n", + " # remaining square\n", + " gr_sq = +gr_cyl_minx & -gr_cyl_maxx & +gr_cyl_miny & -gr_cyl_maxy\n", + " \n", + " r_d = 1.16\n", + " e_d = 2 * r_d / np.sqrt(3)\n", + " r_dt = 0.8\n", + " r_c = 0.18\n", + " l1 = 5.8801\n", + " l2 = 6.505\n", + " l3 = 8.03646\n", + " ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c)\n", + " (ul_u, ul_b, ul_ur, ul_br, ul_bl, ul_ul,\n", + " ul_cyl_bl_in, ul_cyl_bl_out, \n", + " ul_cyl_ul_in, ul_cyl_ul_out,\n", + " ul_cyl_br_in, ul_cyl_br_out,\n", + " ul_cyl_ur_in, ul_cyl_ur_out,\n", + " ul_cyl_l_in, ul_cyl_l_out,\n", + " ul_cyl_r_in, ul_cyl_r_out)= list(ul.get_surfaces().values())\n", + " br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c)\n", + " (br_u, br_b, br_ur, br_br, br_bl, br_ul,\n", + " br_cyl_bl_in, br_cyl_bl_out, \n", + " br_cyl_ul_in, br_cyl_ul_out,\n", + " br_cyl_br_in, br_cyl_br_out,\n", + " br_cyl_ur_in, br_cyl_ur_out,\n", + " br_cyl_l_in, br_cyl_l_out,\n", + " br_cyl_r_in, br_cyl_r_out)= list(br.get_surfaces().values())\n", + " lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + " (lb_r, lb_l, lb_ur, lb_ul, lb_br, lb_bl,\n", + " lb_cyl_lb_in, lb_cyl_lb_out, \n", + " lb_cyl_lu_in, lb_cyl_lu_out,\n", + " lb_cyl_rb_in, lb_cyl_rb_out,\n", + " lb_cyl_ru_in, lb_cyl_ru_out,\n", + " lb_cyl_b_in, lb_cyl_b_out,\n", + " lb_cyl_u_in, lb_cyl_u_out)= list(lb.get_surfaces().values())\n", + " ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + " (ru_r, ru_l, ru_ur, ru_ul, ru_br, ru_bl,\n", + " ru_cyl_lb_in, ru_cyl_lb_out, \n", + " ru_cyl_lu_in, ru_cyl_lu_out,\n", + " ru_cyl_rb_in, ru_cyl_rb_out,\n", + " ru_cyl_ru_in, ru_cyl_ru_out,\n", + " ru_cyl_b_in, ru_cyl_b_out,\n", + " ru_cyl_u_in, ru_cyl_u_out)= list(ru.get_surfaces().values())\n", + " ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip')\n", + " br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip')\n", + " ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip')\n", + " lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip')\n", + " \n", + " gr_ul = -ul_ul & -ul_ur & -eb_maxy & +gr_maxy\n", + " gr_ul_fill = -ul_ul & +ul_bl & +ul_cyl_l_out & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy\n", + " gr_br = +br_br & +br_bl & +eb_miny & -gr_miny\n", + " gr_br_fill = +br_br & -br_ur & +br_cyl_r_out & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny\n", + " gr_lb = +lb_bl & -lb_ul & +eb_minx & -gr_minx\n", + " gr_lb_fill = +lb_bl & +lb_br & +lb_cyl_b_out & +gr_cyl_lb & +gr_minx & -gr_cyl_miny & -gr_cyl_minx\n", + " gr_ru = +ru_br & -ru_ur & -eb_maxx & +gr_maxx\n", + " gr_ru_fill = -ru_ur & -ru_ul & +ru_cyl_u_out & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxy & +gr_cyl_maxx\n", + "\n", + "\n", + " gr_ul_t = -ul_t & +eb_miny & -gr_miny\n", + " gr_br_t = -br_t & -eb_maxy & +gr_maxy\n", + " gr_ru_t = -ru_t & +eb_minx & -gr_minx\n", + " gr_lb_t = -lb_t & -eb_maxx & +gr_maxx \n", + " \n", + " gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill)\n", + " gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t)\n", + " \n", + " # need to decomplexify even #s\n", + " iec_1 = +gr_cyl_miny & +gr_maxx & -eb_maxx & +lb_t & -ru_br\n", + " iec_2 = +gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy & (+ru_ur | +ru_ul | (+ru_cyl_u_in & -ru_cyl_u_out)) & +br_t & -eb_maxx & -eb_maxy\n", + " iec_3 = -gr_cyl_maxx & +gr_maxy & -eb_maxy & +br_t & +ul_ur\n", + " iec_4 = +gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy & (+ul_ul | -ul_bl | (+ul_cyl_l_in & -ul_cyl_l_out)) & +ru_t & +eb_minx & -eb_maxy\n", + " iec_5 = -gr_cyl_maxy & -gr_minx & +eb_minx & +ru_t & +lb_ul\n", + " iec_6 = +gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny & (-lb_bl | -lb_br | (+lb_cyl_b_in & -lb_cyl_b_out)) & +ul_t & +eb_minx & +eb_miny\n", + " iec_7 = +gr_cyl_minx & -gr_miny & +eb_miny & +ul_t & -br_bl\n", + " iec_8 = +gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny & (-br_br | +br_ur | (+br_cyl_r_in & -br_cyl_r_out)) & +lb_t & -eb_maxx & +eb_miny\n", + "\n", + "\n", + " inter_elem_channel = (iec_1, iec_2, iec_3, iec_4, iec_5, iec_6, iec_7, iec_8)\n", + " #inter_elem_channel = (~gr_sq_neg & elem_bound &\n", + " # ~ul & ~br & ~lb & ~ru &\n", + " # +ul_t & +br_t & +ru_t & +lb_t)\n", + " \n", + "\n", + " return fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b41ea3c3", + "metadata": {}, + "outputs": [], + "source": [ + "elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary\n", + "eb_minx, eb_maxx, eb_miny, eb_maxy = list(elem_bound.get_surfaces().values())\n", + "\n", + "gr_sq_neg = openmc.rectangular_prism(7.23646*2, 7.23646*2, corner_radius=0.99) # Graphite\n", + " \n", + "(gr_minx, gr_maxx, gr_miny, gr_maxy, \n", + "gr_cyl_lb, gr_cyl_minx, gr_cyl_miny, \n", + "gr_cyl_ul, gr_cyl_maxy, gr_cyl_br, \n", + "gr_cyl_maxx, gr_cyl_ru) = list(gr_sq_neg.get_surfaces().values())\n", + "\n", + "# the rounded edges themselves\n", + "quarter_ul = -gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy\n", + "quarter_br = -gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny\n", + "quarter_lb = -gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny\n", + "quarter_ru = -gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy\n", + "quarters = (quarter_ul, quarter_br, quarter_lb, quarter_ru)\n", + "\n", + "# remaining square\n", + "gr_sq = +gr_cyl_minx & -gr_cyl_maxx & +gr_cyl_miny & -gr_cyl_maxy\n", + "\n", + "r_d = 1.16\n", + "e_d = 2 * r_d / np.sqrt(3)\n", + "r_dt = 0.8\n", + "r_c = 0.18\n", + "l1 = 5.8801\n", + "l2 = 6.505\n", + "l3 = 8.03646\n", + "ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c)\n", + "(ul_u, ul_b, ul_ur, ul_br, ul_bl, ul_ul,\n", + " ul_cyl_bl_in, ul_cyl_bl_out, \n", + " ul_cyl_ul_in, ul_cyl_ul_out,\n", + " ul_cyl_br_in, ul_cyl_br_out,\n", + " ul_cyl_ur_in, ul_cyl_ur_out,\n", + " ul_cyl_l_in, ul_cyl_l_out,\n", + " ul_cyl_r_in, ul_cyl_r_out)= list(ul.get_surfaces().values())\n", + "br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c)\n", + "(br_u, br_b, br_ur, br_br, br_bl, br_ul,\n", + " br_cyl_bl_in, br_cyl_bl_out, \n", + " br_cyl_ul_in, br_cyl_ul_out,\n", + " br_cyl_br_in, br_cyl_br_out,\n", + " br_cyl_ur_in, br_cyl_ur_out,\n", + " br_cyl_l_in, br_cyl_l_out,\n", + " br_cyl_r_in, br_cyl_r_out)= list(br.get_surfaces().values())\n", + "lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + "(lb_r, lb_l, lb_ur, lb_ul, lb_br, lb_bl,\n", + " lb_cyl_lb_in, lb_cyl_lb_out, \n", + " lb_cyl_lu_in, lb_cyl_lu_out,\n", + " lb_cyl_rb_in, lb_cyl_rb_out,\n", + " lb_cyl_ru_in, lb_cyl_ru_out,\n", + " lb_cyl_b_in, lb_cyl_b_out,\n", + " lb_cyl_u_in, lb_cyl_u_out)= list(lb.get_surfaces().values())\n", + "ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + "(ru_r, ru_l, ru_ur, ru_ul, ru_br, ru_bl,\n", + " ru_cyl_lb_in, ru_cyl_lb_out, \n", + " ru_cyl_lu_in, ru_cyl_lu_out,\n", + " ru_cyl_rb_in, ru_cyl_rb_out,\n", + " ru_cyl_ru_in, ru_cyl_ru_out,\n", + " ru_cyl_b_in, ru_cyl_b_out,\n", + " ru_cyl_u_in, ru_cyl_u_out)= list(ru.get_surfaces().values())\n", + "ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip')\n", + "br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip')\n", + "ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip')\n", + "lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip')\n", + "\n", + "gr_ul = -ul_ul & -ul_ur & -eb_maxy & +gr_maxy\n", + "gr_ul_fill = -ul_ul & +ul_bl & +ul_cyl_l_out & +gr_cyl_ul & -gr_maxy & -gr_cyl_minx & +gr_cyl_maxy\n", + "gr_br = +br_br & +br_bl & +eb_miny & -gr_miny\n", + "gr_br_fill = +br_br & -br_ur & +br_cyl_r_out & +gr_cyl_br & +gr_miny & +gr_cyl_maxx & -gr_cyl_miny\n", + "gr_lb = +lb_bl & -lb_ul & +eb_minx & -gr_minx\n", + "gr_lb_fill = +lb_bl & +lb_br & +lb_cyl_b_out & +gr_cyl_lb & +gr_minx & -gr_cyl_miny & -gr_cyl_minx\n", + "gr_ru = +ru_br & -ru_ur & -eb_maxx & +gr_maxx\n", + "gr_ru_fill = -ru_ur & -ru_ul & +ru_cyl_u_out & +gr_cyl_ru & -gr_maxx & +gr_cyl_maxy & +gr_cyl_maxx\n", + "\n", + "\n", + "gr_ul_t = -ul_t & +eb_miny & -gr_miny\n", + "gr_br_t = -br_t & -eb_maxy & +gr_maxy\n", + "gr_ru_t = -ru_t & +eb_minx & -gr_minx\n", + "gr_lb_t = -lb_t & -eb_maxx & +gr_maxx \n", + "\n", + "gr_corners = (gr_ul, gr_br, gr_lb, gr_ru, gr_ul_fill, gr_br_fill, gr_lb_fill, gr_ru_fill)\n", + "gr_t = (gr_ul_t, gr_br_t, gr_ru_t, gr_lb_t)\n", + "\n", + "# need to decomplexify even #s\n", + "iec_1 = +gr_cyl_miny & +gr_maxx & -eb_maxx & +lb_t & -ru_br\n", + "iec_2 = +gr_cyl_ru & +gr_cyl_maxx & +gr_cyl_maxy & ~ru & +br_t & -eb_maxx & -eb_maxy\n", + "iec_3 = -gr_cyl_maxx & +gr_maxy & -eb_maxy & +br_t & +ul_ur\n", + "iec_4 = +gr_cyl_ul & -gr_cyl_minx & +gr_cyl_maxy & ~ul & +ru_t & +eb_minx & -eb_maxy\n", + "iec_5 = -gr_cyl_maxy & -gr_minx & +eb_minx & +ru_t & +lb_ul\n", + "iec_6 = +gr_cyl_lb & -gr_cyl_minx & -gr_cyl_miny & ~lb & +ul_t & +eb_minx & +eb_miny\n", + "iec_7 = +gr_cyl_minx & -gr_miny & +eb_miny & +ul_t & -br_bl\n", + "iec_8 = +gr_cyl_br & +gr_cyl_maxx & -gr_cyl_miny & ~br & +lb_t & -eb_maxx & +eb_miny\n", + "\n", + "inter_elem_channel = (iec_1, iec_2, iec_3, iec_4, iec_5, iec_6, iec_7, iec_8)\n", + "#inter_elem_channel = (~gr_sq_neg & elem_bound &\n", + "\n", + "c3 = openmc.Cell(fill=moder, region=(gr_sq), name='cr_moderator')\n", + "\n", + "c3_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='cr_moderator_quarter_u')\n", + "c3_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='cr_moderator_quarter_l')\n", + "c3_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='cr_moderator_quarter_b')\n", + "c3_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='cr_moderator_quarter_r')\n", + "\n", + "c3_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='cr_moderator_ul')\n", + "c3_br = openmc.Cell(fill=moder, region=gr_corners[1], name='cr_moderator_br')\n", + "c3_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='cr_moderator_ru')\n", + "c3_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='cr_moderator_lb')\n", + "\n", + "c3_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='cr_moderator_ul_fill')\n", + "c3_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='cr_moderator_br_fill')\n", + "c3_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='cr_moderator_ru_fill')\n", + "c3_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='cr_moderator_lb_fill')\n", + "\n", + "c3_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='cr_moderator_ul_t')\n", + "c3_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='cr_moderator_br_t')\n", + "c3_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='cr_moderator_ru_t')\n", + "c3_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='cr_moderator_lb_t')\n", + "\n", + "#c4 = openmc.Cell(fill=fuel, region=inter_elem_channel, name='cr_fuel_outer')\n", + "c41 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='cr_fuel_outer_1')\n", + "c42 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='cr_fuel_outer_2')\n", + "c43 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='cr_fuel_outer_3')\n", + "c44 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='cr_fuel_outer_4')\n", + "c45 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='cr_fuel_outer_5')\n", + "c46 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='cr_fuel_outer_6')\n", + "c47 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='cr_fuel_outer_7')\n", + "c48 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='cr_fuel_outer_8')\n", + "\n", + "\n", + "#universe_id=3\n", + "cr = openmc.Universe(name='control_rod')\n", + "cr.add_cells([#c3, #c4,\n", + " #c3_quarter_u, c3_quarter_l, c3_quarter_b, c3_quarter_r,\n", + " c3_ul, c3_br, c3_ru, c3_lb, \n", + " #c3_ulf, c3_brf, c3_ruf, c3_lbf, \n", + " c3_ul_t, c3_br_t, c3_ru_t, c3_lb_t,\n", + " c41, c42, c43, c44, c45, c46, c47, c48])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0dbf19f", + "metadata": {}, + "outputs": [], + "source": [ + "cr.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " #color_by='material',\n", + " pixels=(2000,2000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d17047bd", + "metadata": {}, + "outputs": [], + "source": [ + "def control_rod(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel):\n", + " s1 = openmc.ZCylinder(r=4.7625, name='control_rod')\n", + "\n", + " c1 = openmc.Cell(fill=moder, region=-s1, name='control_rod')\n", + " c2 = openmc.Cell(fill=fuel, region=(+s1 & -fuel_hole), name='cr_fuel_inner')\n", + " \n", + " c3 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq), name='cr_moderator')\n", + " \n", + " c3_slab_u = openmc.Cell(fill=moder, region=slabs[0], name='cr_moderator_slab_u')\n", + " c3_slab_l = openmc.Cell(fill=moder, region=slabs[1], name='cr_moderator_slab_l')\n", + " c3_slab_b = openmc.Cell(fill=moder, region=slabs[2], name='cr_moderator_slab_b')\n", + " c3_slab_r = openmc.Cell(fill=moder, region=slabs[3], name='cr_moderator_slab_r')\n", + " \n", + " c3_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='cr_moderator_quarter_u')\n", + " c3_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='cr_moderator_quarter_l')\n", + " c3_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='cr_moderator_quarter_b')\n", + " c3_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='cr_moderator_quarter_r')\n", + " \n", + " c3_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='cr_moderator_ul')\n", + " c3_br = openmc.Cell(fill=moder, region=gr_corners[1], name='cr_moderator_br')\n", + " c3_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='cr_moderator_ru')\n", + " c3_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='cr_moderator_lb')\n", + " \n", + " c3_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='cr_moderator_ul_fill')\n", + " c3_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='cr_moderator_br_fill')\n", + " c3_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='cr_moderator_ru_fill')\n", + " c3_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='cr_moderator_lb_fill')\n", + " \n", + " c3_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='cr_moderator_ul_t')\n", + " c3_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='cr_moderator_br_t')\n", + " c3_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='cr_moderator_ru_t')\n", + " c3_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='cr_moderator_lb_t')\n", + " \n", + " #c4 = openmc.Cell(fill=fuel, region=inter_elem_channel, name='cr_fuel_outer')\n", + " c41 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='cr_fuel_outer_1')\n", + " c42 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='cr_fuel_outer_2')\n", + " c43 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='cr_fuel_outer_3')\n", + " c44 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='cr_fuel_outer_4')\n", + " c45 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='cr_fuel_outer_5')\n", + " c46 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='cr_fuel_outer_6')\n", + " c47 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='cr_fuel_outer_7')\n", + " c48 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='cr_fuel_outer_8')\n", + "\n", + "\n", + " #universe_id=3\n", + " cr = openmc.Universe(name='control_rod')\n", + " cr.add_cells([c1, c2, c3, #c4,\n", + " c3_slab_u, c3_slab_l, c3_slab_b, c3_slab_r,\n", + " c3_quarter_u, c3_quarter_l, c3_quarter_b, c3_quarter_r,\n", + " c3_ul, c3_br, c3_ru, c3_lb, \n", + " c3_ulf, c3_brf, c3_ruf, c3_lbf, \n", + " c3_ul_t, c3_br_t, c3_ru_t, c3_lb_t,\n", + " c41, c42, c43, c44, c45, c46, c47, c48])\n", + " \n", + " return cr\n", + " \n", + "def control_rod_channel(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel):\n", + " c1 = openmc.Cell(fill=fuel, region=(-fuel_hole), name='crc_fuel_inner')\n", + " c2 = openmc.Cell(fill=moder, region=(+fuel_hole & gr_sq), name='crc_moderator')\n", + " \n", + " c2_slab_u = openmc.Cell(fill=moder, region=slabs[0], name='crc_moderator_slab_u')\n", + " c2_slab_l = openmc.Cell(fill=moder, region=slabs[1], name='crc_moderator_slab_l')\n", + " c2_slab_b = openmc.Cell(fill=moder, region=slabs[2], name='crc_moderator_slab_b')\n", + " c2_slab_r = openmc.Cell(fill=moder, region=slabs[3], name='crc_moderator_slab_r')\n", + " \n", + " c2_quarter_u = openmc.Cell(fill=moder, region=quarters[0], name='crc_moderator_quarter_u')\n", + " c2_quarter_l = openmc.Cell(fill=moder, region=quarters[1], name='crc_moderator_quarter_l')\n", + " c2_quarter_b = openmc.Cell(fill=moder, region=quarters[2], name='crc_moderator_quarter_b')\n", + " c2_quarter_r = openmc.Cell(fill=moder, region=quarters[3], name='crc_moderator_quarter_r')\n", + " \n", + " c2_ul = openmc.Cell(fill=moder, region=gr_corners[0], name='crc_moderator_ul')\n", + " c2_br = openmc.Cell(fill=moder, region=gr_corners[1], name='crc_moderator_br')\n", + " c2_ru = openmc.Cell(fill=moder, region=gr_corners[2], name='crc_moderator_ru')\n", + " c2_lb = openmc.Cell(fill=moder, region=gr_corners[3], name='crc_moderator_lb')\n", + " \n", + " c2_ulf = openmc.Cell(fill=moder, region=gr_corners[4], name='crc_moderator_ul_fill')\n", + " c2_brf = openmc.Cell(fill=moder, region=gr_corners[5], name='crc_moderator_br_fill')\n", + " c2_ruf = openmc.Cell(fill=moder, region=gr_corners[6], name='crc_moderator_ru_fill')\n", + " c2_lbf = openmc.Cell(fill=moder, region=gr_corners[7], name='crc_moderator_lb_fill')\n", + " \n", + " c2_ul_t = openmc.Cell(fill=moder, region=gr_t[0], name='crc_moderator_ul_t')\n", + " c2_br_t = openmc.Cell(fill=moder, region=gr_t[1], name='crc_moderator_br_t')\n", + " c2_ru_t = openmc.Cell(fill=moder, region=gr_t[2], name='crc_moderator_ru_t')\n", + " c2_lb_t = openmc.Cell(fill=moder, region=gr_t[3], name='crc_moderator_lb_t')\n", + " \n", + " #c3 = openmc.Cell(fill=fuel, region=inter_elem_channel, name='crc_fuel_outer')\n", + " c31 = openmc.Cell(fill=fuel, region=(inter_elem_channel[0]), name='crc_fuel_outer_1')\n", + " c32 = openmc.Cell(fill=fuel, region=(inter_elem_channel[1]), name='crc_fuel_outer_2')\n", + " c33 = openmc.Cell(fill=fuel, region=(inter_elem_channel[2]), name='crc_fuel_outer_3')\n", + " c34 = openmc.Cell(fill=fuel, region=(inter_elem_channel[3]), name='crc_fuel_outer_4')\n", + " c35 = openmc.Cell(fill=fuel, region=(inter_elem_channel[4]), name='crc_fuel_outer_5')\n", + " c36 = openmc.Cell(fill=fuel, region=(inter_elem_channel[5]), name='crc_fuel_outer_3')\n", + " c37 = openmc.Cell(fill=fuel, region=(inter_elem_channel[6]), name='crc_fuel_outer_7')\n", + " c38 = openmc.Cell(fill=fuel, region=(inter_elem_channel[7]), name='crc_fuel_outer_8')\n", + "\n", + " \n", + "\n", + "\n", + " # universe_id=4\n", + " crc = openmc.Universe(name='control_rod_channel')\n", + " crc.add_cells([c1, c2, #c3, \n", + " c2_slab_u, c2_slab_l, c2_slab_b, c2_slab_r,\n", + " c2_quarter_u, c2_quarter_l, c2_quarter_b, c2_quarter_r,\n", + " c2_ul, c2_br, c2_ru, c2_lb, \n", + " c2_ulf, c2_brf, c2_ruf, c2_lbf, \n", + " c2_ul_t, c2_br_t, c2_ru_t, c2_lb_t,\n", + " c31, c32, c33, c34, c35, c36, c37, c38])\n", + " \n", + " return crc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cea2d476", + "metadata": {}, + "outputs": [], + "source": [ + "fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() \n", + "cr = control_rod(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel)\n", + "fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() \n", + "crc = control_rod_channel(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35f02392", + "metadata": {}, + "outputs": [], + "source": [ + "cr.cells" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb7abdf0", + "metadata": {}, + "outputs": [], + "source": [ + "openmc.Universe(cells=[cr.cells[6724]]).plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " color_by='material',\n", + " pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa83205f", + "metadata": {}, + "outputs": [], + "source": [ + "cr.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " #color_by='material',\n", + " pixels=(2000,2000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37a3fa2a", + "metadata": {}, + "outputs": [], + "source": [ + "crc.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " color_by='material',\n", + " pixels=(400,400))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff53e9bb", + "metadata": {}, + "outputs": [], + "source": [ + "cr.cells" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e965d29", + "metadata": {}, + "outputs": [], + "source": [ + "c1 = cr.cells[129]\n", + "c2 = cr.cells[130]\n", + "c3 = cr.cells[131]\n", + "c4 = cr.cells[132]\n", + "c5 = cr.cells[133]\n", + "c6 = cr.cells[134]\n", + "c7 = cr.cells[135]\n", + "c8 = cr.cells[136]\n", + "ttt = openmc.Universe(cells=[])\n", + "ttt = openmc.Universe(cells=[c1, c2, c3, c4, c5, c6, c7, c8])\n", + "ttt.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd37f04d", + "metadata": {}, + "outputs": [], + "source": [ + "c5 = cr.cells[250]\n", + "c6 = cr.cells[251]\n", + "c7 = cr.cells[252]\n", + "c8 = cr.cells[253]\n", + "ttt = openmc.Universe(cells=[c5, c6, c7, c8])\n", + "ttt.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7721c7b", + "metadata": {}, + "outputs": [], + "source": [ + "c5" + ] + }, + { + "cell_type": "markdown", + "id": "2f5fd590", + "metadata": {}, + "source": [ + "## Figuring out how to decomplexify CR geometry (WIP)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f168514", + "metadata": {}, + "outputs": [], + "source": [ + "r_d = 1.16\n", + "r_dt = 0.8\n", + "e_d = 2 * r_d / np.sqrt(3)\n", + "r_c = 0.18\n", + "l1 = 5.8801\n", + "l2 = 6.505\n", + "l3 = 8.03646\n", + "ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c)\n", + "br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c)\n", + "lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + "ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c)\n", + "gr_sq_neg = openmc.rectangular_prism(7.23646*2, 7.23646*2, corner_radius=0.99) # Graphite square\n", + "elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary\n", + "\n", + "ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip')\n", + "br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip')\n", + "ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip')\n", + "lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip')\n", + "c1 = openmc.Cell(fill=moder, region=ul)\n", + "c2 = openmc.Cell(fill=moder, region=br)\n", + "c3 = openmc.Cell(fill=moder, region=lb)\n", + "c4 = openmc.Cell(fill=moder, region=ru)\n", + "c5 = openmc.Cell(fill=hast, region=gr_sq_neg)\n", + "c6 = openmc.Cell(fill=fuel, region=-ul_t)\n", + "c7 = openmc.Cell(fill=fuel, region=-br_t)\n", + "c8 = openmc.Cell(fill=fuel, region=-ru_t)\n", + "c9 = openmc.Cell(fill=fuel, region=-lb_t)\n", + "c10 = openmc.Cell(fill=moder, region=elem_bound)\n", + "\n", + "tu = openmc.Universe(cells=[c1, c2, c3, c4, c5, c6, c7, c8, c9, c10])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "557bcb8f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67a1317a", + "metadata": {}, + "outputs": [], + "source": [ + "ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c)\n", + "# for x oriented\n", + "(ul_u, ul_b, ul_ur, ul_br, ul_bl, ul_ul,\n", + " ul_cyl_bl_in, ul_cyl_bl_out, \n", + " ul_cyl_ul_in, ul_cyl_ul_out,\n", + " ul_cyl_br_in, ul_cyl_br_out,\n", + " ul_cyl_ur_in, ul_cyl_ur_out,\n", + " ul_cyl_l_in, ul_cyl_l_out,\n", + " ul_cyl_r_in, ul_yl_r_out)= list(ul.get_surfaces().values())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b07e3b36", + "metadata": {}, + "outputs": [], + "source": [ + "c3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "690e90ee", + "metadata": {}, + "outputs": [], + "source": [ + "(lb_r, lb_l, lb_ur, lb_ul, lb_br, lb_bl,\n", + " lb_cyl_lb_in, lb_cyl_lb_out, \n", + " lb_cyl_lu_in, lb_cyl_lu_out,\n", + " lb_cyl_rb_in, lb_cyl_rb_out,\n", + " lb_cyl_ru_in, lb_cyl_ru_out,\n", + " lb_cyl_b_in, lb_cyl_b_out,\n", + " lb_cyl_u_in, lb_yl_u_out)= list(lb.get_surfaces().values())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ea0eaec", + "metadata": {}, + "outputs": [], + "source": [ + "lb.get_surfaces()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f268e09", + "metadata": {}, + "outputs": [], + "source": [ + "#t = openmc.Region.from_expression('(-1448 1449 -1450 1451 1452 -1453 ~((1454 -1460) | (1455 -1461) | (1456 -1462) | (1457 -1463) | (1458 -1464) | (1459 -1465)))', c1.region.get_surfaces())\n", + "#t = openmc.Region.from_expression('(-1448 1449 -1450 1451 1452 -1453 ((-1454 | 1460) (-1455 | 1461) (-1456 | 1462) (-1457 | 1463) (-1458 | 1464) (-1459 | 1465)))', c1.region.get_surfaces())\n", + "t = openmc.Region.from_expression('(-772 773 -774 -775 776 777 (778 -784))', c3.region.get_surfaces())\n", + "\n", + "#t = openmc.Region.from_expression('((-1454 | -1455 | -1456 | -1457 | -1458 | -1459) | 1460 1461 1462 1463 1464 1465)', c1.region.get_surfaces())\n", + "#t = openmc.Region.from_expression('(-1448 1449 -1450 1451 1452 -1453 (-1454 | -1455 | -1456 | -1457 | -1458 | -1459 | 1460 1461 1462 1463 1464 1465))', c1.region.get_surfaces())\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e6ee5dc", + "metadata": {}, + "outputs": [], + "source": [ + "tt = openmc.Cell(fill=moder, region=t)\n", + "ttt = openmc.Universe(cells=[tt])\n", + "ttt.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b22c009", + "metadata": {}, + "outputs": [], + "source": [ + "tu.plot(width=(20,20),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,440),\n", + " #color_by='material',\n", + " pixels=(800,800))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cb184e6", + "metadata": {}, + "outputs": [], + "source": [ + "s1 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=208.28, r2=222.71, name='base_octader')\n", + "s2 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=218.44, r2=215.53, name='smaller_octader')\n", + "s3 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=228.60, r2=193.97, name='smallest_octader')\n", + "oct1 = list((-s1).get_surfaces().values())\n", + "oct2 = list((-s2).get_surfaces().values())\n", + "oct3 = list((-s3).get_surfaces().values())\n", + "octas = (oct1, oct2, oct3)\n", + "oct1_maxy, oct1_miny, oct1_maxx, oct1_minx, oct1_ur, oct1_br, oct1_bl, oct1_ul = octas[0]\n", + "oct2_maxy, oct2_miny, oct2_maxx, oct2_minx, oct2_ur, oct2_br, oct2_bl, oct2_ul = octas[1]\n", + "oct3_maxy, oct3_miny, oct3_maxx, oct3_minx, oct3_ur, oct3_br, oct3_bl, oct3_ul = octas[2]\n", + "\n", + "cr_boundary = openmc.model.rectangular_prism(15.24*2, 15.24*2)\n", + "\n", + "cb_minx, cb_maxx, cb_miny, cb_maxy = list(cr_boundary.get_surfaces().values())\n", + "\n", + "#c1 = openmc.Cell(fill=moder, region=+s3 & -s2)\n", + "c1_ur = openmc.Cell(fill=moder, region=(-oct3_ur & +oct2_ur & -oct2_maxx & -oct2_maxy))\n", + "c1_ul = openmc.Cell(fill=moder, region=(+oct3_ul & -oct2_ul & +oct2_minx & -oct2_maxy))\n", + "c1_bl = openmc.Cell(fill=moder, region=(+oct3_bl & -oct2_bl & +oct2_minx & +oct2_miny))\n", + "c1_br = openmc.Cell(fill=moder, region=(-oct3_br & +oct2_br & -oct2_maxx & +oct2_miny))\n", + "c1 = [c1_ur, c1_ul, c1_bl, c1_br]\n", + "\n", + "#c2 = openmc.Cell(fill=moder, region=-s3)\n", + "c2_r = openmc.Cell(fill=moder, region=(+cb_maxx & +cb_miny & -cb_maxy & -oct3_maxx))\n", + "c2_ur = openmc.Cell(fill=moder, region=(+cb_maxx & +cb_maxy & -oct3_maxx & -oct3_maxy & +oct3_ur))\n", + "c2_u = openmc.Cell(fill=moder, region=(+cb_maxy & +cb_minx & -cb_maxx & -oct3_maxy))\n", + "c2_ul = openmc.Cell(fill=moder, region=(-cb_minx & +cb_maxy & +oct3_minx & -oct3_maxy & -oct3_ul))\n", + "c2_l = openmc.Cell(fill=moder, region=(-cb_minx & +cb_miny & -cb_maxy & +oct3_minx))\n", + "c2_bl = openmc.Cell(fill=moder, region=(-cb_minx & -cb_miny & +oct3_minx & +oct3_miny & -oct3_bl))\n", + "c2_b = openmc.Cell(fill=moder, region=(-cb_miny & +cb_minx & -cb_maxx & +oct3_miny))\n", + "c2_br = openmc.Cell(fill=moder, region=(+cb_maxx & -cb_miny & -oct3_maxx & +oct3_miny & +oct3_br))\n", + "\n", + "\n", + "c2 = [c2_r, c2_ur, c2_u, c2_ul, c2_l, c2_bl, c2_b, c2_br]\n", + "\n", + "#c3 = openmc.Cell(fill=moder, region=-s1 & +s2)\n", + "c3_ur = openmc.Cell(fill=moder, region=(-oct2_ur & +oct1_ur & -oct1_maxx & -oct1_maxy))\n", + "c3_ul = openmc.Cell(fill=moder, region=(+oct2_ul & -oct1_ul & +oct1_minx & -oct1_maxy))\n", + "c3_bl = openmc.Cell(fill=moder, region=(+oct2_bl & -oct1_bl & +oct1_minx & +oct1_miny))\n", + "c3_br = openmc.Cell(fill=moder, region=(-oct2_br & +oct1_br & -oct1_maxx & +oct1_miny))\n", + "c3 = [c3_ur, c3_ul, c3_bl, c3_br]\n", + "\n", + "u = openmc.Universe()\n", + "#u.add_cells(c1)\n", + "#u.add_cells(c2)\n", + "u.add_cells(c3)\n", + "\n", + "u.plot(width=(700,700),\n", + " basis='xy',\n", + " #colors={c1: 'red', c2:'blue', c3:'green'},\n", + " origin=(0.,0.,200),\n", + " #color_by='material',\n", + " pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce941880", + "metadata": {}, + "outputs": [], + "source": [ + "s1 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=208.28, r2=222.71, name='base_octader')\n", + "s2 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=218.44, r2=215.53, name='smaller_octader')\n", + "s3 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=228.60, r2=193.97, name='smallest_octader')\n", + "oct1 = list((-s1).get_surfaces().values())\n", + "oct2 = list((-s2).get_surfaces().values())\n", + "oct3 = list((-s3).get_surfaces().values())\n", + "octas = (oct1, oct2, oct3)\n", + "oct1_maxy, oct1_miny, oct1_maxx, oct1_minx, oct1_ur, oct1_br, oct1_bl, oct1_ul = octas[0]\n", + "oct2_maxy, oct2_miny, oct2_maxx, oct2_minx, oct2_ur, oct2_br, oct2_bl, oct2_ul = octas[1]\n", + "oct3_maxy, oct3_miny, oct3_maxx, oct3_minx, oct3_ur, oct3_br, oct3_bl, oct3_ul = octas[2]\n", + "\n", + "zone_ii_boundary = openmc.ZCylinder(r=256.032, name='iib_boundary')\n", + "\n", + "cr_boundary = openmc.model.rectangular_prism(15.24*2, 15.24*2)\n", + "cb_minx, cb_maxx, cb_miny, cb_maxy = list(cr_boundary.get_surfaces().values())\n", + "\n", + "\n", + "rr = openmc.Cell(fill=moder, region=(+oct3_maxx & +cb_miny & -cb_maxy & -zone_ii_boundary))\n", + "uu = openmc.Cell(fill=moder, region=(+oct3_maxy & +cb_minx & -cb_maxx & -zone_ii_boundary))\n", + "ll = openmc.Cell(fill=moder, region=(-oct3_minx & +cb_miny & -cb_maxy & -zone_ii_boundary))\n", + "bb = openmc.Cell(fill=moder, region=(-oct3_miny & +cb_minx & -cb_maxx & -zone_ii_boundary))\n", + "\n", + "#c1 = openmc.Cell(fill=moder, region=((+s3 & +s2 & +s1) & -zone_ii_boundary))\n", + "ur = openmc.Cell(fill=moder, region=((-oct1_ur | +oct1_maxx | +oct1_maxy) &\n", + " (-oct2_ur | +oct2_maxx | +oct2_maxy) &\n", + " (-oct3_ur | +oct3_maxx | +oct3_maxy) & -zone_ii_boundary\n", + " & +cb_maxx & +cb_maxy))\n", + "ul = openmc.Cell(fill=moder, region=((+oct1_ul | -oct1_minx | +oct1_maxy) &\n", + " (+oct2_ul | -oct2_minx | +oct2_maxy) &\n", + " (+oct3_ul | -oct3_minx | +oct3_maxy) & -zone_ii_boundary\n", + " & -cb_minx & +cb_maxy))\n", + "bl = openmc.Cell(fill=moder, region=((+oct1_bl | -oct1_minx | -oct1_miny) &\n", + " (+oct2_bl | -oct2_minx | -oct2_miny) &\n", + " (+oct3_bl | -oct3_minx | -oct3_miny) & -zone_ii_boundary\n", + " & -cb_minx & -cb_miny))\n", + "br = openmc.Cell(fill=moder, region=((-oct1_br | +oct1_maxx | -oct1_miny) &\n", + " (-oct2_br | +oct2_maxx | -oct2_miny) &\n", + " (-oct3_br | +oct3_maxx | -oct3_miny) & -zone_ii_boundary\n", + " & +cb_maxx & -cb_miny))\n", + "#c3_ul = openmc.Cell(fill=moder, region=(+oct2_ul & -oct1_ul & +oct1_minx & -oct1_maxy))\n", + "#c3_bl = openmc.Cell(fill=moder, region=(+oct2_bl & -oct1_bl & +oct1_minx & +oct1_miny))\n", + "#c3_br = openmc.Cell(fill=moder, region=(-oct2_br & +oct1_br & -oct1_maxx & +oct1_miny))\n", + "#c1_ur = openmc.Cell(fill=moder, region=(-oct3_ur & +oct1_maxx & +oct1_maxy & -zone_ii_boundary))\n", + "#c1 = [c1_ur]\n", + "#c2 = openmc.Cell(fill=moder, region=+s1 & -zone_ii_boundary)\n", + "\n", + "\n", + "#c2 = [c2_r, c2_ur, c2_u, c2_ul, c2_l, c2_bl, c2_b, c2_br]\n", + "\n", + "#c3 = openmc.Cell(fill=moder, region=-s1 & +s2)\n", + "\n", + "u = openmc.Universe()\n", + "#u.add_cells(c1)\n", + "u.add_cells([ur, ul, bl, br, ll, rr, bb, uu])\n", + "#u.add_cells(c3)\n", + "\n", + "u.plot(width=(700,700),\n", + " basis='xy',\n", + " #colors={c1: 'red', c2:'blue', c3:'green'},\n", + " origin=(0.,0.,200),\n", + " #color_by='material',\n", + " pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f439005", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f985e3a", + "metadata": {}, + "outputs": [], + "source": [ + "def shared_root_geometry():\n", + " cr_boundary = openmc.model.rectangular_prism(15.24*2, 15.24*2)\n", + " \n", + " core_base = openmc.ZPlane(z0=0.0, name='core_base')\n", + " core_top = openmc.ZPlane(z0=449.58, name='core_top')\n", + " \n", + " s1 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=208.28, r2=222.71, name='base_octader')\n", + " s2 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=218.44, r2=215.53, name='smaller_octader')\n", + " s3 = openmc.model.IsogonalOctagon(center=(0.0,0.0), r1=228.60, r2=193.97, name='smallest_octader')\n", + " \n", + " oct1 = list((-s1).get_surfaces().values())\n", + " oct2 = list((-s2).get_surfaces().values())\n", + " oct3 = list((-s3).get_surfaces().values())\n", + " zone_i_octas = (oct1, oct2, oct3)\n", + " \n", + " zone_i_boundary = (s1, s2, s3)\n", + "\n", + " zone_ii_boundary = openmc.ZCylinder(r=256.032, name='iib_boundary')\n", + " annulus_boundary = openmc.ZCylinder(r=261.112, name='annulus_boundary')\n", + " lower_plenum_boundary = openmc.ZPlane(z0=-7.62, name='lower_plenum_boundary')\n", + "\n", + " zone_bounds = (cr_boundary, zone_i_boundary, zone_i_octas, zone_ii_boundary)\n", + " \n", + " core_bounds = (annulus_boundary, lower_plenum_boundary, core_base, core_top)\n", + " \n", + " radial_reflector_boundary = openmc.ZCylinder(r=338.328, name='radial_reflector_boundary') \n", + " bottom_reflector_boundary = openmc.ZPlane(z0=-76.2, name='bottom_axial_reflector_boundary')\n", + " top_reflector_boundary = openmc.ZPlane(z0=525.78, name='top_axial_reflector_boundary') \n", + " \n", + " reflector_bounds = (radial_reflector_boundary,\n", + " bottom_reflector_boundary,\n", + " top_reflector_boundary)\n", + " \n", + " radial_vessel_boundary = openmc.ZCylinder(r=343.408, name='radial_vessel_wall', boundary_type='vacuum')\n", + " bottom_vessel_boundary = openmc.ZPlane(z0=-81.28, name='bottom_vessel_wall', boundary_type='vacuum')\n", + " top_vessel_boundary = openmc.ZPlane(z0=530.86, name='top_vessel_wall', boundary_type='vacuum')\n", + " \n", + " vessel_bounds = (radial_vessel_boundary,\n", + " bottom_vessel_boundary,\n", + " top_vessel_boundary)\n", + " \n", + " return zone_bounds, core_bounds, reflector_bounds, vessel_bounds\n", + "\n", + "def cr_lattice(cr_boundary, core_base, core_top):\n", + " fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() \n", + " f = control_rod(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel)\n", + " fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() \n", + " e = control_rod_channel(fuel_hole, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel)\n", + " \n", + " cr = openmc.RectLattice()\n", + " cr.pitch = np.array([15.24, 15.24])\n", + " N = 2 / 2\n", + " cr.lower_left = -1 * cr.pitch * N\n", + " cr.universes = [[f, e],\n", + " [e, f]]\n", + " \n", + " c1 = openmc.Cell(fill=cr, region=(+core_base & -core_top & cr_boundary), name='cr_lattice')\n", + " \n", + " return c1\n", + "\n", + "def main_lattice(zone_i_boundary, zone_i_octas, cr_boundary, core_base, core_top):\n", + " elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry()\n", + " l = zoneIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4)\n", + " elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry()\n", + " z = zoneIIA(elem_bound, gr_sq, slabs, quarters, gr_corners, gr_t, inter_elem_channel, gr_round_4)\n", + " v = void_cell(elem_bound)\n", + " # tres, uno, dos, quatro\n", + " t, u, d, q = graphite_triangles(elem_bound) \n", + " \n", + " main = openmc.RectLattice()\n", + " main.pitch = np.array([10.16, 10.16])\n", + " N = 45 / 2\n", + " main.lower_left = -1 * main.pitch * N\n", + " main.universes = [[v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, d, z, z, z, z, z, z, z, z, z, u, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, v, v, v, d, z, z, z, z, l, l, l, l, l, l, l, l, l, z, z, z, z, u, v, v, v, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, v, d, z, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, z, u, v, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v],\n", + " [v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v],\n", + " [v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v],\n", + " [v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v],\n", + " [v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v],\n", + " [v, v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v, v],\n", + " [v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z],\n", + " [t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v],\n", + " [v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v],\n", + " [v, v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v, v],\n", + " [v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v],\n", + " [v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v],\n", + " [v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v],\n", + " [v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v],\n", + " [v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, v, t, z, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, z, q, v, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, l, l, l, l, l, l, l, l, l, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v],\n", + " [v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, z, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v]]\n", + " \n", + " oct1_maxy, oct1_miny, oct1_maxx, oct1_minx, oct1_ur, oct1_br, oct1_bl, oct1_ul = zone_i_octas[0]\n", + " oct2_maxy, oct2_miny, oct2_maxx, oct2_minx, oct2_ur, oct2_br, oct2_bl, oct2_ul = zone_i_octas[1]\n", + " oct3_maxy, oct3_miny, oct3_maxx, oct3_minx, oct3_ur, oct3_br, oct3_bl, oct3_ul = zone_i_octas[2]\n", + " cb_minx, cb_maxx, cb_miny, cb_maxy = list(cr_boundary.get_surfaces().values())\n", + " \n", + " #c1 = openmc.Cell(fill=moder, region=+s3 & -s2)\n", + " c1_ur = openmc.Cell(fill=main, region=(-oct3_ur & +oct2_ur & -oct2_maxx & -oct2_maxy),\n", + " name='smaller_octader_ur')\n", + " c1_ul = openmc.Cell(fill=main, region=(+oct3_ul & -oct2_ul & +oct2_minx & -oct2_maxy),\n", + " name='smaller_octader_ul')\n", + " c1_bl = openmc.Cell(fill=main, region=(+oct3_bl & -oct2_bl & +oct2_minx & +oct2_miny),\n", + " name='smaller_octader_bl')\n", + " c1_br = openmc.Cell(fill=main, region=(-oct3_br & +oct2_br & -oct2_maxx & +oct2_miny),\n", + " name='smaller_octader_br')\n", + " c1 = [c1_ur, c1_ul, c1_bl, c1_br]\n", + "\n", + " #c2 = openmc.Cell(fill=moder, region=-s3)\n", + " c2_r = openmc.Cell(fill=main, region=(+cb_maxx & +cb_miny & -cb_maxy & -oct3_maxx),\n", + " name='smallest_octader_r')\n", + " c2_ur = openmc.Cell(fill=main, region=(+cb_maxx & +cb_maxy & -oct3_maxx & -oct3_maxy & +oct3_ur),\n", + " name='smallest_octader_ur')\n", + " c2_u = openmc.Cell(fill=main, region=(+cb_maxy & +cb_minx & -cb_maxx & -oct3_maxy),\n", + " name='smallest_octader_u')\n", + " c2_ul = openmc.Cell(fill=main, region=(-cb_minx & +cb_maxy & +oct3_minx & -oct3_maxy & -oct3_ul),\n", + " name='smallest_octader_ul')\n", + " c2_l = openmc.Cell(fill=main, region=(-cb_minx & +cb_miny & -cb_maxy & +oct3_minx),\n", + " name='smallest_octader_l')\n", + " c2_bl = openmc.Cell(fill=main, region=(-cb_minx & -cb_miny & +oct3_minx & +oct3_miny & -oct3_bl),\n", + " name='smallest_octader_bl')\n", + " c2_b = openmc.Cell(fill=main, region=(-cb_miny & +cb_minx & -cb_maxx & +oct3_miny),\n", + " name='smallest_octader_b')\n", + " c2_br = openmc.Cell(fill=main, region=(+cb_maxx & -cb_miny & -oct3_maxx & +oct3_miny & +oct3_br),\n", + " name='smallest_octader_br')\n", + " c2 = [c2_r, c2_ur, c2_u, c2_ul, c2_l, c2_bl, c2_b, c2_br]\n", + "\n", + " #c3 = openmc.Cell(fill=moder, region=-s1 & +s2)\n", + " c3_ur = openmc.Cell(fill=main, region=(-oct2_ur & +oct1_ur & -oct1_maxx & -oct1_maxy),\n", + " name='base_octader_ur')\n", + " c3_ul = openmc.Cell(fill=main, region=(+oct2_ul & -oct1_ul & +oct1_minx & -oct1_maxy),\n", + " name='base_octader_ul')\n", + " c3_bl = openmc.Cell(fill=main, region=(+oct2_bl & -oct1_bl & +oct1_minx & +oct1_miny),\n", + " name='base_octader_bl')\n", + " c3_br = openmc.Cell(fill=main, region=(-oct2_br & +oct1_br & -oct1_maxx & +oct1_miny),\n", + " name='base_octader_br')\n", + " c3 = [c3_ur, c3_ul, c3_bl, c3_br]\n", + " for cells in (c1, c2, c3):\n", + " for cell in cells:\n", + " cell.region = cell.region & +core_base & -core_top \n", + "\n", + " #c1 = openmc.Cell(fill=main, region=(+core_base & -core_top & \n", + " # +zone_i_boundary[2] & -zone_i_boundary[1] & \n", + " # ~cr_boundary), name='main_lattice')\n", + " #c2 = openmc.Cell(fill=main, region=(+core_base & -core_top & \n", + " # -zone_i_boundary[2] & \n", + " # ~cr_boundary), name='main_lattice')\n", + " #c3 = openmc.Cell(fill=main, region=(+core_base & -core_top & \n", + " # -zone_i_boundary[0] & +zone_i_boundary[1] & +zone_i_boundary[2] &\n", + " # ~cr_boundary), name='main_lattice')\n", + " main = c1 + c2 + c3\n", + " #main = [c1] + [c2] + [c3]\n", + " return main\n", + "\n", + "def zoneIIB(zone_i_boundary, zone_i_octas, cr_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top):\n", + " # Large elements\n", + " large_angular_width = 3.538\n", + " large_half_w = large_angular_width / 2\n", + " large_positions = np.linspace(0, 315, 8)\n", + " r_outer = 256.032\n", + " r_big1 = 229.6\n", + " r_big2 = 223.6\n", + " rb_1 = (r_big1, r_outer)\n", + " rb_2 = (r_big2, r_outer)\n", + " big_radii = [rb_1, rb_2] * 4\n", + " small_radii = (207.28, r_outer)\n", + " \n", + " r_hole = 3.0875\n", + " hole_args = ({'x0': 242.679, 'y0': 0.0, 'r': r_hole},\n", + " {'x0': 171.60, 'y0': 171.60, 'r': r_hole},\n", + " {'x0': 0.0, 'y0': 242.679, 'r': r_hole},\n", + " {'x0': -171.60, 'y0': 171.60, 'r': r_hole},\n", + " {'x0': -242.679, 'y0': 0.0, 'r': r_hole},\n", + " {'x0': -171.60, 'y0': -171.60, 'r': r_hole},\n", + " {'x0': 0.0, 'y0': -242.697, 'r': r_hole},\n", + " {'x0': 171.60, 'y0': -171.60, 'r': r_hole})\n", + " \n", + " # Small elements\n", + " small_angular_width = 0.96\n", + " adjacent_angular_offset = 0.675 #27/40\n", + " small_elems_per_octant = 25\n", + " \n", + " elem_cells = []\n", + " for i, pos in enumerate(large_positions):\n", + " pos = np.round(pos, 3)\n", + " r1_big, r2_big = big_radii[i]\n", + " t1_big = pos - large_half_w\n", + " t2_big = pos + large_half_w\n", + " s1 = openmc.model.CylinderSector(r1_big, r2_big, t1_big, t2_big)\n", + " s2 = openmc.ZCylinder(**hole_args[i])\n", + " elem_cells.append(openmc.Cell(fill=moder, region=(-s1 & +s2), name=f'iib_large_element_{pos}'))\n", + " elem_cells.append(openmc.Cell(fill=fuel, region=(-s2),\n", + " name=f'iib_large_element_fuel_hole_{pos}'))\n", + " t1_small = t2_big + adjacent_angular_offset\n", + " r1_small, r2_small = small_radii\n", + " \n", + " # Inter element fuel channel\n", + " s3 = openmc.model.CylinderSector(r1_small, r2_small, t2_big, t1_small)\n", + " cpos = t2_big + (adjacent_angular_offset / 2)\n", + " cpos = np.round(cpos, 3)\n", + " elem_cells.append(openmc.Cell(fill=fuel, region=-s3,\n", + " name=f'inter_element_fuel_channel_{cpos}'))\n", + " \n", + " t4a = t1_big - adjacent_angular_offset\n", + " s4 = openmc.model.CylinderSector(r1_small, r1_big, t4a, t1_small)\n", + " elem_cells.append(openmc.Cell(fill=fuel, region=-s4,\n", + " name=f'inter_element_fuel_channel_{pos}'))\n", + " \n", + " for i in range(0, small_elems_per_octant):\n", + " t2_small = t1_small + small_angular_width\n", + " \n", + " # reflector element\n", + " s5 = openmc.model.CylinderSector(r1_small, r2_small, t1_small, t2_small)\n", + " pos = t2_small - (small_angular_width / 2)\n", + " pos = np.round(pos, 3)\n", + " elem_cells.append(openmc.Cell(fill=moder, region=-s5, name=f'iib_small_element_{pos}'))\n", + " t1_small = t2_small + adjacent_angular_offset\n", + " \n", + " # inter-element fuel channel\n", + " s6 = openmc.model.CylinderSector(r1_small, r2_small, t2_small, t1_small)\n", + " cpos = t2_small + (adjacent_angular_offset/2)\n", + " cpos = np.round(cpos, 3)\n", + " elem_cells.append(openmc.Cell(fill=fuel, region=-s6,\n", + " name=f'inter_element_fuel_channel_{cpos}'))\n", + " \n", + "\n", + " #universe_id=10\n", + " iib = openmc.Universe(name='zone_iib')\n", + " iib.add_cells(elem_cells)\n", + " \n", + " oct1_maxy, oct1_miny, oct1_maxx, oct1_minx, oct1_ur, oct1_br, oct1_bl, oct1_ul = zone_i_octas[0]\n", + " oct2_maxy, oct2_miny, oct2_maxx, oct2_minx, oct2_ur, oct2_br, oct2_bl, oct2_ul = zone_i_octas[1]\n", + " oct3_maxy, oct3_miny, oct3_maxx, oct3_minx, oct3_ur, oct3_br, oct3_bl, oct3_ul = zone_i_octas[2]\n", + " cb_minx, cb_maxx, cb_miny, cb_maxy = list(cr_boundary.get_surfaces().values())\n", + " \n", + " \n", + " c1_r = openmc.Cell(fill=iib, region=(+oct3_maxx & +cb_miny & -cb_maxy & -zone_ii_boundary),\n", + " name='zone_ii_b_r')\n", + " c1_ur = openmc.Cell(fill=iib, region=((-oct1_ur | +oct1_maxx | +oct1_maxy) &\n", + " (-oct2_ur | +oct2_maxx | +oct2_maxy) &\n", + " (-oct3_ur | +oct3_maxx | +oct3_maxy) & -zone_ii_boundary\n", + " & +cb_maxx & +cb_maxy),\n", + " name='zone_ii_b_ur')\n", + " c1_u = openmc.Cell(fill=iib, region=(+oct3_maxy & +cb_minx & -cb_maxx & -zone_ii_boundary),\n", + " name='zone_ii_b_u')\n", + " c1_ul = openmc.Cell(fill=iib, region=((+oct1_ul | -oct1_minx | +oct1_maxy) &\n", + " (+oct2_ul | -oct2_minx | +oct2_maxy) &\n", + " (+oct3_ul | -oct3_minx | +oct3_maxy) & -zone_ii_boundary\n", + " & -cb_minx & +cb_maxy),\n", + " name='zone_ii_b_ul')\n", + " c1_l = openmc.Cell(fill=iib, region=(-oct3_minx & +cb_miny & -cb_maxy & -zone_ii_boundary),\n", + " name='zone_ii_b_l')\n", + " c1_bl = openmc.Cell(fill=iib, region=((+oct1_bl | -oct1_minx | -oct1_miny) &\n", + " (+oct2_bl | -oct2_minx | -oct2_miny) &\n", + " (+oct3_bl | -oct3_minx | -oct3_miny) & -zone_ii_boundary\n", + " & -cb_minx & -cb_miny),\n", + " name='zone_ii_b_bl')\n", + " c1_b = openmc.Cell(fill=iib, region=(-oct3_miny & +cb_minx & -cb_maxx & -zone_ii_boundary),\n", + " name='zone_ii_b_b')\n", + " c1_br = openmc.Cell(fill=iib, region=((-oct1_br | +oct1_maxx | -oct1_miny) &\n", + " (-oct2_br | +oct2_maxx | -oct2_miny) &\n", + " (-oct3_br | +oct3_maxx | -oct3_miny) & -zone_ii_boundary\n", + " & +cb_maxx & -cb_miny),\n", + " name='zone_ii_b_br')\n", + " c1 = openmc.Cell(fill=iib, region=(+zone_i_boundary[0] & +zone_i_boundary[1] & +zone_i_boundary[2] & \n", + " -zone_ii_boundary & \n", + " +core_base & \n", + " -core_top), name='zone_iib')\n", + " #c1 = [c1_r, c1_ur, c1_u, c1_ul, c1_l, c1_bl, c1_b, c1_br]\n", + " #for c in c1:\n", + " # c.region = c.region & +core_base & -core_top\n", + "\n", + " \n", + " return [c1]\n", + "\n", + "def annulus(zone_ii_boundary, annulus_boundary, core_base, core_top):\n", + " annulus_reg = +zone_ii_boundary & -annulus_boundary & +core_base & -core_top\n", + " c1 = openmc.Cell(fill=fuel, region=annulus_reg, name='annulus')\n", + " \n", + " return c1\n", + "\n", + "def lower_plenum(core_base, lower_plenum_boundary, annulus_boundary):\n", + " lower_plenum_reg = -core_base & +lower_plenum_boundary & -annulus_boundary \n", + " c1 = openmc.Cell(fill=fuel, region=lower_plenum_reg, name='lower_plenum')\n", + "\n", + " return c1\n", + "\n", + "def reflectors(annulus_boundary, \n", + " radial_reflector_boundary, \n", + " lower_plenum_boundary,\n", + " bottom_reflector_boundary, \n", + " core_top, \n", + " top_reflector_boundary):\n", + " \n", + " radial_reflector_reg = +annulus_boundary & -radial_reflector_boundary & +bottom_reflector_boundary & -top_reflector_boundary\n", + " bottom_reflector_reg = -annulus_boundary & -lower_plenum_boundary & +bottom_reflector_boundary\n", + " top_reflector_reg = -annulus_boundary & +core_top & -top_reflector_boundary\n", + "\n", + " c1 = openmc.Cell(fill=moder, region=radial_reflector_reg, name='radial_reflector')\n", + " c2 = openmc.Cell(fill=moder, region=bottom_reflector_reg, name='bottom_axial_reflector')\n", + " c3 = openmc.Cell(fill=moder, region=top_reflector_reg, name='top_axial_reflector')\n", + " \n", + " return c1, c2, c3\n", + "\n", + "def vessel(radial_reflector_boundary,\n", + " radial_vessel_boundary,\n", + " bottom_vessel_boundary,\n", + " top_vessel_boundary,\n", + " top_reflector_boundary,\n", + " bottom_reflector_boundary):\n", + " radial_vessel_reg = +radial_reflector_boundary & -radial_vessel_boundary & -top_vessel_boundary & +bottom_vessel_boundary\n", + " bottom_vessel_reg = -radial_reflector_boundary & -bottom_reflector_boundary & +bottom_vessel_boundary\n", + " top_vessel_reg = -radial_reflector_boundary & -top_vessel_boundary & +top_reflector_boundary\n", + "\n", + " \n", + " c1 = openmc.Cell(fill=hast, region=radial_vessel_reg, name='radial_vessel_wall')\n", + " c2 = openmc.Cell(fill=hast, region=bottom_vessel_reg, name='bottom_vessel_wall')\n", + " c3 = openmc.Cell(fill=hast, region=top_vessel_reg, name='top_vessel_wall')\n", + " \n", + " return c1, c2, c3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d440e6d", + "metadata": {}, + "outputs": [], + "source": [ + "t = -zone_i_boundary[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a3a1651", + "metadata": {}, + "outputs": [], + "source": [ + "zone_i_boundary[0].bottom" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93fdcaf0", + "metadata": {}, + "outputs": [], + "source": [ + "maxy, miny, maxx, minx, ur, br, bl, ul = list(t.get_surfaces().values())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdf61881", + "metadata": {}, + "outputs": [], + "source": [ + "zone_bounds, core_bounds, reflector_bounds, vessel_bounds = shared_root_geometry()\n", + "\n", + "cr_boundary, zone_i_boundary, zone_i_octas, zone_ii_boundary = zone_bounds\n", + "annulus_boundary, lower_plenum_boundary, core_base, core_top = core_bounds\n", + "radial_reflector_boundary, bottom_reflector_boundary, top_reflector_boundary = reflector_bounds\n", + "radial_vessel_boundary, bottom_vessel_boundary, top_vessel_boundary= vessel_bounds\n", + "\n", + "#main1, main2, main3 = main_lattice(zone_i_boundary, zone_i_octas, cr_boundary, core_base, core_top)\n", + "main = main_lattice(zone_i_boundary, zone_i_octas, cr_boundary, core_base, core_top)\n", + "cr = cr_lattice(cr_boundary, core_base, core_top)\n", + "iib = zoneIIB(zone_i_boundary, zone_i_octas, cr_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top)\n", + "a = annulus(zone_ii_boundary, annulus_boundary, core_base, core_top)\n", + "lp = lower_plenum(core_base, lower_plenum_boundary, annulus_boundary)\n", + "\n", + "rr, rb, rt = reflectors(annulus_boundary, \n", + " radial_reflector_boundary, \n", + " lower_plenum_boundary,\n", + " bottom_reflector_boundary, \n", + " core_top, \n", + " top_reflector_boundary)\n", + "\n", + "vr, vb, vt = vessel(radial_reflector_boundary,\n", + " radial_vessel_boundary,\n", + " bottom_vessel_boundary,\n", + " top_vessel_boundary,\n", + " top_reflector_boundary,\n", + " bottom_reflector_boundary)\n", + "\n", + "testuniverse = openmc.Universe()\n", + "testuniverse.add_cells([cr, lp, a, rr, rb, rt, vr, vb, vt])\n", + "testuniverse.add_cells(main)\n", + "testuniverse.add_cells(iib)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b888598", + "metadata": {}, + "outputs": [], + "source": [ + "iib.fill.plot(width=(700,700),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,200),\n", + " color_by='material',\n", + " pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a283441", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "testuniverse.plot(width=(700,700),\n", + " basis='xy',\n", + " colors=colormap,\n", + " origin=(0.,0.,200),\n", + " color_by='material',\n", + " pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fab3df4", + "metadata": {}, + "outputs": [], + "source": [ + "#testuniverse.plot(width=(700,650),\n", + "# basis='xz',\n", + "# colors=colormap,\n", + "# origin=(0.,0.,220),\n", + "# color_by='material',\n", + "# pixels=(1000,1000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6b0f4320", + "metadata": {}, + "outputs": [], + "source": [ + "geo = openmc.Geometry()\n", + "geo.root_universe = testuniverse\n", + "geo.remove_redundant_surfaces()\n", + "#geo.export_to_xml()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48d744d2", + "metadata": {}, + "outputs": [], + "source": [ + "geo.export_to_xml()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "524f65b6", + "metadata": {}, + "outputs": [], + "source": [ + "geo.get_all_cells()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec266af0", + "metadata": {}, + "outputs": [], + "source": [ + "plots = openmc.Plots()\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIA-IIA-lower1')\n", + "plot.origin=(215, 0, 10.0)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIA-main')\n", + "plot.origin=(215, 0, 23.0)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIIA-upper1')\n", + "plot.origin=(215, 0, 435)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIA-upper1')\n", + "plot.origin=(215, 0, 420)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIIA-upper2')\n", + "plot.origin=(215, 0, 437)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIA-upper2')\n", + "plot.origin=(215, 0, 439)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIIA-upper3')\n", + "plot.origin=(215, 0, 440)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIA-upper3')\n", + "plot.origin=(215, 0, 448)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='detail-zoneIIA-upper4')\n", + "plot.origin=(215, 0, 442)\n", + "plot.width=(40, 40)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIA-IIA-lower1')\n", + "plot.origin=(0.0, 0, 10.0)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIA-main')\n", + "plot.origin=(0, 0, 23.0)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIIA-upper1')\n", + "plot.origin=(0, 0, 435)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIA-upper1')\n", + "plot.origin=(0, 0, 420)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIIA-upper2')\n", + "plot.origin=(0, 0, 437)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIA-upper2')\n", + "plot.origin=(0, 0, 439)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIIA-upper3')\n", + "plot.origin=(0, 0, 440)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIA-upper3')\n", + "plot.origin=(0, 0, 448)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-zoneIIA-upper4')\n", + "plot.origin=(0, 0, 442)\n", + "plot.width=(600, 600)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xy'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='core-xz-detail-upper')\n", + "plot.origin=(215, 0, 440)\n", + "plot.width=(100, 100)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xz'\n", + "plots.append(plot)\n", + "\n", + "plot = openmc.Plot(name='full-core-xz')\n", + "plot.origin=(0, 0, 200)\n", + "plot.width=(700, 700)\n", + "plot.pixels=(1000,1000)\n", + "plot.color_by='material'\n", + "plot.colors=colormap\n", + "plot.basis='xz'\n", + "plots.append(plot)\n", + "\n", + "plots.export_to_xml()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76ec95ac", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:openmc-env] *", + "language": "python", + "name": "conda-env-openmc-env-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/msbr/openmc_model.py b/examples/msbr/openmc_model.py index f0603b9fd..da3f7a35f 100644 --- a/examples/msbr/openmc_model.py +++ b/examples/msbr/openmc_model.py @@ -34,13 +34,51 @@ mat = openmc.Materials(materials=[fuel, moder, hast]) mat.export_to_xml() -# Geometry +def shared_elem_geometry(elem_type, eb_d, gr_sq_d, gr_sq_r, r_dt, l1, l2, l3, r_es, es_name): + """Surfaces and regions shared by Zone IA and Zone IIA elements. + Specs found in Robertson, 1971 Fig 3.4 (p. 17) and Fig 3.5 (p.18)""" + if elem_type == 'core': + # for main pin section for both I-A and II-A + ul = openmc.ZCylinder(-l1, l2, r_dt, name='corner_ul') + br = openmc.ZCylinder(l1, -l2, r_dt, name='corner_br') + lb = openmc.ZCylinder(-l2, -l1, r_dt, name='corner_lb') + ru = openmc.ZCylinder(l2, l1, r_dt, name='corner_ru') + inter_elem_channel = +ul & +br & +lb & +ru + elif elem_type == 'cr': + # params for control rods + r_d = 1.16 + e_d = 2 * r_d / np.sqrt(3) + r_c = 0.18 + + ul = openmc.model.hexagonal_prism(origin=(-l1, l2), edge_length=e_d, orientation='x', corner_radius=r_c) + br = openmc.model.hexagonal_prism(origin=(l1, -l2), edge_length=e_d, orientation='x',corner_radius=r_c) + lb = openmc.model.hexagonal_prism(origin=(-l2, -l1), edge_length=e_d, orientation='y',corner_radius=r_c) + ru = openmc.model.hexagonal_prism(origin=(l2, l1), edge_length=e_d, orientation='y',corner_radius=r_c) + inter_elem_channel = ~ul & ~br & ~lb & ~ru + + elem_bound = openmc.rectangular_prism(eb_d*2, eb_d*2) # Pin outer boundary + gr_sq_neg = openmc.rectangular_prism(gr_sq_d*2, gr_sq_d*2, corner_radius=gr_sq_r) # Graphite square + + ul_t = openmc.ZCylinder(-l1, -l3, r_dt, name='corner_ul_tip') + br_t = openmc.ZCylinder(l1, l3, r_dt, name='corner_br_tip') + ru_t = openmc.ZCylinder(-l3, l1, r_dt, name='corner_ru_tip') + lb_t = openmc.ZCylinder(l3, -l1, r_dt, name='corner_lb_tip') + + inter_elem_channel = inter_elem_channel & +ul_t & +br_t & +ru_t & +lb_t + + gr_corners = ([ul, 'ul'], [br, 'br'], [ru, 'ru'], [lb, 'lb'], + [ul_t, 'ul_t'], [br_t, 'br_t'], [ru_t, 'ru_t'], [lb_t, 'lb_t']) + + extra_surf = openmc.ZCylinder(r=r_es, name=es_name) + + + return elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, extra_surf + def cr_lattice(cr_boundary, core_base, core_top): - fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel = shared_cr_geometry() - elem_bound = openmc.rectangular_prism(7.62*2, 7.62*2) # Pin outer boundary + elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, fuel_hole = shared_elem_geometry('cr', 7.62, 7.23645, 0.99, 0.8, 5.8801, 6.505, 8.03646, 5.08, 'cr_fuel_hole') - f = control_rod(fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel, moder, fuel) - e = control_rod_channel(fuel_hole, gr_sq_neg, gr_corners, inter_elem_channel, moder, fuel) + f = control_rod(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, fuel_hole, fuel, moder) + e = control_rod_channel(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, fuel_hole, fuel, moder) cr = openmc.RectLattice() cr.pitch = np.array([15.24, 15.24]) @@ -53,8 +91,9 @@ def cr_lattice(cr_boundary, core_base, core_top): return c1 def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): - elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4 = shared_elem_geometry() + elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4 = shared_elem_geometry('core', 5.08, 4.953, 0.46, 0.66802, 4.28498, 4.53898, 5.62102, 2.2225, 'gr_round_4') l = zoneIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, moder, fuel, hast) + elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4 = shared_elem_geometry('core', 5.08, 4.953, 0.46, 0.66802, 4.28498, 4.53898, 5.62102, 2.2225, 'gr_round_4') z = zoneIIA(elem_bound, gr_corners, gr_sq_neg, inter_elem_channel, gr_round_4, moder, fuel) v = void_cell(elem_bound) # tres, uno, dos, quatro @@ -109,8 +148,17 @@ def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): [v, v, v, v, v, v, v, v, v, v, v, t, z, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, z, q, v, v, v, v, v, v, v, v, v, v, v], [v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, l, l, l, l, l, l, l, l, l, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v], [v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, z, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v]] - c1 = openmc.Cell(fill=main, region=(+core_base & -core_top & zone_i_boundary & ~cr_boundary), name='main_lattice') - return c1 + c1 = openmc.Cell(fill=main, region=(+core_base & -core_top & + +zone_i_boundary[2] & -zone_i_boundary[1] & + ~cr_boundary), name='main_lattice_smaller_octader') + c2 = openmc.Cell(fill=main, region=(+core_base & -core_top & + -zone_i_boundary[2] & + ~cr_boundary), name='main_lattice_insite_smallest_octader') + c3 = openmc.Cell(fill=main, region=(+core_base & -core_top & + -zone_i_boundary[0] & +zone_i_boundary[1] & +zone_i_boundary[2] & + ~cr_boundary), name='main_lattice_inside_base_octader_deducted_smaller_smallest') + main = [c1, c2, c3] + return main zone_bounds, core_bounds, reflector_bounds, vessel_bounds = shared_root_geometry() @@ -121,7 +169,7 @@ def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): main = main_lattice(zone_i_boundary, cr_boundary, core_base, core_top) cr = cr_lattice(cr_boundary, core_base, core_top) -iib_m, iib_f = zoneIIB(zone_i_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top, moder, fuel) +iib = zoneIIB(zone_i_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top, moder, fuel) a = annulus(zone_ii_boundary, annulus_boundary, core_base, core_top, fuel) lp = lower_plenum(core_base, lower_plenum_boundary, annulus_boundary, fuel) @@ -142,8 +190,9 @@ def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): hast) geo = openmc.Geometry() -univ = openmc.Universe() -univ.add_cells([cr, main, iib_m, iib_f, lp, a, rr, rb, rt, vr, vb, vt]) +univ = openmc.Universe(cells=[cr, lp, a, rr, rb, rt, vr, vb, vt]) +univ.add_cells(main) +univ.add_cell(iib) geo.root_universe = univ geo.remove_redundant_surfaces() @@ -151,9 +200,9 @@ def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): # Settings settings = openmc.Settings() -settings.particles = 300 -settings.batches = 400 -settings.inactive = 10 +settings.particles = 10000 +settings.batches = 150 +settings.inactive = 25 settings.temperature = {'default': 900, 'method': 'interpolation'} settings.export_to_xml() diff --git a/examples/msbr/openmc_model_optimized.py b/examples/msbr/openmc_model_optimized.py new file mode 100644 index 000000000..fb73f4454 --- /dev/null +++ b/examples/msbr/openmc_model_optimized.py @@ -0,0 +1,471 @@ +import openmc +import numpy as np +from openmc.deplete import CoupledOperator, PredictorIntegrator, CELIIntegrator + +from _core_elements_optimized import * +from _control_rods_optimized import * +from _root_geometry_optimized import * + +# Materials + +T = 900 +fuel = openmc.Material(name='fuel', temperature=T) +fuel.set_density('g/cm3', density=3.35) +fuel.add_components({'Li7': 0.0787474673879085, + 'Be9': 0.0225566879138321, + 'F19': 0.454003012179284, + 'Th232': 0.435579130482336, + 'U233': 0.00911370203663893}, + percent_type='wo') +fuel.depletable = True +fuel.volume = 48710000.0 + +moder = openmc.Material(name='graphite', temperature=T) +moder.set_density('g/cm3', density=1.84) +moder.add_nuclide('C0', 1.000, percent_type='wo') +moder.add_s_alpha_beta('c_Graphite') + +hast = openmc.Material(name='hastelloyN', temperature=T) +hast.set_density('g/cm3', density=8.671) +hast.add_components({'Al27': 0.003, + 'Ni': 0.677, + 'W': 0.250, + 'Cr': 0.070}, + percent_type='wo') + +mat = openmc.Materials(materials=[fuel, moder, hast]) +mat.export_to_xml() + +# Geometry +def cr_lattice(cr_boundary, core_base, core_top): + fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() + f = control_rod(fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, fuel, moder) + fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel = shared_cr_geometry() + e = control_rod_channel(fuel_hole, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, fuel, moder) + + cr = openmc.RectLattice() + cr.pitch = np.array([15.24, 15.24]) + N = 2 / 2 + cr.lower_left = -1 * cr.pitch * N + cr.universes = [[e, f], + [f, e]] + c1 = openmc.Cell(fill=cr, region=(+core_base & -core_top & cr_boundary), name='cr_lattice') + + return c1 + +def main_lattice(zone_i_boundary, cr_boundary, core_base, core_top): + elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry() + l = zoneIA(elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4, moder, fuel, hast) + elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4 = shared_elem_geometry() + z = zoneIIA(elem_bound, gr_sq_neg, gr_corners, gr_t, inter_elem_channel, gr_round_4, moder, fuel) + v = void_cell(elem_bound) + # tres, uno, dos, quatro + t, u, d, q = graphite_triangles(elem_bound, moder, fuel) + + main = openmc.RectLattice() + main.pitch = np.array([10.16, 10.16]) + N = 45 / 2 + main.lower_left = -1 * main.pitch * N + main.universes = [[v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, d, z, z, z, z, z, z, z, z, z, u, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, v, v, v, d, z, z, z, z, l, l, l, l, l, l, l, l, l, z, z, z, z, u, v, v, v, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, v, d, z, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, z, u, v, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v, v], + [v, v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v, v], + [v, v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v, v], + [v, v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v, v], + [v, v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v, v], + [v, v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v, v], + [v, v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v, v], + [v, d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u, v], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [d, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, u], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, v, v, v, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z], + [t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v], + [v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v], + [v, v, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, v, v], + [v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v], + [v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v], + [v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v], + [v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v], + [v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v], + [v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, t, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, q, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, v, t, z, z, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, z, z, q, v, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, l, l, l, l, l, l, l, l, l, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v], + [v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, t, z, z, z, z, z, z, z, z, z, q, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v, v]] + c1 = openmc.Cell(fill=main, region=(+core_base & -core_top & zone_i_boundary & ~cr_boundary), name='main_lattice') + return c1 + +zone_bounds, core_bounds, reflector_bounds, vessel_bounds = shared_root_geometry() + +cr_boundary, zone_i_boundary, zone_ii_boundary = zone_bounds +annulus_boundary, lower_plenum_boundary, core_base, core_top = core_bounds +radial_reflector_boundary, bottom_reflector_boundary, top_reflector_boundary = reflector_bounds +radial_vessel_boundary, bottom_vessel_boundary, top_vessel_boundary= vessel_bounds + +main = main_lattice(zone_i_boundary, cr_boundary, core_base, core_top) +cr = cr_lattice(cr_boundary, core_base, core_top) +iib = zoneIIB(zone_i_boundary, zone_ii_boundary, annulus_boundary, core_base, core_top, moder, fuel) +a = annulus(zone_ii_boundary, annulus_boundary, core_base, core_top, fuel) +lp = lower_plenum(core_base, lower_plenum_boundary, annulus_boundary, fuel) + +rr, rb, rt = reflectors(annulus_boundary, + radial_reflector_boundary, + lower_plenum_boundary, + bottom_reflector_boundary, + core_top, + top_reflector_boundary, + moder) + +vr, vb, vt = vessel(radial_reflector_boundary, + radial_vessel_boundary, + bottom_vessel_boundary, + top_vessel_boundary, + top_reflector_boundary, + bottom_reflector_boundary, + hast) + +geo = openmc.Geometry() +univ = openmc.Universe() +univ.add_cells([cr, main, iib, lp, a, rr, rb, rt, vr, vb, vt]) + +geo.root_universe = univ +geo.remove_redundant_surfaces() +geo.export_to_xml() + +# Settings +settings = openmc.Settings() +settings.particles = 10000 +settings.batches = 150 +settings.inactive = 25 +settings.generations_per_batch = 1 +settings.temperature = {'default': 900, 'method': 'interpolation', 'range': (800, 1000)} +#settings.temperature = {'default': 600, 'method': 'nearest'} +settings.export_to_xml() + +# Plots +plots_3d = False +detail_pixels = (1000, 1000) +full_pixels = (10000, 10000) + +colormap = {moder: 'purple', + hast: 'blue', + fuel: 'yellow'} +## Slice plots +plots = openmc.Plots() + +plot = openmc.Plot(name='serpent-plot-1_full-zoneIA-main') +plot.origin=(0, 0, 150.5) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='serpent-plot-2') +plot.origin=(0, -77.5, 200) +plot.width=(155, 700) +plot.pixels=(1550, 3400) +plot.color_by='material' +plot.colors=colormap +plot.basis='yz' +plots.append(plot) + +plot = openmc.Plot(name='serpent-plot-3') +plot.origin=(0, 0, 155) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='yz' +plots.append(plot) + +plot = openmc.Plot(name='serpent-plot-4') +plot.origin=(16.5, 0, 200) +plot.width=(700, 700) +plot.pixels=(2000, 2000) +plot.color_by='material' +plot.colors=colormap +plot.basis='yz' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIA-IIA-lower1') +plot.origin=(215, 0, 10.0) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIA-main') +plot.origin=(215, 0, 23.0) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIIA-upper1') +plot.origin=(215, 0, 435) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIA-upper1') +plot.origin=(215, 0, 420) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIIA-upper2') +plot.origin=(215, 0, 437) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIA-upper2') +plot.origin=(215, 0, 439) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIIA-upper3') +plot.origin=(215, 0, 440) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIA-upper3') +plot.origin=(215, 0, 448) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='detail-zoneIIA-upper4') +plot.origin=(215, 0, 442) +plot.width=(40, 40) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIA-IIA-lower1') +plot.origin=(0.0, 0, 10.0) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIIA-upper1') +plot.origin=(0, 0, 435) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIA-upper1') +plot.origin=(0, 0, 420) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIIA-upper2') +plot.origin=(0, 0, 437) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIA-upper2') +plot.origin=(0, 0, 439) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIIA-upper3') +plot.origin=(0, 0, 440) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIA-upper3') +plot.origin=(0, 0, 448) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='full-zoneIIA-upper4') +plot.origin=(0, 0, 442) +plot.width=(600, 600) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xy' +plots.append(plot) + +plot = openmc.Plot(name='core-xz-detail-upper') +plot.origin=(215, 0, 440) +plot.width=(100, 100) +plot.pixels=detail_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xz' +plots.append(plot) + +plot = openmc.Plot(name='full-core-xz') +plot.origin=(0, 0, 200) +plot.width=(700, 700) +plot.pixels=full_pixels +plot.color_by='material' +plot.colors=colormap +plot.basis='xz' +plots.append(plot) + +plots.export_to_xml() + +DEPLETE = False +if DEPLETE: + fiss_q = { + "Am243": 212952778.98135212, + "Cm246": 220179493.9541502, + "U239": 196182721.4907133, + "Np239": 198519429.0542043, + "Th233": 185840570.73897627, + "Cf251": 223166396.6983342, + "Am242": 215146730.16368726, + "Cm245": 214624022.18345505, + "Cf254": 230600814.3619568, + "Am241": 211216798.63643932, + "Th232": 197108389.42449385, + "Cm240": 219583368.40646642, + "Th231": 186918512.14598972, + "Bk246": 224446497.874413, + "Cm247": 218956599.9139631, + "U238": 206851381.70909396, + "Bk250": 225432928.78068554, + "U230": 198841127.68309468, + "Cf249": 221434495.10716867, + "U234": 200632850.9958874, + "Cm250": 219425761.1783332, + "Th229": 192235847.44789958, + "Cm241": 219075406.6897822, + "Pu237": 210593272.23024797, + "Am240": 215272544.02927735, + "Cm249": 218622037.52325428, + "Ac226": 183632605.3770991, + "Cf250": 229685291.02082983, + "Th228": 189488754.50737157, + "Cf248": 229015120.40511796, + "Ac227": 183458264.80025893, + "Pu241": 211237715.32232296, + "Pu240": 208612566.66049656, + "Cf252": 230239896.94703457, + "U231": 197643438.24939737, + "Cm242": 212786491.32857716, + "Bk245": 225023484.65451327, + "Np235": 199435370.72904894, + "Pu243": 207499380.63776916, + "Pu239": 208018532.78140113, + "Am242_m1": 215145370.5791048, + "Pu236": 208679081.72160652, + "Bk249": 224740691.06136644, + "Np236": 198952718.20228392, + "Np234": 200175925.99275926, + "U237": 196429642.96756968, + "Cf253": 231148831.53210822, + "U236": 203404311.87546986, + "Es254": 232527659.46555784, + "Ac225": 183891658.531768, + "Cm243": 213375296.0362017, + "Bk248": 224456537.88363716, + "Cm244": 217926766.88448203, + "Pu242": 212072186.50565082, + "U241": 198266755.4887299, + "Np237": 205370480.34853214, + "Th234": 186385345.82281572, + "Pa231": 194099942.4938497, + "Pa230": 194744699.33621296, + "Pa233": 194162901.71835947, + "U240": 207137940.30569986, + "U233": 199796183.56054175, + "Pu246": 208860847.72193536, + "Pa232": 193654730.8348164, + "U232": 193044277.35730234, + "Am244_m1": 213894761.9301219, + "Pu244": 208427244.82356748, + "Np238": 208699370.90691367, + "Bk247": 224883761.19281054, + "Am244": 213894761.9301219, + "Pa229": 194955644.11334947, + "Cm248": 221723145.3723629, + "Fm255": 238051756.2074275, + "Cf246": 229074942.12674516, + "Th230": 188666101.25156796, + "Pu238": 209540012.5125772, + "U235": 202270000.0, + "Th227": 190640950.14927194 + } + model = openmc.Model.from_xml() + op = CoupledOperator(model, chain_file='chain_endfb71_pwr.xml', fission_q=fiss_q) + timesteps = [3] * 12 + integrator = PredictorIntegrator(op, timesteps, timestep_units='d', power=2.25e9) + integrator.integrate() +