Added parametric fuel pin size and driving script for 37-element CANDU bundle

1 files changed, 17 insertions, 8 deletions

diff --git a/candu-37.py b/candu-37.py
index 08a90e1..e9642e7 100644
--- a/candu-37.py
+++ b/candu-37.py
@@ -1,10 +1,19 @@
 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. <connor.moore@ontariotechu.net>
 
+### 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 ###
 
@@ -18,7 +27,7 @@ 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), 2.2, pp. 11
+#> 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)")
@@ -73,8 +82,9 @@ materials.export_to_xml()
 fuel_region_list = []
 clad_region_list = []
 
-r_fuel = 0.60
-r_clad = 0.64
+#> 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.
 
@@ -161,7 +171,7 @@ 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=5.16890, n=1),
+        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])
         )
@@ -174,7 +184,7 @@ settings.export_to_xml()
 
 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)
-V_mod = np.pi*pt_inner.r**2 - (V_fuel + V_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)
@@ -183,11 +193,10 @@ 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"Flow area, cm² = {V_mod}")
 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}")
-print(f"{V_fuel, V_clad}")
+print(f"DTU ratio = {DTU_ratio}")
 
 
 #> Run the model!