diff options
Diffstat (limited to 'hpmr.py')
| -rwxr-xr-x | hpmr.py | 116 |
1 files changed, 64 insertions, 52 deletions
@@ -9,9 +9,6 @@ import time # HP-MR Model in OpenMC # C.Moore, <connor.moore@psi.ch> -# For LCLRS -#openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml" - ############ COMMAND LINE ARGUMENTS ############ cli_parser = argparse.ArgumentParser( prog="HP-RM TRISO Homogenization Study", @@ -28,9 +25,7 @@ if cli_args.lclrs: openmc.config["cross_sections"] = "/libs/endfb81_official/cross_sections.xml" ############ INTPUT-OUTPUT MAPPING ############ -path = sys.argv[0] - -with open(path, "rb") as file: +with open(sys.argv[0], "rb") as file: contents = file.read() filehash = hashlib.md5(contents).hexdigest() @@ -430,13 +425,13 @@ mat_core_b4c_central.set_density("g/cc",1.25) ############ GEOMETRY ############ # All measurements are in cm unless specified -geo_constant_pin_pitch = 2.3 -geo_constant_lattice_pitch = 26.752/3**0.5 -geo_constant_triso_radius = 1.0 -geo_constant_triso_packing = 0.40 -geo_constant_height = 1.0 -geo_constant_z_max = openmc.ZPlane(z0=+geo_constant_height,boundary_type="reflective") -geo_constant_z_min = openmc.ZPlane(z0=-geo_constant_height,boundary_type="reflective") +geo_const_pin_pitch = 2.3 +geo_const_lattice_pitch = 26.752/3**0.5 +geo_const_triso_radius = 1.0 +geo_const_triso_packing = 0.40 +geo_const_height = 1.0 +geo_const_z_max = openmc.ZPlane(z0=+geo_const_height,boundary_type="reflective") +geo_const_z_min = openmc.ZPlane(z0=-geo_const_height,boundary_type="reflective") geo_none_universe = openmc.Universe(name="Empty Universe",cells=[openmc.Cell(name="Empty Cell",fill=None)]) @@ -452,14 +447,14 @@ geo_triso_cells = [openmc.Cell(name="TRISO Kernel UCO",fill=mat_triso_uco,region geo_triso_universe = openmc.Universe(name="TRISO Kernel Universe",cells=geo_triso_cells) ## TRISO compact ## -geo_triso_compact_region = -openmc.ZCylinder(r=geo_constant_triso_radius,x0=0.0,y0=0.0)&+geo_constant_z_min&-geo_constant_z_max +geo_triso_compact_region = -openmc.ZCylinder(r=geo_const_triso_radius,x0=0.0,y0=0.0)&+geo_const_z_min&-geo_const_z_max geo_triso_compact_domain = (10,10,10) # Generate packing geo_triso_compact_packing = openmc.model.pack_spheres( radius=geo_triso_spheres[4].r, region=geo_triso_compact_region, - pf=geo_constant_triso_packing + pf=geo_const_triso_packing ) # Create the TRISO particles @@ -468,8 +463,8 @@ geo_triso_compact_trisos = [openmc.model.TRISO(geo_triso_spheres[4].r, geo_triso # Support with a lattice for domain decomposition geo_triso_compact_lattice = openmc.model.create_triso_lattice( trisos=geo_triso_compact_trisos, - lower_left=(-geo_constant_triso_radius,-geo_constant_triso_radius,-geo_constant_height), - pitch=(geo_constant_triso_radius*2,geo_constant_triso_radius*2,geo_constant_height*2), + lower_left=(-geo_const_triso_radius,-geo_const_triso_radius,-geo_const_height), + pitch=(geo_const_triso_radius*2,geo_const_triso_radius*2,geo_const_height*2), shape=geo_triso_compact_domain, background=mat_graphite_matrix ) @@ -483,19 +478,19 @@ geo_triso_compact_universe = openmc.Universe(name="TRISO Compact Universe",cells ## Homogenized TRISO Compact ## geo_vwh_fracs = [] geo_rpt_fracs = [] -geo_triso_rpt_region = -openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0)&+geo_constant_z_min&-geo_constant_z_max +geo_triso_rpt_region = -openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0)&+geo_const_z_min&-geo_const_z_max geo_n_triso = len(geo_triso_compact_packing) print("\n TRISO Stacking Information from Packing Calculation") print(f" -> Total TRISO pebbles packed: {geo_n_triso}") print(f" -> Total TRISO volume in compact: {geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3:.4f}") -print(f" -> Minimum RPT cylinder radius: {(geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_constant_height))**0.5:.4f}") -print(f" -> Maximum RRPT ring radius: {(1.15**2 - geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_constant_height))**0.5:.4f}") +print(f" -> Minimum RPT cylinder radius: {(geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_const_height))**0.5:.4f}") +print(f" -> Maximum RRPT ring radius: {(1.15**2 - geo_n_triso*4/3*geo_triso_spheres[4].r**3/(2*geo_const_height))**0.5:.4f}") print(" TRISO Stacking calculation complete!\n") # Get volume fractions for mixing -geo_triso_compact_vwh_v_total = math.pi*geo_constant_triso_radius**2*2*geo_constant_height +geo_triso_compact_vwh_v_total = math.pi*geo_const_triso_radius**2*2*geo_const_height geo_triso_compact_vwh_v_matrix = geo_triso_compact_vwh_v_total - geo_n_triso * 4/3*math.pi*geo_triso_spheres[4].r**3 -geo_triso_compact_rpt_v_total = math.pi*cli_args.radius**2*2*geo_constant_height +geo_triso_compact_rpt_v_total = math.pi*cli_args.radius**2*2*geo_const_height geo_triso_compact_rpt_v_matrix = geo_triso_compact_rpt_v_total - geo_n_triso * 4/3*math.pi*geo_triso_spheres[4].r**3 geo_v_layers = [] @@ -558,9 +553,9 @@ for mat, tsl_list in mat_tsl_saved.items(): mat._sab = tsl_list # RRPT transformation -geo_rrpt_r_out = math.sqrt(4*geo_n_triso*geo_triso_spheres[4].r**3/(6*geo_constant_height) + cli_args.radius**2) +geo_rrpt_r_out = math.sqrt(4*geo_n_triso*geo_triso_spheres[4].r**3/(6*geo_const_height) + cli_args.radius**2) geo_rrpt_region = +openmc.ZCylinder(r=cli_args.radius,x0=0.0,y0=0.0) & -openmc.ZCylinder(r=geo_rrpt_r_out,x0=0.0,y0=0.0) \ - & +geo_constant_z_min & -geo_constant_z_max + & +geo_const_z_min & -geo_const_z_max # Cast into cells and a universe geo_triso_vwh_fuel_cell = openmc.Cell(name="TRISO VWH Fuel Cell",fill=mat_triso_vwh,region=geo_triso_compact_region) @@ -609,7 +604,7 @@ geo_mod_universe = openmc.Universe(name="Moderator Pin Universe",cells=geo_mod_c geo_assembly_lattice = openmc.HexLattice() geo_assembly_lattice.center=(0,0) -geo_assembly_lattice.pitch=(geo_constant_pin_pitch,) +geo_assembly_lattice.pitch=(geo_const_pin_pitch,) geo_assembly_lattice.outer=openmc.Universe(cells=[openmc.Cell(fill=mat_graphite_matrix)]) # Create the lattice structure @@ -629,17 +624,23 @@ geo_assembly_lattice.universes=geo_assembly_lattice_map # Bound it to a unit assembly geo_assembly_region = openmc.model.HexagonalPrism( - edge_length=geo_constant_lattice_pitch, + edge_length=geo_const_lattice_pitch, orientation="y", origin=(0.0,0.0), boundary_type="transmission" if cli_args.core_level else "reflective" ) # Clip the lattice to the space and make a universe -geo_assembly_cell = openmc.Cell(name="Unit Assembly Cell",fill=geo_assembly_lattice,region=-geo_assembly_region&+geo_constant_z_min&-geo_constant_z_max) +geo_assembly_cell = openmc.Cell(name="Unit Assembly Cell",fill=geo_assembly_lattice,region=-geo_assembly_region&+geo_const_z_min&-geo_const_z_max) geo_assembly_universe = openmc.Universe(cells=[geo_assembly_cell]) ## Core-level Geometry ## +geo_const_core_z_max = openmc.ZPlane(z0=80) +geo_const_core_z_min = openmc.ZPlane(z0=-80) +geo_const_refl_z_max = openmc.ZPlane(z0=100, boundary_type="vacuum") +geo_const_refl_z_min = openmc.ZPlane(z0=-100, boundary_type="vacuum") +geo_const_core_cyl = openmc.ZCylinder(r=115, x0=0.0, y0=0.0, boundary_type="vacuum") + # Reflector region cell and universe geo_refl_cell = openmc.Cell(name="Core Reflector Cell",fill=mat_core_be_reflector) geo_refl_universe = openmc.Universe(name="Core Reflector Universe",cells=[geo_refl_cell]) @@ -652,41 +653,54 @@ geo_cd_universe = openmc.Universe(name="Core Control Drum Universe",cells=[geo_c geo_core_rad_lattice=openmc.HexLattice() geo_core_rad_lattice.orientation="x" geo_core_rad_lattice.center=(0,0) -geo_core_rad_lattice.pitch=(geo_constant_lattice_pitch*3**0.5,) +geo_core_rad_lattice.pitch=(geo_const_lattice_pitch*3**0.5,) geo_core_rad_lattice.outer=geo_none_universe # Assemble lattice map geo_core_rad_lattice_map=[] -geo_core_rad_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i,_ in enumerate(range(24))]) -#geo_core_rad_lattice_map.append([geo_assembly_universe]*18) -geo_core_rad_lattice_map.append([geo_cd_universe if i%3==0 else geo_assembly_universe for i,_ in enumerate(range(18))]) -#print(f"{type(geo_core_rad_lattice_map[1]) = }") +geo_core_rad_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)]) +geo_core_rad_lattice_map.append([geo_cd_universe if i%3==0 else geo_assembly_universe for i in range(18)]) geo_core_rad_lattice_map.append([geo_assembly_universe]*12) geo_core_rad_lattice_map.append([geo_assembly_universe]*6) geo_core_rad_lattice_map.append([geo_refl_universe]) geo_core_rad_lattice.universes=geo_core_rad_lattice_map +# Cell and universe +geo_core_rad_cell=openmc.Cell(name="Core Radial Slice Cell",fill=geo_core_rad_lattice) +geo_core_rad_universe = openmc.Universe(cells=[geo_core_rad_cell]) -# Reflector lattice map + +# Axial Reflector lattice geo_core_rad_refl_lattice=openmc.HexLattice() geo_core_rad_refl_lattice.center=(0,0) -geo_core_rad_refl_lattice.pitch=(geo_constant_lattice_pitch,) +geo_core_rad_refl_lattice.orientation="x" +geo_core_rad_refl_lattice.pitch=(geo_const_lattice_pitch*3**0.5,) geo_core_rad_refl_lattice.outer=geo_none_universe -geo_core_rad_refl_lattice.universes=[[geo_refl_universe]*n for n in [24,18,12,6,1]] -geo_core_rad_cell=openmc.Cell(name="Core Radial Slice Cell",fill=geo_core_rad_lattice) -geo_core_rad_universe = openmc.Universe(cells=[geo_core_rad_cell]) +# Axial Reflector lattice map +geo_core_rad_refl_lattice_map = [] +geo_core_rad_refl_lattice_map.append([geo_cd_universe if i%4==2 else geo_none_universe if i%4==0 else geo_refl_universe for i in range(24)]) +geo_core_rad_refl_lattice_map.append([geo_cd_universe if i%3==0 else geo_refl_universe for i in range(18)]) +geo_core_rad_refl_lattice_map.append([geo_refl_universe]*12) +geo_core_rad_refl_lattice_map.append([geo_refl_universe]*6) +geo_core_rad_refl_lattice_map.append([geo_refl_universe]) +geo_core_rad_refl_lattice.universes = geo_core_rad_refl_lattice_map + +geo_core_refl_top_cell = openmc.Cell(name="Top Axial Reflector Cell",fill=geo_core_rad_refl_lattice, + region=-geo_const_core_cyl & +geo_const_core_z_max & -geo_const_refl_z_max) +geo_core_refl_bot_cell = openmc.Cell(name="Bottom Axial Reflector Cell",fill=geo_core_rad_refl_lattice, + region=-geo_const_core_cyl & -geo_const_core_z_min & +geo_const_refl_z_min) # 3D Stacking Lattice geo_core_axial_lattice = openmc.RectLattice() -geo_core_axial_lattice.lower_left = (-sys.float_info.max/2, -sys.float_info.max/2, -int(160//geo_constant_height)/2*geo_constant_height) -#geo_core_axial_lattice.center=(0.0,0.0,0.0) -geo_core_axial_lattice.pitch=(sys.float_info.max, sys.float_info.max, geo_constant_height) -geo_core_axial_lattice.universes = [[[geo_core_rad_universe]] for _ in range(int(160//geo_constant_height))] +geo_core_axial_lattice.lower_left = (-sys.float_info.max/2, -sys.float_info.max/2, -int(160//geo_const_height)/2*geo_const_height) +geo_core_axial_lattice.pitch=(sys.float_info.max, sys.float_info.max, geo_const_height) +geo_core_axial_lattice.universes = [[[geo_core_rad_universe]] for _ in range(int(160//geo_const_height))] geo_core_axial_lattice.outer=geo_none_universe -geo_core_cell = openmc.Cell(name="Core 3D Cell",fill=geo_core_axial_lattice) -geo_core_universe = openmc.Universe(cells=[geo_core_cell]) +geo_core_cell = openmc.Cell(name="Core 3D Cell",fill=geo_core_axial_lattice, region=-geo_const_core_cyl & +geo_const_core_z_min & -geo_const_core_z_max) +geo_core_universe = openmc.Universe(cells=[geo_core_cell, geo_core_refl_top_cell, geo_core_refl_bot_cell]) +#geo_core_universe = openmc.Universe(cells=[geo_core_refl_bot_cell]) ## Final Geometry Universe geometry = openmc.Geometry(geo_core_universe if cli_args.core_level else geo_assembly_universe) @@ -696,17 +710,15 @@ settings = openmc.Settings() settings.particles = 100000 settings.batches = 200 settings.inactive = 50 +#settings.verbosity = 10 -# Source sampling -#set_source_pts = [] -#for loc in geo_triso_compact_packing[::100]: -# point = openmc.stats.Point(xyz=loc) -# set_source_pts.append(openmc.IndependentSource(space=point)) +set_dyn_radial = geo_const_core_cyl.r if cli_args.core_level else geo_const_lattice_pitch +set_dyn_axial = 180 if cli_args.core_level else geo_const_height set_source_pts = openmc.IndependentSource( space=openmc.stats.Box( - lower_left = (-geo_constant_lattice_pitch, -geo_constant_lattice_pitch, -geo_constant_height), - upper_right = (geo_constant_lattice_pitch, geo_constant_lattice_pitch, geo_constant_height) + lower_left = (-set_dyn_radial, -set_dyn_radial, -set_dyn_axial), + upper_right = (set_dyn_radial, set_dyn_radial, set_dyn_axial) )) @@ -714,8 +726,8 @@ settings.source = set_source_pts # Mesh for Shannon Entropy set_entropy_mesh = openmc.RegularMesh() -set_entropy_mesh.lower_left=(-geo_constant_lattice_pitch,-geo_constant_lattice_pitch,-geo_constant_height) -set_entropy_mesh.upper_right=(+geo_constant_lattice_pitch,+geo_constant_lattice_pitch,+geo_constant_height) +set_entropy_mesh.lower_left=(-set_dyn_radial,-set_dyn_radial,-set_dyn_axial) +set_entropy_mesh.upper_right=(+set_dyn_radial,+set_dyn_radial,+set_dyn_axial) set_entropy_mesh.dimension=(30,30,30) settings.entropy_mesh=set_entropy_mesh |