#!/usr/bin/python3 import openmc import math import argparse import hashlib import sys import time # KRUSTY Model in OpenMC # C.Moore, ############ COMMAND LINE ARGUMENTS ############ cli_parser = argparse.ArgumentParser( prog="KRUSTY Heat Pipe Study", description="Program to investigate the reactivity effect of heat pipe saturation on the KRUSTY reactor." ) #cli_parser.add_argument("-t","--technique",default="none",choices=["none","vwh","rpt","rrpt"],help="What homogenization technique to use (if any) [none/vwh/rpt/rrpt], default none") cli_parser.add_argument("-n","--na-height",default="0.8",type=float,help="Radius for RPT or RRPT homogenization [cm], default 0.8 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 ############ mat_list = [] mat_temp_global = 900 mat_temp_fuel = 900 ## 100 - Structural Materials ## mat_haynes230 = openmc.Material(100,name='Structural HAYNES 230') mat_haynes230.add_element('Ni',57,'wo') mat_haynes230.add_element('Cr',22,'wo') mat_haynes230.add_element('W',14,'wo') mat_haynes230.add_element('Mo', 2,'wo') mat_haynes230.add_element('Fe', 1.875,'wo') mat_haynes230.add_element('Co', 3.125,'wo') mat_haynes230.temperature = mat_temp_global mat_haynes230.set_density('g/cm3',8.97) mat_list.append(mat_haynes230) mat_ss316 = openmc.Material(101,name='Structural SS316') mat_ss316.add_element('C', 9.12187E-02+1.01457E-03, 'ao') mat_ss316.add_nuclide('Si28', 1.10957E+00, 'ao') mat_ss316.add_nuclide('Si29', 5.61820E-02, 'ao') mat_ss316.add_nuclide('Si30', 3.72944E-02, 'ao') mat_ss316.add_nuclide('Cr50', 7.94742E-01, 'ao') mat_ss316.add_nuclide('Cr52', 1.53260E+01, 'ao') mat_ss316.add_nuclide('Cr53', 1.73763E+00, 'ao') mat_ss316.add_nuclide('Cr54', 4.32580E-01, 'ao') mat_ss316.add_nuclide('Mn55', 1.77454E+00, 'ao') mat_ss316.add_nuclide('Fe54', 4.11552E+00, 'ao') mat_ss316.add_nuclide('Fe56', 6.39783E+01, 'ao') mat_ss316.add_nuclide('Fe57', 1.46483E+00, 'ao') mat_ss316.add_nuclide('Fe58', 1.95315E-01, 'ao') mat_ss316.add_nuclide('Ni58', 5.91824E+00, 'ao') mat_ss316.add_nuclide('Ni60', 2.26262E+00, 'ao') mat_ss316.add_nuclide('Ni61', 9.79589E-02, 'ao') mat_ss316.add_nuclide('Ni62', 3.11211E-01, 'ao') mat_ss316.add_nuclide('Ni64', 7.88874E-02, 'ao') mat_ss316.add_nuclide('Mo92', 3.21318E-02, 'ao') mat_ss316.add_nuclide('Mo94', 2.00286E-02, 'ao') mat_ss316.add_nuclide('Mo95', 3.44709E-02, 'ao') mat_ss316.add_nuclide('Mo96', 3.61161E-02, 'ao') mat_ss316.add_nuclide('Mo97', 2.06783E-02, 'ao') mat_ss316.add_nuclide('Mo98', 5.22477E-02, 'ao') mat_ss316.add_nuclide('Mo100', 2.08518E-02, 'ao') mat_ss316.temperature = mat_temp_global mat_ss316.set_density('atom/b-cm', 8.58870E-02) mat_list.append(mat_ss316) mat_air = openmc.Material(102,name="Core Air") mat_air.add_nuclide('H1', 1.41615E+01, 'ao') mat_air.add_nuclide('H2', 2.12456E-03, 'ao') mat_air.add_element('C', 1.17053E-02+1.30189E-04, 'ao') mat_air.add_nuclide('N14', 6.15481E+01, 'ao') mat_air.add_nuclide('N15', 2.28569E-01, 'ao') mat_air.add_nuclide('O16', 2.36230E+01+4.73602E-02, 'ao') mat_air.add_nuclide('O17', 9.47204E-03, 'ao') mat_air.add_nuclide('Ar36', 1.23939E-03, 'ao') mat_air.add_nuclide('Ar38', 2.31686E-04, 'ao') mat_air.add_nuclide('Ar40', 3.66299E-01, 'ao') mat_air.temperature = mat_temp_global mat_air.set_density('g/cm3', 0.5) mat_list.append(mat_air) ## 200 - Core Materials ## mat_fuel = openmc.Material(200,name='Fuel U-10Mo') mat_fuel.add_element('C', 6.19354E-01+6.88872E-03,'ao') mat_fuel.add_nuclide('Mo92', 2.49068E+00,'ao') mat_fuel.add_nuclide('Mo94', 1.55250E+00,'ao') mat_fuel.add_nuclide('Mo95', 2.67188E+00,'ao') mat_fuel.add_nuclide('Mo96', 2.79947E+00,'ao') mat_fuel.add_nuclide('Mo97', 1.60285E+00,'ao') mat_fuel.add_nuclide('Mo98', 4.04995E+00,'ao') mat_fuel.add_nuclide('Mo100', 1.61627E+00,'ao') mat_fuel.add_nuclide('U234', 8.45363E-01,'ao') mat_fuel.add_nuclide('U235', 7.69214E+01,'ao') mat_fuel.add_nuclide('U236', 3.81533E-01,'ao') mat_fuel.add_nuclide('U238', 4.44178E+00,'ao') mat_fuel.temperature = mat_temp_fuel mat_fuel.set_density('atom/b-cm', 4.96130E-02) mat_list.append(mat_fuel) mat_hp_na = openmc.Material(201,name="Heat Pipe Sodium Coolant") mat_hp_na.add_element("Na",1) mat_hp_na.temperature = mat_temp_global mat_hp_na.set_density("g/cc",0.763) mat_list.append(mat_hp_na) mat_beo = openmc.Material(202,name="BeO Reflector") mat_beo.add_nuclide('Be9', 4.99528E+01,'ao') mat_beo.add_nuclide('O16', 4.99042E+01+1.00046E-01,'ao') mat_beo.add_nuclide('O17', 2.00099E-02,'ao') mat_beo.add_nuclide('S32', 2.24316E-02,'ao') mat_beo.add_nuclide('S33', 1.77054E-04,'ao') mat_beo.add_nuclide('S34', 9.93896E-04,'ao') mat_beo.add_nuclide('S36', 4.72141E-06,'ao') mat_beo.temperature = mat_temp_global mat_beo.set_density('atom/b-cm', 0.135648) mat_beo.add_s_alpha_beta("c_Be_in_BeO") mat_beo.add_s_alpha_beta("c_O_in_BeO") mat_list.append(mat_beo) mat_b4c = openmc.Material(203,name="B4C Control Rods") mat_b4c.add_nuclide('B10',4,'ao') mat_b4c.add_element('C',1,'ao') mat_b4c.temperature = mat_temp_global mat_b4c.set_density('g/cm3',2.52) mat_list.append(mat_b4c) materials = openmc.Materials(mat_list) materials.export_to_xml() ############ GEOMETRY ############ geo_hp_pitch = 5.2 geo_z_abs_bot = openmc.ZPlane(z0=0.0,boundary_type="vacuum") geo_z_core_bot = openmc.ZPlane(z0=10.16) geo_z_core_top = openmc.ZPlane(z0=35.16) geo_z_abs_top = openmc.ZPlane(z0=45.32,boundary_type="vacuum") geo_z_region = +geo_z_abs_bot & -geo_z_abs_top ## Central Air Cylinder ## geo_air_cyl = openmc.ZCylinder(r=1.9939, x0=0.0, y0=0.0) geo_air_region = -geo_air_cyl & +geo_z_core_bot geo_air_cell = openmc.Cell(name="Inner Core Air", region=geo_air_region, fill=mat_air) ## Fuel ## geo_fuel_cyl = openmc.ZCylinder(r=5.4991, x0=0.0, y0=0.0) geo_fuel_region = -geo_fuel_cyl & +geo_air_cyl & -geo_z_core_top & +geo_z_core_bot geo_fuel_cell = openmc.Cell(name="Fuel Block", region=geo_fuel_region, fill=mat_fuel) ## BeO Inner Reflectors ## # Bottom Reflector geo_beo_lower_region = -geo_fuel_cyl & -geo_z_core_bot geo_beo_lower_cell = openmc.Cell(name="Lower BeO Reflector", region=geo_beo_lower_region, fill=mat_beo) # Top reflector geo_beo_upper_region = +geo_air_cyl & -geo_fuel_cyl & +geo_z_core_top geo_beo_upper_cell = openmc.Cell(name="Upper BeO Reflector", region=geo_beo_upper_region, fill=mat_beo) ## Radial Layers ## # Air, SS316, Air, BeO geo_rad_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [12.7, 13.3, 14.11, 38.40]] geo_rad_cyl[3].boundary_type="vacuum" geo_rad_cells = [openmc.Cell(name="Intermediate Core Air",region=+geo_fuel_cyl&-geo_rad_cyl[0],fill=mat_air), openmc.Cell(name="Core SS316 Cell",region=+geo_rad_cyl[0]&-geo_rad_cyl[1],fill=mat_ss316), openmc.Cell(name="Outer Core Air",region=+geo_rad_cyl[1]&-geo_rad_cyl[2],fill=mat_air), openmc.Cell(name="Outer BeO Reflector",region=+geo_rad_cyl[2],fill=mat_beo)] # Export simplified universe without heat pipes geo_rad_universe = openmc.Universe(name="Simplified Radial Universe",cells=[geo_air_cell,geo_fuel_cell,geo_beo_lower_cell,geo_beo_upper_cell]+geo_rad_cells) ## Heat Pipes ## geo_hp_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [1.092, 1.270]] geo_hp_unit_cell = [openmc.Cell(name="Heat Pipe Coolant Cell",region=-geo_hp_cyl[0],fill=mat_hp_na), openmc.Cell(name="Heat Pipe Wall Cell",region=+geo_hp_cyl[0],fill=mat_haynes230)] #openmc.Cell(name="Heat Pipe Surrounding Cell",region=+geo_hp_cyl[1],fill=mat_air)] geo_hp_unit_universe = openmc.Universe(name="Heat Pipe Universe",cells=geo_hp_unit_cell) # Pattern around the core geo_hp_steps = [math.pi*theta/180 for theta in range(0,360,45)] geo_hp_regions = [openmc.ZCylinder(r=geo_hp_cyl[1].r, x0=geo_hp_pitch*math.cos(theta), y0=geo_hp_pitch*math.sin(theta) ) for theta in geo_hp_steps] geo_hp_cells = [openmc.Cell(name=f"Heat Pipe {num}",region=-space,fill=geo_hp_unit_universe) for num,space in enumerate(geo_hp_regions, start=1)] for cell, region in zip(geo_hp_cells, geo_hp_regions): cell.translation = (region.x0, region.y0, 0.0) geo_hp_universe = openmc.Universe(name="Heat Pipes Universe",cells=geo_hp_cells) # Combine universes geo_hp_comb_cell = openmc.Cell(name="Combined Heat Pipes Cell",region=openmc.Union(-cyl for cyl in geo_hp_regions),fill=geo_hp_universe) geo_rad_comb_cell = openmc.Cell(name="Combined Radial Layers Cell",region=~openmc.Union(-cyl for cyl in geo_hp_regions),fill=geo_rad_universe) geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_rad_comb_cell,geo_hp_comb_cell]) geo_clip_cell = openmc.Cell(name="Clipped Root Universe Cell",region=-geo_rad_cyl[3],fill=geo_root_universe) geo_clip_universe = openmc.Universe(name="Clipped Root Universe",cells=[geo_clip_cell]) geo_object = openmc.Geometry(geo_clip_universe) geo_object.export_to_xml() ## Settings ## set_obj = openmc.Settings() set_obj.particles = 10000 set_obj.batches = 1500 set_obj.inactive = 100 set_obj.temperature = {"method": "interpolation"} set_src = openmc.IndependentSource() set_src.space = openmc.stats.CylindricalIndependent( r = openmc.stats.Uniform(geo_air_cyl.r,geo_fuel_cyl.r), z = openmc.stats.Uniform(geo_z_core_bot.z0,geo_z_core_top.z0), phi = openmc.stats.Uniform(0.0,2*math.pi), origin = (0.0,0.0,geo_z_core_bot.z0) ) set_src.track = [(1,1,i) for i in range(1,20)] set_src.max_events = 2000 set_obj.export_to_xml() # Run! openmc.run()