path: root/hpmr.py

#!/bin/python3
import openmc
import math
import argparse as ap
import hashlib as hl
import sys
import time

def id_file():
    path = sys.argv[0]

    with open(path, "rb") as file:
        contents = file.read()

    filehash = hl.md5(contents).hexdigest()

    print("="*60)
    print(f"Run date:       {time.ctime()}")
    print(f"Input file MD5: {filehash}")
    print("="*60)


############ COMMAND LINE ARGUMENTS ############ 
id_file()

cli_parser = ap.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="none",choices=["none","vwh","rpt","rrpt"],help="What homogenization technique to use (if any) [none/vwh/rpt/rrpt], default none")
cli_parser.add_argument("-r","--radius",default="0.5",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()

############ MATERIALS ############ 
# All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
# <https://www.tandfonline.com/doi/full/10.1080/00295639.2024.2375175>
# 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
# <https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-15870Rev2.pdf>
mat_list = []

## Matrix Graphite ##
mat_graphite_matrix = openmc.Material(name="Graphite Matrix (inside/outside)")
mat_graphite_matrix.add_element("C",1.0)
mat_graphite_matrix.add_nuclide("B10",0.3e-6)
mat_graphite_matrix.set_density("g/cc",1.806)
mat_list.append(mat_graphite_matrix)

## TRISO Pins ##
mat_triso_uco = openmc.Material(name="TRISO Pin UCO Kernel")
mat_triso_uco.add_element("U",1.0,enrichment=19.95,enrichment_type="wo")
mat_triso_uco.add_element("C",1.0)
mat_triso_uco.add_element("O",1.0)
mat_triso_uco.set_density("g/cc",10.744)
mat_triso_uco.add_s_alpha_beta("c_Graphite")
mat_list.append(mat_triso_uco)

mat_triso_buffer = openmc.Material(name="TRISO Pin Buffer")
mat_triso_buffer.add_element("C",1.0)
mat_triso_buffer.add_s_alpha_beta("c_Graphite")
mat_triso_buffer.set_density("g/cc",1.04)
mat_list.append(mat_triso_buffer)

mat_triso_pyc1 = openmc.Material(name="TRISO Pin PyC1")
mat_triso_pyc1.add_element("C",1.0)
mat_triso_pyc1.add_s_alpha_beta("c_Graphite")
mat_list.append(mat_triso_pyc1)

mat_triso_sic = openmc.Material(name="TRISO Pin SiC")
mat_triso_sic.add_element("Si",1.0)
mat_triso_sic.add_element("C",1.0)
mat_triso_sic.add_s_alpha_beta("c_Graphite")
mat_triso_sic.set_density("g/cc",3.171)
mat_list.append(mat_triso_sic)

mat_triso_pyc2 = openmc.Material(name="TRISO Pin PyC2")
mat_triso_pyc2.add_element("C",1.0)
mat_triso_pyc2.add_s_alpha_beta("c_Graphite")
mat_triso_pyc2.set_density("g/cc",1.882)
mat_list.append(mat_triso_pyc2)

## Moderator Pins ##
mat_mod_helium = openmc.Material(name="Moderator Pin Helium")
mat_mod_helium.add_element("He",1.0)
mat_mod_helium.set_density("g/cc",0.18e-3)
mat_list.append(mat_mod_helium)

# Composition from PNNL-15870r2
mat_mod_ss316 = openmc.Material(name="Moderator Pin SS316")
mat_mod_ss316.add_element("C",0.000800,"wo")
mat_mod_ss316.add_element("Mn",0.020000,"wo")
mat_mod_ss316.add_element("P",0.000450,"wo")
mat_mod_ss316.add_element("S",0.000300,"wo")
mat_mod_ss316.add_element("Si",0.010000,"wo")
mat_mod_ss316.add_element("Cr",0.170000,"wo")
mat_mod_ss316.add_element("Ni",0.120000,"wo")
mat_mod_ss316.add_element("Mo",0.025000,"wo")
mat_mod_ss316.add_element("Fe",0.653450,"wo")
mat_mod_ss316.set_density("g/cc",7.67)
mat_list.append(mat_mod_ss316)

mat_mod_yh2 = openmc.Material(name="Moderater Pin YH2")
mat_mod_yh2.add_element("Y",1.0)
mat_mod_yh2.add_element("H",2.0)
mat_mod_yh2.add_s_alpha_beta("c_H_in_YH2")
mat_mod_yh2.add_s_alpha_beta("c_Y_in_YH2")
mat_mod_yh2.set_density("g/cc",4.3)
mat_list.append(mat_mod_yh2)

## Heat Pipes ##
mat_pipe_helium = openmc.Material(name="Heat Pipe Helium")
mat_pipe_helium.add_element("He",1.0)
mat_pipe_helium.set_density("g/cc",0.18e-3)
mat_list.append(mat_pipe_helium)

# Composition from PNNL-15870r2
mat_pipe_ss316 = openmc.Material(name="Heat Pipe SS316")
mat_pipe_ss316.add_element("C",0.000800,"wo")
mat_pipe_ss316.add_element("Mn",0.020000,"wo")
mat_pipe_ss316.add_element("P",0.000450,"wo")
mat_pipe_ss316.add_element("S",0.000300,"wo")
mat_pipe_ss316.add_element("Si",0.010000,"wo")
mat_pipe_ss316.add_element("Cr",0.170000,"wo")
mat_pipe_ss316.add_element("Ni",0.120000,"wo")
mat_pipe_ss316.add_element("Mo",0.025000,"wo")
mat_pipe_ss316.add_element("Fe",0.653450,"wo")
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])

## VWH TRISO Compact ##
geo_vwh_fracs = []
geo_n_triso = len(geo_triso_compact_packing)

# Get volume fractions for mixing
geo_triso_compact_v_tot = math.pi*geo_constant_triso_radius**2*2*geo_constant_height

geo_triso_compact_v_uco = geo_n_triso*4/3*math.pi*geo_triso_spheres[0].r**3
geo_vwh_fracs.append(geo_triso_compact_v_uco/geo_triso_compact_v_tot)

geo_triso_compact_v_buffer = geo_n_triso*4/3*math.pi*geo_triso_spheres[1].r**3 - geo_triso_compact_v_uco
geo_vwh_fracs.append(geo_triso_compact_v_buffer/geo_triso_compact_v_tot)

geo_triso_compact_v_pyc1 = geo_n_triso*4/3*math.pi*geo_triso_spheres[2].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer)
geo_vwh_fracs.append(geo_triso_compact_v_pyc1/geo_triso_compact_v_tot)

geo_triso_compact_v_sic = geo_n_triso*4/3*math.pi*geo_triso_spheres[3].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1)
geo_vwh_fracs.append(geo_triso_compact_v_sic/geo_triso_compact_v_tot)

geo_triso_compact_v_pyc2 = geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1 + geo_triso_compact_v_sic)
geo_vwh_fracs.append(geo_triso_compact_v_pyc2/geo_triso_compact_v_tot)

geo_triso_compact_v_graphite = geo_triso_compact_v_tot - geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3
geo_vwh_fracs.append(geo_triso_compact_v_graphite/geo_triso_compact_v_tot)

# 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=[]

# 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"
        )

# RPT Homogenization

# Restore TSL afterwards
for mat, tsl_list in mat_tsl_saved.items():
    mat._sab = tsl_list

mat_triso_vwh.add_s_alpha_beta("c_Graphite")

# Add to list
mat_list.append(mat_triso_vwh)
materials = openmc.Materials(mat_list)

# Cast into cells and a universe
geo_triso_vwh_fuel_cell = openmc.Cell(name="TRISO VWH 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])

## 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_triso_vwh_universe]*len(ring[1::2])
        geo_assembly_lattice_map.append(ring)
    else:
        ring = [geo_mod_universe]*(36-6*i)
        ring[::2] = [geo_triso_vwh_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 = 100
settings.inactive = 50

# 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

## Export ##
materials.export_to_xml()
geometry.export_to_xml()
settings.export_to_xml()

# Run!
openmc.run()