#!/bin/python3 import openmc import math import argparse import hashlib import sys import time ############ COMMAND LINE ARGUMENTS ############ cli_parser = argparse.ArgumentParser( prog="HP-RM TRISO Homogenization Study", description="Program to investigate various homogenization techniques for TRISO compacts in the HP-MR benchmark." ) cli_parser.add_argument("-t","--technique",default="vwh",choices=["none","vwh","rpt","rrpt"],help="What homogenization technique to use (if any) [none/vwh/rpt/rrpt], default VWH") cli_parser.add_argument("-r","--radius",default="0.8",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.5 cm. In the case of RRPT, this is the inner radius.") cli_args = cli_parser.parse_args() ############ INTPUT-OUTPUT MAPPING ############ path = sys.argv[0] with open(path, "rb") as file: contents = file.read() filehash = hashlib.md5(contents).hexdigest() print("="*60) print(f"Run date CET: {time.ctime()}") print(f"Input file MD5: {filehash}") print(f"Input arguments: {cli_args}") print("="*60) ############ MATERIALS ############ # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB # # Materials marked PNNL-15870 are from the 2nd revision of the report "Compendium of Material # Composition Data for Radiation Transport Modeling " by US Dept. of Homeland Security and PNNL # mat_list = [] mat_temp_global = 700.0 mat_temp_fuel = 1000.0 ## 100s - Matrix Graphite ## mat_graphite_matrix = openmc.Material(100,name="Graphite Matrix (between pins)") mat_graphite_matrix.add_nuclide("C12",0.9999997) mat_graphite_matrix.add_nuclide("B10",3e-7) mat_graphite_matrix.set_density("g/cc",1.806) mat_graphite_matrix.temperature = mat_temp_global mat_list.append(mat_graphite_matrix) ## 200s - TRISO Pins ## mat_triso_uco = openmc.Material(200,name="TRISO Pin UCO Kernel") mat_triso_uco.add_nuclide("U235",0.068794) mat_triso_uco.add_nuclide("U238",0.27604) mat_triso_uco.add_nuclide("C12",0.13793) mat_triso_uco.add_nuclide("O16",0.51724) mat_triso_uco.set_density("g/cc",10.744) mat_triso_uco.temperature = mat_temp_fuel mat_triso_uco.add_s_alpha_beta("c_Graphite") mat_list.append(mat_triso_uco) mat_triso_buffer = openmc.Material(201,name="TRISO Pin Buffer") mat_triso_buffer.add_nuclide("C12",1.0) mat_triso_buffer.set_density("g/cc",1.04) mat_triso_buffer.temperature = mat_temp_fuel mat_triso_buffer.add_s_alpha_beta("c_Graphite") mat_list.append(mat_triso_buffer) mat_triso_pyc1 = openmc.Material(202,name="TRISO Pin PyC1") mat_triso_pyc1.add_nuclide("C12",1.0) mat_triso_pyc1.set_density("g/cc",1.882) mat_triso_pyc1.temperature = mat_temp_fuel mat_triso_pyc1.add_s_alpha_beta("c_Graphite") mat_list.append(mat_triso_pyc1) mat_triso_sic = openmc.Material(203,name="TRISO Pin SiC") mat_triso_sic.add_nuclide("Si28",0.4611) mat_triso_sic.add_nuclide("Si29",0.0234) mat_triso_sic.add_nuclide("Si30",0.0154) mat_triso_sic.add_nuclide("C12",0.5) mat_triso_sic.set_density("g/cc",3.171) mat_triso_sic.temperature = mat_temp_fuel mat_triso_sic.add_s_alpha_beta("c_SiC") mat_list.append(mat_triso_sic) mat_triso_pyc2 = openmc.Material(204,name="TRISO Pin PyC2") mat_triso_pyc2.add_nuclide("C12",1.0) mat_triso_pyc2.add_s_alpha_beta("c_Graphite") mat_triso_pyc2.temperature = mat_temp_fuel mat_triso_pyc2.set_density("g/cc",1.882) mat_list.append(mat_triso_pyc2) ## 300 - Moderator Pins ## mat_mod_helium = openmc.Material(300,name="Moderator Pin Helium") mat_mod_helium.add_element("He4",1.0) mat_mod_helium.set_density("g/cc",0.18e-3) mat_list.append(mat_mod_helium) # From the offical Serpent input mat_mod_ss316 = openmc.Material(301,name="Moderator Pin SS316") mat_mod_ss316.add_nuclide('C12', 0.001901) mat_mod_ss316.add_nuclide('Si28', 0.0092693) mat_mod_ss316.add_nuclide('Si29', 4.7251e-4) mat_mod_ss316.add_nuclide('Si30', 3.1166e-4) mat_mod_ss316.add_nuclide('P31', 4.1322e-4) mat_mod_ss316.add_nuclide('S32', 2.471e-4) mat_mod_ss316.add_nuclide('S33', 1.9511e-6) mat_mod_ss316.add_nuclide('S34', 1.1056e-5) mat_mod_ss316.add_nuclide('S36', 3.0016e-8) mat_mod_ss316.add_nuclide('Cr50', 7.9116e-3) mat_mod_ss316.add_nuclide('Cr52', 0.15257) mat_mod_ss316.add_nuclide('Cr53', 0.01730) mat_mod_ss316.add_nuclide('Cr54', 4.3063e-3) mat_mod_ss316.add_nuclide('Mn55', 0.01028) mat_mod_ss316.add_nuclide('Fe54', 0.039029) mat_mod_ss316.add_nuclide('Fe56', 0.61213) mat_mod_ss316.add_nuclide('Fe57', 0.014144) mat_mod_ss316.add_nuclide('Fe58', 1.3343e-3) mat_mod_ss316.add_nuclide('Ni58', 0.077516) mat_mod_ss316.add_nuclide('Ni60', 0.029859) mat_mod_ss316.add_nuclide('Ni61', 1.2981e-3) mat_mod_ss316.add_nuclide('Ni62', 4.1389e-3) mat_mod_ss316.add_nuclide('Ni64', 1.0544e-3) mat_mod_ss316.add_nuclide('Mo92', 2.1259e-3) mat_mod_ss316.add_nuclide('Mo94', 1.3299e-3) mat_mod_ss316.add_nuclide('Mo95', 2.303e-3) mat_mod_ss316.add_nuclide('Mo96', 2.4191e-3) mat_mod_ss316.add_nuclide('Mo97', 1.3903e-3) mat_mod_ss316.add_nuclide('Mo98', 3.5249e-3) mat_mod_ss316.add_nuclide('Mo100',1.4135e-3) mat_mod_ss316.set_density("g/cc",7.67) mat_list.append(mat_mod_ss316) # Official Serpent smeared composition mat_mod_ss316_he_smear = openmc.Material(302,name="Moderator Pin SS316/He Smear") mat_mod_ss316_he_smear.add_nuclide('He4',1.983e-5) mat_mod_ss316_he_smear.add_nuclide('C12',4.349e-5) mat_mod_ss316_he_smear.add_nuclide('Si28',2.121e-4) mat_mod_ss316_he_smear.add_nuclide('Si29',1.081e-5) mat_mod_ss316_he_smear.add_nuclide('Si30',7.13e-6) mat_mod_ss316_he_smear.add_nuclide('P31',9.454e-6) mat_mod_ss316_he_smear.add_nuclide('S32',5.653e-6) mat_mod_ss316_he_smear.add_nuclide('S33',4.464e-8) mat_mod_ss316_he_smear.add_nuclide('S34',2.53e-7) mat_mod_ss316_he_smear.add_nuclide('S36',6.867e-10) mat_mod_ss316_he_smear.add_nuclide('Cr50',1.81e-4) mat_mod_ss316_he_smear.add_nuclide('Cr52',3.491e-3) mat_mod_ss316_he_smear.add_nuclide('Cr53',3.958e-4) mat_mod_ss316_he_smear.add_nuclide('Cr54',9.852e-5) mat_mod_ss316_he_smear.add_nuclide('Mn55',2.352e-4) mat_mod_ss316_he_smear.add_nuclide('Fe54',8.929e-4) mat_mod_ss316_he_smear.add_nuclide('Fe56',0.014) mat_mod_ss316_he_smear.add_nuclide('Fe57',3.236e-4) mat_mod_ss316_he_smear.add_nuclide('Fe58',3.053e-5) mat_mod_ss316_he_smear.add_nuclide('Ni58',1.773e-3) mat_mod_ss316_he_smear.add_nuclide('Ni60',6.831e-4) mat_mod_ss316_he_smear.add_nuclide('Ni61',2.97e-5) mat_mod_ss316_he_smear.add_nuclide('Ni62',9.469e-5) mat_mod_ss316_he_smear.add_nuclide('Ni64',2.412e-5) mat_mod_ss316_he_smear.add_nuclide('Mo92',4.864e-5) mat_mod_ss316_he_smear.add_nuclide('Mo94',3.043e-5) mat_mod_ss316_he_smear.add_nuclide('Mo95',5.269e-5) mat_mod_ss316_he_smear.add_nuclide('Mo96',5.535e-5) mat_mod_ss316_he_smear.add_nuclide('Mo97',3.181e-5) mat_mod_ss316_he_smear.add_nuclide('Mo98',8.065e-5) mat_mod_ss316_he_smear.add_nuclide('Mo100',3.234e-5) mat_mod_ss316_he_smear.set_density("atom/b-cm",0.0229) mat_mod_ss316_he_smear.temperature = mat_temp_global mat_list.add_nuclide(mat_mod_ss316_he_smear) # YH2 from official Serpent mat_mod_yh2 = openmc.Material(303,name="Moderater Pin YH2") mat_mod_yh2.add_nuclide("Y89",0.357142857) mat_mod_yh2.add_nuclide("H1",0.642857143) mat_mod_yh2.set_density("g/cc",4.085) mat_mod_yh2.temperature = mat_temp_global mat_mod_yh2.add_s_alpha_beta("c_H_in_YH2") mat_mod_yh2.add_s_alpha_beta("c_Y_in_YH2") mat_list.append(mat_mod_yh2) ## 400s - Heat Pipes ## mat_pipe_helium = openmc.Material(name="Heat Pipe Helium") mat_pipe_helium.add_element("He4",1.0) mat_pipe_helium.set_density("g/cc",0.18e-3) mat_list.append(mat_pipe_helium) # From the offical Serpent input mat_pipe_ss316 = openmc.Material(401,name="Heat Pipe SS316") mat_pipe_ss316.add_nuclide('C12', 0.001901) mat_pipe_ss316.add_nuclide('Si28', 0.0092693) mat_pipe_ss316.add_nuclide('Si29', 4.7251e-4) mat_pipe_ss316.add_nuclide('Si30', 3.1166e-4) mat_pipe_ss316.add_nuclide('P31', 4.1322e-4) mat_pipe_ss316.add_nuclide('S32', 2.471e-4) mat_pipe_ss316.add_nuclide('S33', 1.9511e-6) mat_pipe_ss316.add_nuclide('S34', 1.1056e-5) mat_pipe_ss316.add_nuclide('S36', 3.0016e-8) mat_pipe_ss316.add_nuclide('Cr50', 7.9116e-3) mat_pipe_ss316.add_nuclide('Cr52', 0.15257) mat_pipe_ss316.add_nuclide('Cr53', 0.01730) mat_pipe_ss316.add_nuclide('Cr54', 4.3063e-3) mat_pipe_ss316.add_nuclide('Mn55', 0.01028) mat_pipe_ss316.add_nuclide('Fe54', 0.039029) mat_pipe_ss316.add_nuclide('Fe56', 0.61213) mat_pipe_ss316.add_nuclide('Fe57', 0.014144) mat_pipe_ss316.add_nuclide('Fe58', 1.3343e-3) mat_pipe_ss316.add_nuclide('Ni58', 0.077516) mat_pipe_ss316.add_nuclide('Ni60', 0.029859) mat_pipe_ss316.add_nuclide('Ni61', 1.2981e-3) mat_pipe_ss316.add_nuclide('Ni62', 4.1389e-3) mat_pipe_ss316.add_nuclide('Ni64', 1.0544e-3) mat_pipe_ss316.add_nuclide('Mo92', 2.1259e-3) mat_pipe_ss316.add_nuclide('Mo94', 1.3299e-3) mat_pipe_ss316.add_nuclide('Mo95', 2.303e-3) mat_pipe_ss316.add_nuclide('Mo96', 2.4191e-3) mat_pipe_ss316.add_nuclide('Mo97', 1.3903e-3) mat_pipe_ss316.add_nuclide('Mo98', 3.5249e-3) mat_pipe_ss316.add_nuclide('Mo100',1.4135e-3) mat_pipe_ss316.set_density("g/cc",7.67) mat_list.append(mat_pipe_ss316) mat_pipe_k_liquid = openmc.Material(name="Heat Pipe Potassium (Liquid)") mat_pipe_k_liquid.add_element("K",1.0) mat_pipe_k_liquid.set_density("g/cc",0.705) mat_list.append(mat_pipe_k_liquid) mat_pipe_wick = openmc.Material.mix_materials( name = "Heat Pipe Wick", materials = [mat_pipe_k_liquid, mat_pipe_ss316], fracs = [0.7, 0.3], percent_type = 'vo' ) mat_list.append(mat_pipe_wick) mat_pipe_k_gas = openmc.Material(name="Heat Pipe Potassium (Gaseous)") mat_pipe_k_gas.add_element("K",1.0) mat_pipe_k_gas.set_density("g/cc",1.11e-4) mat_list.append(mat_pipe_k_gas) ############ 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") ### Inidividual Pin Universes ### ## TRISO fuel kernel ## geo_triso_spheres = [openmc.Sphere(r=rad*1e-2, x0=0.0, y0=0.0) for rad in [2.125, 3.125, 3.525, 3.875, 4.275]] geo_triso_cells = [openmc.Cell(name="TRISO Kernel UCO",fill=mat_triso_uco,region=-geo_triso_spheres[0]), openmc.Cell(name="TRISO Kernel Buffer",fill=mat_triso_buffer,region=+geo_triso_spheres[0]&-geo_triso_spheres[1]), openmc.Cell(name="TRISO Kernel PyC1",fill=mat_triso_pyc1,region=+geo_triso_spheres[1]&-geo_triso_spheres[2]), openmc.Cell(name="TRISO Kernel SiC",fill=mat_triso_sic,region=+geo_triso_spheres[2]&-geo_triso_spheres[3]), openmc.Cell(name="TRISO Kernel PyC2",fill=mat_triso_pyc2,region=+geo_triso_spheres[3]&-geo_triso_spheres[4])] 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_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 ) # Create the TRISO particles geo_triso_compact_trisos = [openmc.model.TRISO(geo_triso_spheres[4].r, geo_triso_universe, loc) for loc in geo_triso_compact_packing] # 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), shape=geo_triso_compact_domain, background=mat_graphite_matrix ) # Cast into cells and a universe geo_triso_compact_fuel_cell = openmc.Cell(name="TRISO Compact Cell",fill=geo_triso_compact_lattice,region=geo_triso_compact_region) geo_triso_compact_surrounding_cell = openmc.Cell(name="TRISO Compact Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_compact_region) geo_triso_compact_universe = openmc.Universe(name="TRISO Compact Universe",cells=[geo_triso_compact_fuel_cell,geo_triso_compact_surrounding_cell]) ## 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_n_triso = len(geo_triso_compact_packing) # 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_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_matrix = geo_triso_compact_rpt_v_total - geo_n_triso * 4/3*math.pi*geo_triso_spheres[4].r**3 geo_v_layers = [] geo_temp_prev_volume = 0.0 for layer in geo_triso_spheres: geo_temp_cumulative_volume = 4/3*math.pi*layer.r**3 geo_temp_shell_volume = geo_temp_cumulative_volume-geo_temp_prev_volume geo_v_layers.append(geo_temp_shell_volume*geo_n_triso) geo_temp_prev_volume=geo_temp_cumulative_volume geo_kernel_fracs = [v / sum(geo_v_layers) for v in geo_v_layers] geo_vwh_fracs = [v / geo_triso_compact_vwh_v_total for v in (geo_v_layers + [geo_triso_compact_vwh_v_matrix])] geo_rpt_fracs = [v / geo_triso_compact_rpt_v_total for v in (geo_v_layers + [geo_triso_compact_rpt_v_matrix])] # Strip TSL temporarily for mixing mat_tsl_saved = {} for mat in mat_list: if hasattr(mat, "_sab") and mat._sab: mat_tsl_saved[mat]=list(mat._sab) mat._sab=[] # TRISO only VWH mat_kernel_vwh = openmc.Material.mix_materials( materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2], fracs=geo_kernel_fracs, percent_type="vo", name="VWH kernel" ) mat_kernel_vwh.add_s_alpha_beta("c_Graphite") mat_list.append(mat_kernel_vwh) # Simple VWH mat_triso_vwh = openmc.Material.mix_materials( materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2, mat_graphite_matrix], fracs=geo_vwh_fracs, percent_type="vo", name="VWH Homogenized TRISO Compact" ) mat_triso_vwh.add_s_alpha_beta("c_Graphite") mat_list.append(mat_triso_vwh) # RPT Homogenization if cli_args.technique=="rpt": # There is a minimum viable radius for RPT, so only calculate this if needed to aovid errors mat_triso_rpt = openmc.Material.mix_materials( materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2, mat_graphite_matrix], fracs=geo_rpt_fracs, percent_type="vo", name="RPT Homogenized TRISO Compact" ) mat_triso_rpt.add_s_alpha_beta("c_Graphite") else: mat_triso_rpt = openmc.Material(name="PLACEHOLDER//DO NOT USE") mat_triso_rpt.add_nuclide("U235",1.0) mat_triso_rpt.set_density("g/cc",100.0) mat_list.append(mat_triso_rpt) # Restore TSL for other materials 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_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 # 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) geo_triso_vwh_surrounding_cell = openmc.Cell(name="TRISO VWH Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_compact_region) geo_triso_vwh_universe = openmc.Universe(name="TRISO VWH Universe",cells=[geo_triso_vwh_fuel_cell,geo_triso_vwh_surrounding_cell]) geo_triso_rpt_fuel_cell = openmc.Cell(name="TRISO RPT Fuel Cell",fill=mat_triso_rpt,region=geo_triso_rpt_region) geo_triso_rpt_surrounding_cell = openmc.Cell(name="TRISO RPT Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_rpt_region) geo_triso_rpt_universe = openmc.Universe(name="TRISO RPT Universe",cells=[geo_triso_rpt_fuel_cell,geo_triso_rpt_surrounding_cell]) geo_triso_rrpt_fuel_cell = openmc.Cell(name="TRISO RRPT Fuel Cell",fill=mat_kernel_vwh,region=geo_rrpt_region) geo_triso_rrpt_surrouding_cell = openmc.Cell(name="TRISO RRPT Surrounding Cell",fill=mat_graphite_matrix,region=~geo_rrpt_region) geo_triso_rrpt_universe = openmc.Universe(name="TRISO RRPT Universe",cells=[geo_triso_rrpt_fuel_cell,geo_triso_rrpt_surrouding_cell]) # Dynamically pick the correct one geo_homo_map = { "none":geo_triso_compact_universe, "vwh":geo_triso_vwh_universe, "rpt":geo_triso_rpt_universe, "rrpt":geo_triso_rrpt_universe } geo_dyn_universe = geo_homo_map[cli_args.technique] ## Heat Pipe ## geo_pipe_cylinders = [openmc.ZCylinder(r=rad, x0=0.0, y0=0.0) for rad in [0.80, 0.90, 0.97, 1.05, 1.07]] geo_pipe_cells = [openmc.Cell(name="Heat Pipe K Gas",fill=mat_pipe_k_gas,region=-geo_pipe_cylinders[0]), openmc.Cell(name="Heat Pipe Wick",fill=mat_pipe_wick,region=+geo_pipe_cylinders[0]&-geo_pipe_cylinders[1]), openmc.Cell(name="Heat Pipe K Layer",fill=mat_pipe_k_liquid,region=+geo_pipe_cylinders[1]&-geo_pipe_cylinders[2]), openmc.Cell(name="Heat Pipe SS316 Envelope",fill=mat_pipe_ss316,region=+geo_pipe_cylinders[2]&-geo_pipe_cylinders[3]), openmc.Cell(name="Heat Pipe Helium Gap",fill=mat_pipe_helium,region=+geo_pipe_cylinders[3]&-geo_pipe_cylinders[4]), openmc.Cell(name="Heat Pipe Graphite Matrix",fill=mat_graphite_matrix,region=+geo_pipe_cylinders[4])] geo_pipe_universe = openmc.Universe(name="Heat Pipe Universe",cells=geo_pipe_cells) ## Moderator Pin ## geo_mod_cylinders = [openmc.ZCylinder(r=rad, x0=0.0, y0=0.0) for rad in [0.825, 0.875, 0.900, 0.920]] geo_mod_cells = [openmc.Cell(name="Moderator Pin YH2",fill=mat_mod_yh2,region=-geo_mod_cylinders[0]), openmc.Cell(name="Moderator Pin Inner Helium Gap",fill=mat_mod_helium,region=+geo_mod_cylinders[0]&-geo_mod_cylinders[1]), openmc.Cell(name="Moderator Pin SS316 Envelope",fill=mat_mod_ss316,region=+geo_mod_cylinders[1]&-geo_mod_cylinders[2]), openmc.Cell(name="Moderator Pin Outer Helium Gap",fill=mat_mod_helium,region=+geo_mod_cylinders[2]&-geo_mod_cylinders[3]), openmc.Cell(name="Moderator Pin Graphite Matrix",fill=mat_graphite_matrix,region=+geo_mod_cylinders[3])] geo_mod_universe = openmc.Universe(name="Moderator Pin Universe",cells=geo_mod_cells) ### Assembly-level Geometry ### geo_assembly_lattice = openmc.HexLattice() geo_assembly_lattice.center=(0,0) geo_assembly_lattice.pitch=(geo_constant_pin_pitch,) geo_assembly_lattice.outer=openmc.Universe(cells=[openmc.Cell(fill=mat_graphite_matrix)]) # Create the lattice structure geo_assembly_lattice_map = [] for i in range(0,7): if i%2 == 0: ring = [geo_pipe_universe]*(36-6*i) ring[1::2] = [geo_dyn_universe]*len(ring[1::2]) geo_assembly_lattice_map.append(ring) else: ring = [geo_mod_universe]*(36-6*i) ring[::2] = [geo_dyn_universe]*len(ring[::2]) geo_assembly_lattice_map.append(ring) geo_assembly_lattice_map[-1] = [geo_pipe_universe] 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, orientation="y", origin=(0.0,0.0), boundary_type="reflective" ) # Clip the lattice to the space and make it root 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_universe = openmc.Universe(cells=[geo_assembly_cell]) geometry = openmc.Geometry(geo_assembly_universe) ############ SETTINGS ############ settings = openmc.Settings() settings.particles = 10000 settings.batches = 400 settings.inactive = 80 # Source sampling set_source_pts = [] for loc in geo_triso_compact_packing: point = openmc.stats.Point(xyz=loc) set_source_pts.append(openmc.IndependentSource(space=point)) 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.dimension=(30,30,30) settings.entropy_mesh=set_entropy_mesh settings.temperature = {'method': "interpolation"} ## Export ## materials = openmc.Materials(mat_list) materials.export_to_xml() geometry.export_to_xml() settings.export_to_xml() # Run! openmc.run()