import openmc import numpy as np import argparse as ap import csv #> Apollonian 52-element CANDU fuel bundle in OpenMC for AnnCon2026. #> Connor Moore, March 2026. ### Argparser for parametric analysis ### parser = ap.ArgumentParser(prog="52-Element Apollonian Bundle Analysis", description="Variable pin radius for calculating ratios and criticality.") parser.add_argument("-fmr","--fuel-moderator-ratio",help="Ratio of the radius of fuel to whole circle. Default is 0.75",default="0.75") parser.add_argument("-bt","--bundle-type",help="Apollonian bundle type [normal/alt], default=normal",default="normal") args = parser.parse_args() ### Material Definitions ### #> Materials marked PNNL-15870 are from the 2021 revision of "Compendium of Material #> Composition Data for Radiation Transport Modeling" published by the U.S. Department #> of Homeland Security and Pacific Northwest National Laboratory. #> #> The natural UO2 fuel mat_fuel = openmc.Material(name="Natural Uranium Fuel (UO2)") mat_fuel.add_element("U", 1.0, enrichment=0.71) mat_fuel.add_element("O", 2.0) mat_fuel.set_density("g/cc", 10.6) #> Density of 10.6 g/cc is from The Essential CANDU, Ch. 17 (Fuel), sec. 2.2, pp. 11 #> Pressure tube and calandria tube use Zircaloy-2 mat_zircaloy_2 = openmc.Material(name="Zircaloy-2 (PNNL-15870)") mat_zircaloy_2.add_nuclide("O16", 0.001194,"wo") mat_zircaloy_2.add_nuclide("O17", 0.000000,"wo") mat_zircaloy_2.add_nuclide("O18", 0.000003,"wo") mat_zircaloy_2.add_element("Cr", 0.000997,"wo") mat_zircaloy_2.add_element("Fe", 0.000997,"wo") mat_zircaloy_2.add_element("Ni", 0.000499,"wo") mat_zircaloy_2.add_nuclide("Zr90", 0.498109,"wo") mat_zircaloy_2.add_nuclide("Zr91", 0.109835,"wo") mat_zircaloy_2.add_nuclide("Zr92", 0.169730,"wo") mat_zircaloy_2.add_nuclide("Zr94", 0.175752,"wo") mat_zircaloy_2.add_nuclide("Zr96", 0.028918,"wo") mat_zircaloy_2.add_element("Sn", 0.013962,"wo") mat_zircaloy_2.set_density("g/cc",6.56) #> The fuel cladding uses Zircaloy-4 mat_zircaloy_4 = openmc.Material(name="Zircaloy-4 (PNNL-15870") mat_zircaloy_4.add_nuclide("O16", 0.001193,"wo") mat_zircaloy_4.add_nuclide("O17", 0.000000,"wo") mat_zircaloy_4.add_nuclide("O18", 0.000003,"wo") mat_zircaloy_4.add_element("Cr", 0.000997,"wo") mat_zircaloy_4.add_element("Fe", 0.001993,"wo") mat_zircaloy_4.add_nuclide("Zr90", 0.497860,"wo") mat_zircaloy_4.add_nuclide("Zr91", 0.109780,"wo") mat_zircaloy_4.add_nuclide("Zr92", 0.169646,"wo") mat_zircaloy_4.add_nuclide("Zr94", 0.175665,"wo") mat_zircaloy_4.add_nuclide("Zr96", 0.028904,"wo") mat_zircaloy_4.add_element("Sn", 0.013955,"wo") mat_zircaloy_4.set_density("g/cc",6.56) #> Heavy water is used for moderation and cooling mat_d2o = openmc.Material(name="Heavy Water (PNNL-15870)") mat_d2o.add_nuclide("H2", 0.201133,"wo") mat_d2o.add_nuclide("O16", 0.796703,"wo") mat_d2o.add_nuclide("O17", 0.000323,"wo") mat_d2o.add_nuclide("O18", 0.001842,"wo") mat_d2o.add_s_alpha_beta("c_D_in_D2O") mat_d2o.set_density("g/cc",1.1044) materials = openmc.Materials([mat_fuel, mat_zircaloy_2, mat_zircaloy_4, mat_d2o]) materials.export_to_xml() ### Geometry Definition ### fuel_region_list = [] clad_region_list = [] fuel_rad_list = [] fuel_area_list = [] clad_area_list = [] clad_circumference_list = [] #> Define a function to create a fuel "pin" using a radius. def make_apollonian_pin(rg: float, rfm: float, x0: float, y0: float) -> None: fuel_surf = openmc.ZCylinder(r=rg*rfm, x0=x0, y0=y0) #> Hard-coded cladding thickness of 0.4 mm clad_surf = openmc.ZCylinder(r=fuel_surf.r+0.04, x0=x0, y0=y0) fuel_region = -fuel_surf clad_region = +fuel_surf & -clad_surf fuel_region_list.append(fuel_region) clad_region_list.append(clad_region) fuel_area_list.append(np.pi*fuel_surf.r**2) clad_area_list.append(np.pi*(clad_surf.r**2 - fuel_surf.r**2)) clad_circumference_list.append(2*np.pi*clad_surf.r) return #> Import gasket points from file filename = "gasket-geometry/apollonian-gasket-52-normal.csv" if args.bundle_type=="normal" else "gasket-geometry/apollonian-gasket-52-alt.csv" rfm = eval(args.fuel_moderator_ratio) with open(filename, "r", newline="", encoding="utf-8") as file: reader = csv.reader(file) try: header = next(reader) except StopIteration: header = [] for row in reader: make_apollonian_pin(float(row[2]),rfm,float(row[0]),float(row[1])) #> Combine the regions to make a fuel cell and a cladding cell fuel_cell = openmc.Cell(name="UO2 Fuel Regions", region=openmc.Union(fuel_region_list), fill=mat_fuel) clad_cell = openmc.Cell(name="Zircaloy-4 Cladding Regions", region=openmc.Union(clad_region_list), fill=mat_zircaloy_4) #> Add the pressure tube pt_inner = openmc.ZCylinder(r=5.16890, x0=0.0, y0=0.0) pt_outer = openmc.ZCylinder(r=5.60320, x0=0.0, y0=0.0) pt_region = +pt_inner & -pt_outer pt_cell = openmc.Cell(name="Pressure Tube", region=pt_region, fill=mat_zircaloy_2) #> Pack the fuel with D2O pt_d2o_region = ~(fuel_cell.region | clad_cell.region) & -pt_inner pt_d2o_cell = openmc.Cell(name="D2O Coolant", region=pt_d2o_region, fill=mat_d2o) #> Add the calandria tube ct_inner = openmc.ZCylinder(r=6.44780, x0=0.0, y0=0.0) ct_outer = openmc.ZCylinder(r=6.58750, x0=0.0, y0=0.0) ct_region = +ct_inner & -ct_outer ct_cell = openmc.Cell(name="Calandria Tube", region=ct_region, fill=mat_zircaloy_2) #> The space between the calandria tube and pressure tube is considered void pt_ct_region = +pt_outer & -ct_inner pt_ct_cell = openmc.Cell(name="Annulus Gap", region=pt_ct_region, fill=None) #> The lattice pitch for the assembly is 28.575, so apply a reflecting boundary outer_boundary = openmc.model.RectangularPrism(width=28.575, height=28.575, origin=(0.0, 0.0), boundary_type="reflective") ct_d2o_region = +ct_outer & -outer_boundary ct_d2o_cell = openmc.Cell(name="D2O Moderator", region=ct_d2o_region, fill=mat_d2o) universe = openmc.Universe(cells=[fuel_cell, clad_cell, pt_cell, pt_d2o_cell, ct_cell, pt_ct_cell, ct_d2o_cell]) geometry = openmc.Geometry(universe) geometry.export_to_xml() ### Settings definition ### settings = openmc.Settings() settings.particles = 10000 settings.batches = 400 settings.inactive = 80 #> Set up the source to sample inside the pressure tube region uniformly settings.source = openmc.IndependentSource() settings.source.space = openmc.stats.CylindricalIndependent( r = openmc.stats.PowerLaw(a=0, b=pt_inner.r, n=1), phi = openmc.stats.Uniform(0, 2*np.pi), z = openmc.stats.Discrete([0.0], [1.0]) ) settings.source.energy = openmc.stats.Watt() settings.export_to_xml() ### Area and DTU Calculations ### V_fuel = sum(fuel_area_list) V_clad = sum(clad_area_list) C_clad = sum(clad_circumference_list) V_mod = np.pi*pt_inner.r**2 - (V_fuel + V_clad) + (28.575**2 - np.pi*ct_outer.r**2) d_fuel = mat_fuel.get_nuclide_atom_densities() N_fuel = sum(density for nuclide,density in d_fuel.items() if "U" in nuclide) N_mod = mat_d2o.get_nuclide_atom_densities()["H2"] DTU_ratio = (V_mod * N_mod) / (V_fuel * N_fuel) ##> Print these for the final table print(f"Flow area, cm² = {V_mod}") print(f"Cladding circumference, cm = {C_clad}") print(f"Fuel mass per length, g/cm = {V_fuel*mat_fuel.density}") print(f"Cladding mass per length, g/cm = {V_clad*mat_zircaloy_4.density}") print(f"DTU ratio = {DTU_ratio}") #> Run the model! openmc.run()