Added cladding circumference, 52-element (and alt 51-element) Apollonian bunddles

1 files changed, 197 insertions, 0 deletions

diff --git a/apollonian-52.py b/apollonian-52.py
new file mode 100644
index 0000000..2b7fcd8
--- /dev/null
+++ b/apollonian-52.py
@@ -0,0 +1,197 @@
+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. <connor.moore@ontariotechu.net>
+
+### 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.
+#> <https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-15870Rev2.pdf>
+
+#> 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 = 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 = 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 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()