diff options
Diffstat (limited to 'krusty.py')
| -rwxr-xr-x | krusty.py | 243 |
1 files changed, 243 insertions, 0 deletions
diff --git a/krusty.py b/krusty.py new file mode 100755 index 0000000..f010689 --- /dev/null +++ b/krusty.py @@ -0,0 +1,243 @@ +#!/usr/bin/python3 +import openmc +import math +import argparse +import hashlib +import sys +import time + +# KRUSTY Model in OpenMC +# C.Moore, <connor.moore@psi.ch> + +############ 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() |