diff --git a/core_design/correction_factor.py b/core_design/correction_factor.py index 5713339..a63ad7b 100644 --- a/core_design/correction_factor.py +++ b/core_design/correction_factor.py @@ -206,15 +206,67 @@ def corrected_keff_2d(depletion_2d_results_file, total_height, core_radius=None) f"{estimated_total_leakage_bol_pct:.5f}" if idx == 0 and not np.isnan(estimated_total_leakage_bol_pct) else "" ]) + # plt.figure() + # plt.plot(time_steps, keff_2d_values, marker='o', linestyle='-', color='r', label='keff_2D') + # plt.plot(time_steps, keff_2d_corrected_values, marker='o', linestyle='-', color='g', label='corrected_keff_2D') + # plt.xlabel('Time [days]') + # plt.ylabel('k-effective') + # plt.title('Comparison of keff_2D and corrected_keff_2D vs. Time') + # plt.grid(True) + # plt.legend() + # plt.savefig('keff_comparison_vs_Time.png') + # plt.show() + + + # Plot only the operating-period portion of the depletion history. + # This does not change the depletion calculation or cycle-length result. + plot_limit_days = 2000.0 + + plot_indices = [ + i for i, t in enumerate(time_steps) + if t <= plot_limit_days + ] + + # Include the first point after the limit so the downward trend is visible. + if plot_indices: + last_index = min(plot_indices[-1] + 2, len(time_steps)) + else: + last_index = min(2, len(time_steps)) + plt.figure() - plt.plot(time_steps, keff_2d_values, marker='o', linestyle='-', color='r', label='keff_2D') - plt.plot(time_steps, keff_2d_corrected_values, marker='o', linestyle='-', color='g', label='corrected_keff_2D') + + plt.plot( + time_steps[:last_index], + keff_2d_values[:last_index], + marker='o', + linestyle='-', + color='r', + label='keff_2D' + ) + + plt.plot( + time_steps[:last_index], + keff_2d_corrected_values[:last_index], + marker='o', + linestyle='-', + color='g', + label='corrected_keff_2D' + ) + + plt.axhline( + y=1.0, + color='k', + linestyle='--', + label='k = 1' + ) + plt.xlabel('Time [days]') plt.ylabel('k-effective') plt.title('Comparison of keff_2D and corrected_keff_2D vs. Time') plt.grid(True) plt.legend() - plt.savefig('keff_comparison_vs_Time.png') + plt.tight_layout() + plt.savefig('keff_comparison_vs_Time.png', dpi=300) plt.show() cycle_length = None @@ -246,4 +298,67 @@ def corrected_keff_2d(depletion_2d_results_file, total_height, core_radius=None) estimated_axial_leakage_bol_pct, bol_total_non_leakage_probability, estimated_total_leakage_bol_pct - ) \ No newline at end of file + ) + + + # This function is for steady state calculations + def corrected_keff_steady_state(statepoint_file, total_height, core_radius=None): + geometry = openmc.Geometry.from_xml() + root_universe = geometry.root_universe + + group_edges = np.array([ + 1e-5, 6.7e-2, 3.2e-1, 1, 4, 9.88, + 4.81e1, 4.54e2, 4.9e4, 1.83e5, 8.21e5, 4e7 + ]) + + groups = openmc.mgxs.EnergyGroups(group_edges) + + mgxs_lib = openmc.mgxs.Library(geometry) + mgxs_lib.energy_groups = groups + mgxs_lib.mgxs_types = [ + 'absorption', + 'diffusion-coefficient', + 'transport', + 'scatter matrix', + 'total', + 'scatter' + ] + mgxs_lib.domain_type = 'universe' + mgxs_lib.domains = [root_universe] + mgxs_lib.build_library() + + with openmc.StatePoint(statepoint_file) as sp: + mgxs_lib.load_from_statepoint(sp) + + keff_2d = sp.keff.nominal_value + keff_2d_uncertainty = sp.keff.std_dev + + abs_xs_mg = mgxs_lib.get_mgxs(root_universe, 'absorption') + trans_xs_mg = mgxs_lib.get_mgxs(root_universe, 'transport') + + abs_xs_array = abs_xs_mg.get_xs( + nuclide='total', + mgxs_type='absorption', + collapse=True + ) + + trans_xs_array = trans_xs_mg.get_xs( + nuclide='total', + mgxs_type='transport', + collapse=True + ) + + abs_xs_1g = float(np.mean(abs_xs_array)) + trans_xs_1g = float(np.mean(trans_xs_array)) + + diffcoeff_1g = 1.0 / (3.0 * trans_xs_1g) + diffusion_length_squared = diffcoeff_1g / abs_xs_1g + + extrapolated_height = total_height + (2.0 * diffcoeff_1g) + buckling_axial = (np.pi / extrapolated_height) ** 2 + p_nl_axial = 1.0 / (1.0 + diffusion_length_squared * buckling_axial) + + keff_3d_corrected = p_nl_axial * keff_2d + keff_3d_corrected_uncertainty = p_nl_axial * keff_2d_uncertainty + + return keff_2d, keff_3d_corrected, p_nl_axial \ No newline at end of file diff --git a/core_design/openmc_template_GCMR.py b/core_design/openmc_template_GCMR.py index a22a42e..aa817fd 100644 --- a/core_design/openmc_template_GCMR.py +++ b/core_design/openmc_template_GCMR.py @@ -103,6 +103,90 @@ def create_multiregion_pin_universe(radii, materials, active_core_maxz, active_c cells.append(openmc.Cell(region=+surfs[-1] & -active_core_maxz & +active_core_minz, fill=outer_material)) return openmc.Universe(cells=cells) + + def create_shutdown_pin_universe( + params, + rod_radius, + rod_name, + active_core_maxz, + active_core_minz, + absorber_material, + coolant_material, + outer_material + ): + """ + Simplified movable shutdown channel. + + ARO: + Rod withdrawn; channel contains helium. + + ARI: + Enriched B4C rod inserted. + """ + + if rod_radius <= 0.0: + raise ValueError( + f"{rod_name} radius must be greater than zero." + ) + + max_radius = 0.5 * params['Lattice Pitch'] + + if rod_radius >= max_radius: + raise ValueError( + f"{rod_name} radius ({rod_radius} cm) must be smaller " + f"than half the lattice pitch ({max_radius} cm)." + ) + + rod_surface = openmc.ZCylinder( + r=rod_radius, + name=f'{rod_name}_surface' + ) + + inside_region = ( + -rod_surface + & -active_core_maxz + & +active_core_minz + ) + + outside_region = ( + +rod_surface + & -active_core_maxz + & +active_core_minz + ) + + if params['Shutdown Margin Calc']: + print( + f">>> {rod_name}: INSERTED B4C, " + f"radius = {rod_radius} cm" + ) + inner_material = absorber_material + inner_name = f'{rod_name}_inserted' + else: + print( + f">>> {rod_name}: WITHDRAWN Helium, " + f"radius = {rod_radius} cm" + ) + inner_material = coolant_material + inner_name = f'{rod_name}_withdrawn' + + inner_cell = openmc.Cell( + name=inner_name, + fill=inner_material, + region=inside_region + ) + + outer_cell = openmc.Cell( + name=f'{rod_name}_outer_graphite', + fill=outer_material, + region=outside_region + ) + + return openmc.Universe( + name=f'{rod_name}_universe', + cells=[inner_cell, outer_cell] + ) + + def create_assembly(num_rings, lattice_pitch, inner_fill, fuel_pin , moderator_pin, outer_ring=None, simplified_output=True): # Create a hexagonal lattice for the assembly assembly = openmc.HexLattice() @@ -120,10 +204,45 @@ def create_assembly(num_rings, lattice_pitch, inner_fill, fuel_pin , moderator_p # Initialize the count of fuel cells fuel_cells = 1 # Loop to create the rings of fuel pins around the center - for n in range(1, num_rings-1): - ring_cells = 6*n - rings.insert(0, [fuel_pin]*ring_cells) - fuel_cells += ring_cells + # for n in range(1, num_rings-1): + # ring_cells = 6*n + # rings.insert(0, [fuel_pin]*ring_cells) + # fuel_cells += ring_cells + + for n in range(1, num_rings - 1): + ring_cells = 6 * n + ring = [fuel_pin] * ring_cells + + if ( + shutdown_pin is not None + and shutdown_ring is not None + and n == shutdown_ring + and number_of_shutdown_rods > 0 + ): + if ring_cells % number_of_shutdown_rods != 0: + raise ValueError( + f"Cannot distribute {number_of_shutdown_rods} shutdown rods " + f"uniformly over assembly ring {n}, which contains " + f"{ring_cells} positions." + ) + + spacing = ring_cells // number_of_shutdown_rods + + for rod_number in range(number_of_shutdown_rods): + index = rod_number * spacing + ring[index] = shutdown_pin + + fuel_cells += ring_cells - number_of_shutdown_rods + + print( + f">>> Added {number_of_shutdown_rods} shutdown positions " + f"to assembly ring {n}" + ) + + else: + fuel_cells += ring_cells + + rings.insert(0, ring) if outer_ring: rings.insert(0, outer_ring) @@ -250,6 +369,52 @@ def create_drums_universe_CGMR(params, absorber_thickness, drum_radius, params['Hex Lattice Radius'] = params['Lattice Pitch'] /np.sqrt(3) # Define the boundary of the hexagonal prism with the given edge length hex_boundary = openmc.model.hexagonal_prism(edge_length= params['Hex Lattice Radius']) + + # Central-assembly shutdown pin + central_shutdown_pin_universe = create_shutdown_pin_universe( + params=params, + rod_radius=params['Central Shutdown Rod Radius'], + rod_name='central_shutdown_rod', + active_core_maxz=active_core_maxz, + active_core_minz=active_core_minz, + absorber_material=control_drum_absorber, + coolant_material=coolant, + outer_material=materials_database[params['Moderator']] + ) + + central_shutdown_lattice_hex = openmc.Universe( + name='central_shutdown_lattice_hex', + cells=[ + openmc.Cell( + name='central_shutdown_lattice_cell', + fill=central_shutdown_pin_universe, + region=hex_boundary + ) + ] + ) + + # Surrounding-assembly shutdown pin + surrounding_shutdown_pin_universe = create_shutdown_pin_universe( + params=params, + rod_radius=params['Surrounding Shutdown Rod Radius'], + rod_name='surrounding_shutdown_rod', + active_core_maxz=active_core_maxz, + active_core_minz=active_core_minz, + absorber_material=control_drum_absorber, + coolant_material=coolant, + outer_material=materials_database[params['Moderator']] + ) + + surrounding_shutdown_lattice_hex = openmc.Universe( + name='surrounding_shutdown_lattice_hex', + cells=[ + openmc.Cell( + name='surrounding_shutdown_lattice_cell', + fill=surrounding_shutdown_pin_universe, + region=hex_boundary + ) + ] + ) # Create an instance of HexLattice fuel_lattice = openmc.HexLattice() # Set the center of the hexagonal lattice @@ -281,9 +446,75 @@ def create_drums_universe_CGMR(params, absorber_thickness, drum_radius, # Sec. 3 : Fuel ASSEMBLY # ************************************************************************************************************************** - assembly_universe, assembly_fuel_cells = create_assembly(params['Assembly Rings'] , params['Lattice Pitch'],\ - openmc.Universe(cells=[openmc.Cell(fill= materials_database[params['Moderator']])]),\ - fuel_lattice_hex, booster_lattice_hex, outer_ring=None, simplified_output=False) + # assembly_universe, assembly_fuel_cells = create_assembly(params['Assembly Rings'] , params['Lattice Pitch'],\ + # openmc.Universe(cells=[openmc.Cell(fill= materials_database[params['Moderator']])]),\ + # fuel_lattice_hex, booster_lattice_hex, outer_ring=None, simplified_output=False) + + assembly_inner_fill = openmc.Universe( + cells=[ + openmc.Cell( + fill=materials_database[params['Moderator']] + ) + ] + ) + + assembly_universe, assembly_fuel_cells = create_assembly( + num_rings=params['Assembly Rings'], + lattice_pitch=params['Lattice Pitch'], + inner_fill=assembly_inner_fill, + fuel_pin=fuel_lattice_hex, + moderator_pin=booster_lattice_hex, + outer_ring=None, + simplified_output=False + ) + + central_shutdown_assembly_universe, central_shutdown_fuel_cells = ( + create_assembly( + num_rings=params['Assembly Rings'], + lattice_pitch=params['Lattice Pitch'], + inner_fill=assembly_inner_fill, + fuel_pin=fuel_lattice_hex, + moderator_pin=booster_lattice_hex, + shutdown_pin=central_shutdown_lattice_hex, + shutdown_ring=params['Central Shutdown Rod Ring'], + number_of_shutdown_rods=params[ + 'Central Shutdown Rod Count' + ], + outer_ring=None, + simplified_output=False + ) + ) + + surrounding_shutdown_assembly_universe, surrounding_shutdown_fuel_cells = ( + create_assembly( + num_rings=params['Assembly Rings'], + lattice_pitch=params['Lattice Pitch'], + inner_fill=assembly_inner_fill, + fuel_pin=fuel_lattice_hex, + moderator_pin=booster_lattice_hex, + shutdown_pin=surrounding_shutdown_lattice_hex, + shutdown_ring=params['Surrounding Shutdown Rod Ring'], + number_of_shutdown_rods=params[ + 'Surrounding Shutdown Rod Count' + ], + outer_ring=None, + simplified_output=False + ) + ) + + print( + "Normal assembly fuel positions:", + assembly_fuel_cells + ) + print( + "Central shutdown assembly fuel positions:", + central_shutdown_fuel_cells + ) + print( + "Surrounding shutdown assembly fuel positions:", + surrounding_shutdown_fuel_cells + ) + if params['plotting'] == "Y": # plotting @@ -370,11 +601,35 @@ def create_drums_universe_CGMR(params, absorber_thickness, drum_radius, active_core.pitch = (params['Assembly FTF'],) active_core.outer = openmc.Universe(cells=[openmc.Cell(fill= materials_database[params['Radial Reflector']])]) # reflector Area - rings = [[assembly_universe]] + # rings = [[assembly_universe]] + # assembly_number = 1 + # for n in range(1, params['Core Rings']-1): + # ring_cells = 6*n + # rings.insert(0, [assembly_universe]*ring_cells) + # assembly_number += ring_cells + + # Center: one assembly with 12 large shutdown rods + rings = [[central_shutdown_assembly_universe]] + assembly_number = 1 - for n in range(1, params['Core Rings']-1): - ring_cells = 6*n - rings.insert(0, [assembly_universe]*ring_cells) + + for n in range(1, params['Core Rings'] - 1): + ring_cells = 6 * n + + if n == 1: + # First core ring: six assemblies, each containing + # six smaller shutdown rods. + rings.insert( + 0, + [surrounding_shutdown_assembly_universe] * ring_cells + ) + else: + # All remaining inner-core assemblies are normal. + rings.insert( + 0, + [assembly_universe] * ring_cells + ) + assembly_number += ring_cells rings.insert(0, flatten_list([[ca] + [ea]*( params['Core Rings']-2)\ @@ -388,23 +643,63 @@ def create_drums_universe_CGMR(params, absorber_thickness, drum_radius, active_core_cell = openmc.Cell(fill=active_core, region=-outer_surface & -active_core_maxz & +active_core_minz) active_core_universe = openmc.Universe(cells=[active_core_cell]) - if params['plotting'] == "Y": - create_universe_plot(materials_database, active_core_universe, - plot_width = 2.2 * params['Core Radius'], - num_pixels = 500, - font_size = 32, - title = "Core", - fig_size = 8, - output_file_name = "Core.png") if params['plotting'] == "Y": - create_universe_plot(materials_database, active_core_universe, - plot_width = 0.5 * params['Assembly FTF'] * params['Core Rings'] , - num_pixels = 500, - font_size = 32, - title = "Core", - fig_size = 8, - output_file_name = "Core (zoomed in).png") + core_state = ( + 'shutdown_ARI' + if params['Shutdown Margin Calc'] + else 'operation_ARO' + ) + + print( + "Saving GCMR plot:", + f"Core_{core_state}.png" + ) + + create_universe_plot( + materials_database, + active_core_universe, + plot_width=2.2 * params['Core Radius'], + num_pixels=2000, + font_size=32, + title=f"Core - {core_state}", + fig_size=8, + output_file_name=f"Core_{core_state}.png" + ) + + create_universe_plot( + materials_database, + active_core_universe, + plot_width=( + 0.5 + * params['Assembly FTF'] + * params['Core Rings'] + ), + num_pixels=2000, + font_size=32, + title=f"Core - {core_state} - Zoomed", + fig_size=8, + output_file_name=( + f"Core_{core_state}_zoomed.png" + ) + ) + # if params['plotting'] == "Y": + # create_universe_plot(materials_database, active_core_universe, + # plot_width = 2.2 * params['Core Radius'], + # num_pixels = 500, + # font_size = 32, + # title = "Core", + # fig_size = 8, + # output_file_name = "Core.png") + + # if params['plotting'] == "Y": + # create_universe_plot(materials_database, active_core_universe, + # plot_width = 0.5 * params['Assembly FTF'] * params['Core Rings'] , + # num_pixels = 500, + # font_size = 32, + # title = "Core", + # fig_size = 8, + # output_file_name = "Core (zoomed in).png") # ************************************************************************************************************************** # Sec. 6 : VOLUME INFO for Depletion @@ -412,8 +707,48 @@ def create_drums_universe_CGMR(params, absorber_thickness, drum_radius, # The compact fuel volume defined earlier has a height of 4 params['Lattice Compact Volume'] = cylinder_volume(params['Compact Fuel Radius'], 4) - outer_fuel_ring_count = 6 * (params['Core Rings'] - 1) # corner + edge assemblies not counted in assembly_number - core_fuel_cells = (assembly_number + outer_fuel_ring_count) * assembly_fuel_cells + # outer_fuel_ring_count = 6 * (params['Core Rings'] - 1) # corner + edge assemblies not counted in assembly_number + # core_fuel_cells = (assembly_number + outer_fuel_ring_count) * assembly_fuel_cells + + outer_fuel_ring_count = 6 * ( + params['Core Rings'] - 1 + ) + + central_shutdown_assembly_count = 1 + surrounding_shutdown_assembly_count = 6 + + # assembly_number includes: + # 1 center + 6 first-ring assemblies + all other inner assemblies + normal_inner_assembly_count = ( + assembly_number + - central_shutdown_assembly_count + - surrounding_shutdown_assembly_count + ) + + core_fuel_cells = ( + normal_inner_assembly_count + * assembly_fuel_cells + + + central_shutdown_assembly_count + * central_shutdown_fuel_cells + + + surrounding_shutdown_assembly_count + * surrounding_shutdown_fuel_cells + + + outer_fuel_ring_count + * assembly_fuel_cells + ) + + print("Normal inner assemblies:", normal_inner_assembly_count) + print( + "Central shutdown assemblies:", + central_shutdown_assembly_count + ) + print( + "Surrounding shutdown assemblies:", + surrounding_shutdown_assembly_count + ) + print("Total core fuel positions:", core_fuel_cells) core_compact_volume = cylinder_volume(params['Compact Fuel Radius'], params['Active Height']) * core_fuel_cells core_triso_number = core_compact_volume / params['Lattice Compact Volume'] * compact_triso_particles_number kernel_volume = sphere_volume(params['Fuel Pin Radii'][0]) diff --git a/core_design/openmc_template_LTMR.py b/core_design/openmc_template_LTMR.py index 50ec373..8d3b2ad 100644 --- a/core_design/openmc_template_LTMR.py +++ b/core_design/openmc_template_LTMR.py @@ -68,6 +68,96 @@ def create_pin_regions(params, pin_type): return regions +def create_shutdown_rod_universe(params, materials_database): + """ + Simplified shutdown channel: + + ARO: + Absorber withdrawn; channel contains NaK. + + ARI: + Enriched B4C absorber inserted. + """ + + absorber_radius = params['Shutdown Rod Absorber Radius'] + clad_radius = params['Shutdown Rod Clad Radius'] + + if absorber_radius <= 0.0: + raise ValueError( + "Shutdown Rod Absorber Radius must be greater than zero." + ) + + if absorber_radius >= clad_radius: + raise ValueError( + "Shutdown Rod Absorber Radius must be smaller than " + "Shutdown Rod Clad Radius." + ) + + if clad_radius >= params['Fuel Pin Radii'][-1]: + raise ValueError( + "Shutdown Rod Clad Radius must be smaller than " + "the existing pin outer radius." + ) + + absorber_surface = openmc.ZCylinder( + r=absorber_radius, + name='shutdown_absorber_surface' + ) + + clad_surface = openmc.ZCylinder( + r=clad_radius, + name='shutdown_clad_surface' + ) + + absorber_material = materials_database[ + params['Shutdown Rod Absorber'] + ] + + clad_material = materials_database[ + params['Shutdown Rod Cladding'] + ] + + coolant_material = materials_database[ + params['Coolant'] + ] + + if params['Shutdown Margin Calc']: + print(">>> SHUTDOWN RODS: INSERTED (B4C)") + inner_material = absorber_material + inner_name = 'shutdown_absorber_inserted' + else: + print(">>> SHUTDOWN RODS: WITHDRAWN (NaK)") + inner_material = coolant_material + inner_name = 'shutdown_channel_empty' + + inner_cell = openmc.Cell( + name=inner_name, + fill=inner_material, + region=-absorber_surface + ) + + clad_cell = openmc.Cell( + name='shutdown_rod_cladding', + fill=clad_material, + region=+absorber_surface & -clad_surface + ) + + outside_coolant_cell = openmc.Cell( + name='shutdown_rod_outer_coolant', + fill=coolant_material, + region=+clad_surface + ) + + shutdown_rod_universe = openmc.Universe( + name='shutdown_rod_universe', + cells=[ + inner_cell, + clad_cell, + outside_coolant_cell, + ] + ) + + return shutdown_rod_universe def _get_valid_drum_counts(): return [6, 12, 18, 24, 30, 36] @@ -231,20 +321,55 @@ def create_drums_universe(params, control_drum_absorber_material, control_drum_r return drums -def create_assembly_universe(params, fuel_pin_universe, moderator_pin_universe, pin_pitch, reflector_material, outer_coolant_universe): +# def create_assembly_universe(params, fuel_pin_universe, moderator_pin_universe, pin_pitch, reflector_material, outer_coolant_universe): +# """ +# Creating the universe of the fuel assembly +# @ In, params, dict, The parameters that are used to "fill in" input files with placeholders. +# @ In, fuel_pin_universe, openmc.universe.Universe +# @ In, moderator_pin_universe, openmc.universe.Universe +# @ In, pin_pitch, float, the center-to-center distance between adjacent fuel/moderator pins +# @ In, reflector_material, openmc.material.Material, the material of the outer radial reflector +# @ In, outer_coolant_universe, openmc.universe.Universe, the OpenMC universe of the coolant in the assembly +# @ out, assembly_universe, openmc.universe.Universe, the fuel assembly universe +# """ + +# assembly = openmc.HexLattice() +# assembly.center = (0., 0.) +# assembly.pitch = (pin_pitch,) +# assembly.outer = outer_coolant_universe + +# rings = copy.deepcopy(params['Pins Arrangement']) +# rings = rings[-params['Number of Rings per Assembly']:] + +# for i in range(len(rings)): +# for j in range(len(rings[i])): +# if rings[i][j] == 'FUEL': +# rings[i][j] = fuel_pin_universe +# elif rings[i][j] == 'MODERATOR': +# rings[i][j] = moderator_pin_universe + +# assembly.universes = rings + +# hex_edge_length = calculate_hex_edge_length(params) +# assembly_boundary = openmc.model.hexagonal_prism( +# edge_length=hex_edge_length, +# corner_radius=params['Fuel Pin Radii'][-1] + params["Pin Gap Distance"] +# ) + +# fuel_assembly_cell = openmc.Cell(fill=assembly, region=assembly_boundary) +# reflector_cell = openmc.Cell(fill=reflector_material, region=~assembly_boundary) + +# assembly_universe = openmc.Universe(cells=[fuel_assembly_cell, reflector_cell]) + +# return assembly_universe + +def create_assembly_universe(params, fuel_pin_universe, moderator_pin_universe, shutdown_rod_universe, pin_pitch, reflector_material, outer_coolant_universe): """ - Creating the universe of the fuel assembly - @ In, params, dict, The parameters that are used to "fill in" input files with placeholders. - @ In, fuel_pin_universe, openmc.universe.Universe - @ In, moderator_pin_universe, openmc.universe.Universe - @ In, pin_pitch, float, the center-to-center distance between adjacent fuel/moderator pins - @ In, reflector_material, openmc.material.Material, the material of the outer radial reflector - @ In, outer_coolant_universe, openmc.universe.Universe, the OpenMC universe of the coolant in the assembly - @ out, assembly_universe, openmc.universe.Universe, the fuel assembly universe + Creating the universe of the fuel assembly. """ assembly = openmc.HexLattice() - assembly.center = (0., 0.) + assembly.center = (0.0, 0.0) assembly.pitch = (pin_pitch,) assembly.outer = outer_coolant_universe @@ -255,21 +380,47 @@ def create_assembly_universe(params, fuel_pin_universe, moderator_pin_universe, for j in range(len(rings[i])): if rings[i][j] == 'FUEL': rings[i][j] = fuel_pin_universe + elif rings[i][j] == 'MODERATOR': rings[i][j] = moderator_pin_universe + elif rings[i][j] == 'SHUTDOWN': + rings[i][j] = shutdown_rod_universe + + else: + raise ValueError( + f"Unknown pin label at ring {i}, position {j}: " + f"{rings[i][j]!r}" + ) + assembly.universes = rings hex_edge_length = calculate_hex_edge_length(params) + assembly_boundary = openmc.model.hexagonal_prism( edge_length=hex_edge_length, - corner_radius=params['Fuel Pin Radii'][-1] + params["Pin Gap Distance"] + corner_radius=( + params['Fuel Pin Radii'][-1] + + params['Pin Gap Distance'] + ) + ) + + fuel_assembly_cell = openmc.Cell( + fill=assembly, + region=assembly_boundary ) - fuel_assembly_cell = openmc.Cell(fill=assembly, region=assembly_boundary) - reflector_cell = openmc.Cell(fill=reflector_material, region=~assembly_boundary) + reflector_cell = openmc.Cell( + fill=reflector_material, + region=~assembly_boundary + ) - assembly_universe = openmc.Universe(cells=[fuel_assembly_cell, reflector_cell]) + assembly_universe = openmc.Universe( + cells=[ + fuel_assembly_cell, + reflector_cell + ] + ) return assembly_universe @@ -457,6 +608,8 @@ def build_openmc_model_LTMR(params): reflector = materials_database[params['Radial Reflector']] control_drum_absorber = materials_database[params['Control Drum Absorber']] control_drum_reflector = materials_database[params['Control Drum Reflector']] + shutdown_rod_absorber = materials_database[params['Shutdown Rod Absorber']] + shutdown_rod_cladding = materials_database[params['Shutdown Rod Cladding']] # ************************************************************************************************************************** # Sec. 1.2 : Pin Cell Universes and Coolant @@ -537,6 +690,21 @@ def build_openmc_model_LTMR(params): output_file_name="moderator_pin_universe.png" ) + # Create shutdown-rod universe + shutdown_rod_universe = create_shutdown_rod_universe(params,materials_database) + + if params['plotting'] == "Y": + create_universe_plot( + materials_database, + shutdown_rod_universe, + plot_width=2.2 * params['Fuel Pin Radii'][-1], + num_pixels=500, + font_size=32, + title="Shutdown Rod Universe", + fig_size=8, + output_file_name="shutdown_rod_universe.png" + ) + # Coolant universe coolant_cell = openmc.Cell(fill=coolant) coolant_universe = openmc.Universe(cells=(coolant_cell,)) @@ -552,6 +720,7 @@ def build_openmc_model_LTMR(params): params, fuel_pin_universe, moderator_pin_universe, + shutdown_rod_universe pin_pitch, reflector, coolant_universe @@ -572,7 +741,7 @@ def build_openmc_model_LTMR(params): all_materials = ( fuel_materials + moderator_materials - + [coolant, reflector, control_drum_absorber, control_drum_reflector] + + [coolant, reflector, control_drum_absorber, control_drum_reflector, shutdown_rod_absorber, shutdown_rod_cladding] ) # Remove None materials while preserving their deterministic order @@ -603,15 +772,44 @@ def build_openmc_model_LTMR(params): core_geometry.export_to_xml() + # if params['plotting'] == "Y": + # drum_state_label = "shutdown" if params['Shutdown Margin Calc'] else "operation" + # create_universe_plot( + # materials_database, + # core_geometry, + # plot_width=2.01 * params['Core Radius'], + # num_pixels=2000, + # font_size=32, + # title="Reactor Core", + # fig_size=8, + # output_file_name=f"core_{drum_state_label}.png" + # ) + if params['plotting'] == "Y": - drum_state_label = "shutdown" if params['Shutdown Margin Calc'] else "operation" + + print( + "CORE PLOT: Shutdown Margin Calc =", + params['Shutdown Margin Calc'] + ) + + drum_state_label = ( + "shutdown" + if params['Shutdown Margin Calc'] + else "operation" + ) + + print( + "CORE PLOT FILE =", + f"core_{drum_state_label}.png" + ) + create_universe_plot( materials_database, core_geometry, plot_width=2.01 * params['Core Radius'], num_pixels=2000, font_size=32, - title="Reactor Core", + title=f"Reactor Core - {drum_state_label}", fig_size=8, output_file_name=f"core_{drum_state_label}.png" ) @@ -657,7 +855,7 @@ def build_openmc_model_LTMR(params): if 'Particles' in params.keys(): settings.particles = int(params['Particles']) else: - settings.particles = 1000 + settings.particles = 10000 if params['Isothermal Temperature Coefficients']: settings.temperature = { diff --git a/core_design/pins_arrangement.py b/core_design/pins_arrangement.py index c0456d1..12cfbac 100644 --- a/core_design/pins_arrangement.py +++ b/core_design/pins_arrangement.py @@ -2,6 +2,7 @@ # Define placeholders for the variables fuel_pin = 'FUEL' moderator_pin = 'MODERATOR' +shutdown_pin = 'SHUTDOWN' # Your list structure with placeholders LTMR_pins_arrangement = [ @@ -48,7 +49,7 @@ [moderator_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin] * 6, [fuel_pin, fuel_pin, moderator_pin, fuel_pin, fuel_pin, fuel_pin, moderator_pin, fuel_pin] * 6, -[moderator_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin] * 6, +[shutdown_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin] * 6, [fuel_pin, fuel_pin, moderator_pin, fuel_pin, moderator_pin, fuel_pin] * 6, [moderator_pin, fuel_pin, fuel_pin, fuel_pin, fuel_pin] * 6, diff --git a/core_design/utils.py b/core_design/utils.py index 9b63794..8d135ad 100644 --- a/core_design/utils.py +++ b/core_design/utils.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt import matplotlib.patches as mpatches from core_design.correction_factor import corrected_keff_2d +#from core_design.correction_factor import corrected_keff_steady_state #Use this line and comment the previuos line if you run steady state from core_design.peaking_factor import compute_pin_peaking_factors import pandas @@ -364,6 +365,32 @@ def run_depletion_analysis(params): params['Mass U238'] = mass_U238 params['Uranium Mass'] = (mass_U235 + mass_U238) / 1000 +# Use this function and comment the previous one if you run steady state +# def run_steady_state_analysis(params): +# import glob + +# openmc.run() + +# statepoint_file = sorted(glob.glob("statepoint.*.h5"))[-1] + +# keff_2d, keff_3d_corrected, p_nl_axial = corrected_keff_steady_state( +# statepoint_file, +# params['Active Height'] + 2 * params['Axial Reflector Thickness'], +# core_radius=params.get('Core Radius', np.nan) +# ) + +# params['keff 2D'] = [float(keff_2d)] +# params['keff 3D (2D corrected)'] = [float(keff_3d_corrected)] +# params['Depletion Time Steps'] = [0.0] + +# params['BOL Axial Non-Leakage Probability'] = p_nl_axial +# params['Estimated Axial Leakage (%)'] = (1.0 - p_nl_axial) * 100.0 + +# params['Fuel Lifetime'] = np.nan +# params['Mass U235'] = np.nan +# params['Mass U238'] = np.nan +# params['Uranium Mass'] = np.nan + def _sum_nuclide_mass(materials, nuclide): total_mass = 0.0 @@ -484,6 +511,7 @@ def _run_isothermal_temperature_coefficients(build_openmc_model, params): openmc_plugin = watts.PluginOpenMC(build_openmc_model, show_stderr=True) openmc_plugin(params, function=lambda: run_depletion_analysis(params)) + #openmc_plugin(params, function=lambda: run_steady_state_analysis(params)) #Use this line and comment the previuos line if you run steady state params['keff 2D high temp'] = params['keff 2D'] params['keff 3D (2D corrected) high temp'] = params['keff 3D (2D corrected)'] @@ -491,6 +519,7 @@ def _run_isothermal_temperature_coefficients(build_openmc_model, params): openmc_plugin = watts.PluginOpenMC(build_openmc_model, show_stderr=True) openmc_plugin(params, function=lambda: run_depletion_analysis(params)) + #openmc_plugin(params, function=lambda: run_steady_state_analysis(params)) #Use this line and comment the previuos line if you run steady state params['keff 2D ARO'] = params['keff 2D'] params['keff 3D (2D corrected) ARO'] = params['keff 3D (2D corrected)'] @@ -544,6 +573,7 @@ def run_openmc(build_openmc_model, heat_flux_monitor, params): params['Common Temperature'] = params['Cold Shutdown Temperature'] openmc_plugin = watts.PluginOpenMC(build_openmc_model, show_stderr=True) openmc_plugin(params, function=lambda: run_depletion_analysis(params)) + #openmc_plugin(params, function=lambda: run_steady_state_analysis(params)) #Use this line and comment the previuos line if you run steady state params['keff 2D ARI'] = params['keff 2D'] params['keff 3D (2D corrected) ARI'] = params['keff 3D (2D corrected)'] @@ -551,6 +581,7 @@ def run_openmc(build_openmc_model, heat_flux_monitor, params): params['Common Temperature'] = original_common_temperature openmc_plugin = watts.PluginOpenMC(build_openmc_model, show_stderr=True) openmc_plugin(params, function=lambda: run_depletion_analysis(params)) + #openmc_plugin(params, function=lambda: run_steady_state_analysis(params)) #Use this line and comment the previuos line if you run steady state params['keff 2D ARO'] = params['keff 2D'] params['keff 3D (2D corrected) ARO'] = params['keff 3D (2D corrected)'] @@ -583,6 +614,7 @@ def run_openmc(build_openmc_model, heat_flux_monitor, params): openmc_plugin = watts.PluginOpenMC(build_openmc_model, show_stderr=True) openmc_plugin(params, function=lambda: run_depletion_analysis(params)) + #openmc_plugin(params, function=lambda: run_steady_state_analysis(params)) #Use this line and comment the previuos line if you run steady state params['keff 2D ARO'] = params['keff 2D'] params['keff 3D (2D corrected) ARO'] = params['keff 3D (2D corrected)'] diff --git a/examples/watts_exec_GCMR_Design_A.py b/examples/watts_exec_GCMR_Design_A.py index 16f5a8d..51bbc80 100644 --- a/examples/watts_exec_GCMR_Design_A.py +++ b/examples/watts_exec_GCMR_Design_A.py @@ -44,7 +44,7 @@ def update_params(updates): 'reactor type': "GCMR", # LTMR or GCMR 'TRISO Fueled': "Yes", 'Fuel': 'UCO', - 'Enrichment': 0.1975, # The enrichment is a fraction. It has to be between 0 and 1 + 'Enrichment': 0.12, # The enrichment is a fraction. It has to be between 0 and 1 'UO2 atom fraction': 0.7, # Mixing UO2 and UC by atom fraction 'Radial Reflector': 'Graphite', 'Axial Reflector': 'Graphite', @@ -67,14 +67,24 @@ def update_params(updates): 'Fuel Pin Materials': ['UCO', 'buffer_graphite', 'PyC', 'SiC', 'PyC'], 'Fuel Pin Radii': [0.0250, 0.0350, 0.0390, 0.0425, 0.0465], # cm # https://art.inl.gov/NRC%20Training%202019/04_TRISO_Fuel.pdf 'Compact Fuel Radius': 0.6225, # cm # The radius of the area that is occupied by the TRISO particles (fuel compact/ fuel element) - 'Packing Fraction': 0.3, + 'Packing Fraction': 0.4, # Coolant channel and booster dimensions 'Coolant Channel Radius': 0.35, # cm - 'Moderator Booster Radii': [0.55], # cm + 'Moderator Booster Radii': [0.5], # cm 'Lattice Pitch': 2.25, 'Assembly Rings': 6, 'Core Rings': 5, + + # Central assembly + 'Central Shutdown Rod Radius': 0.85, # cm + 'Central Shutdown Rod Ring': 2, + 'Central Shutdown Rod Count': 12, + + # Six assemblies surrounding the center + 'Surrounding Shutdown Rod Radius': 0.45, # cm + 'Surrounding Shutdown Rod Ring': 2, + 'Surrounding Shutdown Rod Count': 2, }) params['Assembly FTF'] = params['Lattice Pitch']*(params['Assembly Rings']-1)*np.sqrt(3) # if unspecified, radial reflector thickness defaults to just cover the drums, @@ -82,7 +92,7 @@ def update_params(updates): # params['Radial Reflector Thickness'] = 27.393 # cm # radial reflector # params['Axial Reflector Thickness'] = params['Radial Reflector Thickness'] # cm # params['Core Radius'] = params['Assembly FTF']*params['Core Rings'] + params['Radial Reflector Thickness'] -params['Active Height'] = 250 +params['Active Height'] = 200 # ************************************************************************************************************************** # Sec. 3: Control Drums diff --git a/examples/watts_exec_LTMR.py b/examples/watts_exec_LTMR.py index 50cbaa9..75c8854 100644 --- a/examples/watts_exec_LTMR.py +++ b/examples/watts_exec_LTMR.py @@ -45,7 +45,7 @@ def update_params(updates): 'reactor type': "LTMR", # LTMR or GCMR 'TRISO Fueled': "No", 'Fuel': 'UZrH_alloy', - 'Enrichment': 0.1975, # Fraction between 0 and 1 + 'Enrichment': 0.05, # Fraction between 0 and 1 "H_Zr_ratio": 1.6, # Proportion of hydrogen to zirconium atoms 'U_met_wo': 0.3, # Weight ratio of Uranium to total fuel weight (less than 1) 'er_wo': 0, # Erbium (burnable poison) @@ -70,14 +70,14 @@ def update_params(updates): 'Moderator Pin Radii': [1.5367, 1.5875], # [params['Moderator Pin Inner Radius'], params['Fuel Pin Radii'][-1]] "Pin Gap Distance": 0.1, # cm 'Pins Arrangement': LTMR_pins_arrangement, - 'Number of Rings per Assembly': 12, # the number of rings can be 12 or lower as long as the heat flux criteria is not violated + 'Number of Rings per Assembly': 14, # the number of rings can be 12 or lower as long as the heat flux criteria is not violated 'Radial Reflector Thickness': 14, # cm }) params['Lattice Apothem'] = calculate_hex_apothem(params) params['Lattice Radius'] = params['Lattice Apothem'] params['Assembly FTF'] = 2 * params['Lattice Apothem'] -params['Active Height'] = 78.4 +params['Active Height'] = 120 params['Axial Reflector Thickness'] = params['Radial Reflector Thickness'] # cm params['Fuel Pin Count'] = calculate_pins_in_assembly(params, "FUEL") params['Moderator Pin Count'] = calculate_pins_in_assembly(params, "MODERATOR") @@ -89,11 +89,21 @@ def update_params(updates): # ************************************************************************************************************************** update_params({ - 'Number of Drums': 12, + 'Number of Drums': 6, # When the user does not specify the drum radius, the code automatically sets it to the largest allowable value that avoids drum overlap #'Drum Radius': 9.016, #, # cm 'Drum Absorber Thickness': 1, # cm 'Drum Absorber Arc Degrees': 120, + 'Drum Height': params['Active Height'] + 2*params['Axial Reflector Thickness'], + 'Shutdown Rod Absorber': 'B4C_enriched', + 'Shutdown Rod Cladding': 'SS304', + + # Must fit inside the existing pin envelope + 'Shutdown Rod Absorber Radius': 1.30, # cm + 'Shutdown Rod Clad Radius': 1.50, # cm + + # Number of moderator positions changed to shutdown channels + 'Number of Shutdown Rods': 6, }) update_ltmr_reflector_geometry_from_drums(params)