#!/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=-10.16,boundary_type="vacuum") geo_z_core_bot = openmc.ZPlane(z0=0.00) geo_z_core_top = openmc.ZPlane(z0=25) geo_z_abs_top = openmc.ZPlane(z0=35.16,boundary_type="vacuum") geo_z_region = +geo_z_abs_bot & -geo_z_abs_top ## Heat Pipes ## geo_hp_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [1.092, 1.270]] # Pattern around the core geo_hp_steps = [math.pi*(theta+90)/180 for theta in range(0,360,45)] # Define inner and outer regions geo_hp_inner_regions = [openmc.ZCylinder(r=geo_hp_cyl[0].r, x0=geo_hp_pitch*math.cos(theta), y0=geo_hp_pitch*math.sin(theta) ) for theta in geo_hp_steps] geo_hp_outer_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_outer_union = openmc.Union(-reg for reg in geo_hp_outer_regions) # Fill cells for each pipe geo_hp_inner_cells = [openmc.Cell(name=f"Heat Pipe {num} Inner Cell",region=-reg&geo_z_region,fill=mat_hp_na) for num,reg in enumerate(geo_hp_inner_regions)] geo_hp_outer_cells = [openmc.Cell(name=f"Heat Pipe {num} Outer Cell",region=-reg_out&+reg_in&geo_z_region,fill=mat_haynes230) for num,(reg_out,reg_in) in enumerate(zip(geo_hp_outer_regions,geo_hp_inner_regions))] ## Core Air Cylinder ## geo_air_inner_cyl = openmc.ZCylinder(r=1.9939, x0=0.0, y0=0.0) geo_air_inner_region = +geo_z_core_bot & -geo_z_abs_top & -geo_air_inner_cyl geo_air_inner_cell = openmc.Cell(name="Inner Air Cell",region=geo_air_inner_region,fill=mat_air) ## Core Fuel Block ## geo_fuel_cyl = openmc.ZCylinder(r=5.4991, x0=0.0, y0=0.0) geo_fuel_region = +geo_z_core_bot & -geo_z_core_top & -geo_fuel_cyl & +geo_air_inner_cyl & ~geo_hp_outer_union #geo_fuel_region = +geo_z_abs_bot & -geo_z_abs_top & -geo_fuel_cyl geo_fuel_cell = openmc.Cell(name="Core Fuel Cell",region=geo_fuel_region,fill=mat_fuel) ## BeO Reflectors ## # Lower Reflector geo_refl_lower_region = +geo_z_abs_bot & -geo_z_core_bot & -geo_fuel_cyl & ~geo_hp_outer_union geo_refl_lower_cell = openmc.Cell(name="Lower BeO Reflector Cell",region=geo_refl_lower_region,fill=mat_beo) # Upper Reflector geo_refl_upper_region = +geo_z_core_top & -geo_z_abs_top & +geo_air_inner_cyl & -geo_fuel_cyl & ~geo_hp_outer_union geo_refl_upper_cell = openmc.Cell(name="Upper BeO Reflector Cell",region=geo_refl_upper_region,fill=mat_beo) ## Radial Core Surroundings ## geo_rad_cyl = [openmc.ZCylinder(r=dia/2, x0=0.0, y0=0.0) for dia in [12.70, 13.30, 14.11, 38.40]] geo_rad_cyl[-1].boundary_type="vacuum" geo_rad_cells = [openmc.Cell(name="Intermediate Air Cell",region=+geo_fuel_cyl&-geo_rad_cyl[0]&geo_z_region,fill=mat_air), openmc.Cell(name="Steel Annulus Cell",region=+geo_rad_cyl[0]&-geo_rad_cyl[1]&geo_z_region,fill=mat_ss316), openmc.Cell(name="Outer Air Cell",region=+geo_rad_cyl[1]&-geo_rad_cyl[2]&geo_z_region,fill=mat_air), openmc.Cell(name="Outer BeO Reflector Cell",region=+geo_rad_cyl[2]&-geo_rad_cyl[3]&geo_z_region,fill=mat_beo)] geo_root_universe = openmc.Universe(name="Root Universe",cells=[geo_fuel_cell,geo_air_inner_cell,geo_refl_upper_cell,geo_refl_lower_cell]+geo_rad_cells+geo_hp_inner_cells+geo_hp_outer_cells) geo_obj = openmc.Geometry(geo_root_universe) geo_obj.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.Point((0.0,0.0,30.0)) set_src.space = openmc.stats.CylindricalIndependent( r = openmc.stats.Uniform(geo_air_inner_cyl.r,geo_fuel_cyl.r-1.0), #r = openmc.stats.Uniform(0,2), 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,0.0) ) #set_src.space = openmc.stats.Box(lower_left=(-2,-2,-2),upper_right=(2,2,6)) ## DEBUG SOURCE #for i in range(10000): # sample_r = set_src.space.r.sample() # sample_z = set_src.space.z.sample() # sample_phi = set_src.space.phi.sample() # #print(f"{sample_r = }, {sample_z = }, {sample_phi = }") # x = sample_r*math.cos(sample_phi[0]) # y = sample_r*math.sin(sample_phi[0]) # # cell = geo_obj.find((x,y,sample_z))[-1] # print(f"Particle {i}: ({x}, {y}, {sample_z}) -> Cell: {cell}") # #set_obj.verbosity=10 set_obj.export_to_xml() # Run! openmc.run()