import openmc import numpy as np import argparse as ap #> Standard 37-element CANDU fuel bundle in OpenMC for AnnCon2026. This will #> serve as a basis for a comparison with the Apollonian gasket verison. #> Connor Moore, March 2026. ### Argparser for parametric analysis ### parser = ap.ArgumentParser(prog="37-Element Bundle Analysis", description="Variable pin radius for calculating ratios and criticality.", epilog="Have fun and be yourself!") parser.add_argument("-rf","--radius-fuel",help="Radius of the fuel pin [cm], rf ∈ (0.05,0.8] cm. Default is 0.6 cm.",default="0.6") 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 ### #> The general geometry of the 37-element bundle is from IGE335 by #> Dr. Alain Hebert and contributors, part of the DRAGON project. #> fuel_region_list = [] clad_region_list = [] #> From the Essential CANDU r_fuel = eval(args.radius_fuel) r_clad = r_fuel + 0.04 #> Define a function to create a fuel "pin" using a radius. def make_pin(rf: float, rc: float, x0: float, y0: float) -> None: fuel_surf = openmc.ZCylinder(r=rf, x0=x0, y0=y0) #> Hard-coded cladding thickness of 0.4 mm clad_surf = openmc.ZCylinder(r=rc, 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) return #> Also define a function for making rings of pins def make_ring(n: int, r: float, alpha: float) -> None: for i in range(0,n): theta = (i*(360/n) + alpha)*np.pi/180 x = r*np.cos(theta) y = r*np.sin(theta) #print((float(x),float(y))) make_pin(rf=r_fuel, rc=r_clad, x0=x, y0=y) return #> Start with the innermost pin make_pin(rf=r_fuel, rc=r_clad, x0=0.0, y0=0.0) #> Add a ring of 6 elements around it. The radius is 1.4885 with an angle of 0.0 make_ring(n=6, r=1.4885, alpha=0.0) #> The next ring has 12 pins at a radius of 2.8755, with an angle of 15 degrees make_ring(n=12, r=2.8755, alpha=15.0) #> And the last one has 18 pins, a radius of 4.3305, and an angle of 0.0 make_ring(n=18, r=4.3305, alpha=0.0) #> 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 = 200 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 = len(fuel_region_list)*np.pi*r_fuel**2 V_clad = len(clad_region_list)*np.pi*(r_clad**2 - r_fuel**2) C_clad = len(clad_region_list)*2*np.pi*r_clad 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 ciricumference, 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()