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

4 files changed, 48 insertions, 8 deletions

diff --git a/.gitignore b/.gitignore
index 6141a39..5e3ad64 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 *.h5
 *.pkl
 *.csv
+*-results/
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!
diff --git a/clean.sh b/clean.sh
new file mode 100755
index 0000000..ca8b779
--- /dev/null
+++ b/clean.sh
@@ -0,0 +1,5 @@
+#!/bin/bash
+
+# Clean up the workspace
+
+rm -f *.h5 *.xml
diff --git a/driver-37.sh b/driver-37.sh
new file mode 100755
index 0000000..31a81f5
--- /dev/null
+++ b/driver-37.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+# Driver script for parametric pin size analysis of the 
+# 37-element CANDU fuel bundle. AnnCon 2026 paper.
+# Connor Moore, 2026. <connor.moore@ontariotechu.net>
+
+green="\033[32m"
+reset="\033[0m"
+
+mkdir -p 37-results/
+touch 37-results/37_keff_table
+
+for r in $(seq 0.1 0.05 0.8); do
+    echo -e ${green}Running rf = $r cm... ${reset}
+    filename=37-results/$r.log
+    python3 candu-37.py -rf $r > $filename
+    keff=$(cat $filename | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "+" -f 1)
+    pm=$(cat $filename | grep "Combined k-effective" | cut -d "=" -f 2 | cut -d "-" -f 2)
+    flow_area=$(cat $filename | grep "Flow area" | cut -d "=" -f 2)
+    fuel_mass=$(cat $filename | grep "Fuel mass" | cut -d "=" -f 2)
+    clad_mass=$(cat $filename | grep "Cladding mass" | cut -d "=" -f 2)
+    dtu=$(cat $filename | grep "DTU ratio" | cut -d "=" -f 2)
+
+    echo $keff +/- $pm, $flow_area cm² flow, $fuel_mass g/cm fuel, $clad_mass g/cm clad, $dtu DTU ratio | cowsay
+done