Added simple VWH scheme for comparison

1 files changed, 66 insertions, 7 deletions

diff --git a/hpmr.py b/hpmr.py
index c9e8f18..5d07188 100755
--- a/hpmr.py
+++ b/hpmr.py
@@ -1,5 +1,6 @@
 #!/bin/python3
 import openmc
+import math
 
 ############ MATERIALS ############ 
 # All materials are specified in the report on HEAT PIPE MICROREACTOR MODELING WITH BLUECRAB
@@ -55,6 +56,7 @@ mat_mod_helium.add_element("He",1.0)
 mat_mod_helium.set_density("g/cc",0.18e-3)
 mat_list.append(mat_mod_helium)
 
+# Composition from PNNL-15870r2
 mat_mod_ss316 = openmc.Material(name="Moderator Pin SS316")
 mat_mod_ss316.add_element("C",0.000800,"wo")
 mat_mod_ss316.add_element("Mn",0.020000,"wo")
@@ -82,6 +84,7 @@ mat_pipe_helium.add_element("He",1.0)
 mat_pipe_helium.set_density("g/cc",0.18e-3)
 mat_list.append(mat_pipe_helium)
 
+# Composition from PNNL-15870r2
 mat_pipe_ss316 = openmc.Material(name="Heat Pipe SS316")
 mat_pipe_ss316.add_element("C",0.000800,"wo")
 mat_pipe_ss316.add_element("Mn",0.020000,"wo")
@@ -113,9 +116,6 @@ mat_pipe_k_gas.add_element("K",1.0)
 mat_pipe_k_gas.set_density("g/cc",1.11e-4)
 mat_list.append(mat_pipe_k_gas)
 
-## Export ##
-materials = openmc.Materials(mat_list)
-materials.export_to_xml()
 
 
 ############ GEOMETRY ############ 
@@ -169,6 +169,62 @@ geo_triso_compact_surrounding_cell = openmc.Cell(name="TRISO Compact Surrounding
 
 geo_triso_compact_universe = openmc.Universe(name="TRISO Compact Universe",cells=[geo_triso_compact_fuel_cell,geo_triso_compact_surrounding_cell])
 
+## VWH TRISO Compact ##
+geo_vwh_fracs = []
+geo_n_triso = len(geo_triso_compact_packing)
+print(f'{geo_n_triso = }')
+
+# Get volume fractions for mixing
+geo_triso_compact_v_tot = math.pi*geo_constant_triso_radius**2*2*geo_constant_height
+
+geo_triso_compact_v_uco = geo_n_triso*4/3*math.pi*geo_triso_spheres[0].r**3
+geo_vwh_fracs.append(geo_triso_compact_v_uco/geo_triso_compact_v_tot)
+
+geo_triso_compact_v_buffer = geo_n_triso*4/3*math.pi*geo_triso_spheres[1].r**3 - geo_triso_compact_v_uco
+geo_vwh_fracs.append(geo_triso_compact_v_buffer/geo_triso_compact_v_tot)
+
+geo_triso_compact_v_pyc1 = geo_n_triso*4/3*math.pi*geo_triso_spheres[2].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer)
+geo_vwh_fracs.append(geo_triso_compact_v_pyc1/geo_triso_compact_v_tot)
+
+geo_triso_compact_v_sic = geo_n_triso*4/3*math.pi*geo_triso_spheres[3].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1)
+geo_vwh_fracs.append(geo_triso_compact_v_sic/geo_triso_compact_v_tot)
+
+geo_triso_compact_v_pyc2 = geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3 - (geo_triso_compact_v_uco + geo_triso_compact_v_buffer + geo_triso_compact_v_pyc1 + geo_triso_compact_v_sic)
+geo_vwh_fracs.append(geo_triso_compact_v_pyc2/geo_triso_compact_v_tot)
+
+geo_triso_compact_v_graphite = geo_triso_compact_v_tot - geo_n_triso*4/3*math.pi*geo_triso_spheres[4].r**3
+geo_vwh_fracs.append(geo_triso_compact_v_graphite/geo_triso_compact_v_tot)
+
+# Strip TSL temporarily for mixing
+mat_tsl_saved = {}
+for mat in mat_list:
+    if hasattr(mat, "_sab") and mat._sab:
+        mat_tsl_saved[mat]=list(mat._sab)
+        mat._sab=[]
+
+# Now do a simple volume homogenization
+mat_triso_vwh = openmc.Material.mix_materials(
+        materials=[mat_triso_uco, mat_triso_buffer, mat_triso_pyc1, mat_triso_sic, mat_triso_pyc2, mat_graphite_matrix],
+        fracs=geo_vwh_fracs,
+        percent_type="vo",
+        name="VWH Homogenized TRISO Compact"
+        )
+
+# Restore TSL afterwards
+for mat, tsl_list in mat_tsl_saved.items():
+    mat._sab = tsl_list
+
+mat_triso_vwh.add_s_alpha_beta("c_Graphite")
+
+# Add to list
+mat_list.append(mat_triso_vwh)
+
+# Cast into cells and a universe
+geo_triso_vwh_fuel_cell = openmc.Cell(name="TRISO VWH Cell",fill=mat_triso_vwh,region=geo_triso_compact_region)
+geo_triso_vwh_surrounding_cell = openmc.Cell(name="TRISO VWH Surrounding Cell",fill=mat_graphite_matrix,region=~geo_triso_compact_region)
+
+geo_triso_vwh_universe = openmc.Universe(name="TRISO VWH Universe",cells=[geo_triso_vwh_fuel_cell,geo_triso_vwh_surrounding_cell])
+
 ## Heat Pipe ##
 geo_pipe_cylinders = [openmc.ZCylinder(r=rad, x0=0.0, y0=0.0) for rad in [0.80, 0.90, 0.97, 1.05, 1.07]]
 geo_pipe_cells = [openmc.Cell(name="Heat Pipe K Gas",fill=mat_pipe_k_gas,region=-geo_pipe_cylinders[0]),
@@ -201,11 +257,11 @@ geo_assembly_lattice_map = []
 for i in range(0,7):
     if i%2 == 0:
         ring = [geo_pipe_universe]*(36-6*i)
-        ring[1::2] = [geo_triso_compact_universe]*len(ring[1::2])
+        ring[1::2] = [geo_triso_vwh_universe]*len(ring[1::2])
         geo_assembly_lattice_map.append(ring)
     else:
         ring = [geo_mod_universe]*(36-6*i)
-        ring[::2] = [geo_triso_compact_universe]*len(ring[::2])
+        ring[::2] = [geo_triso_vwh_universe]*len(ring[::2])
         geo_assembly_lattice_map.append(ring)
 
 geo_assembly_lattice_map[-1] = [geo_pipe_universe]
@@ -222,7 +278,6 @@ geo_assembly_region = openmc.model.HexagonalPrism(
 
 # Clip the lattice to the space and make it root
 geo_assembly_cell = openmc.Cell(name="Unit Assembly Cell",fill=geo_assembly_lattice,region=-geo_assembly_region&+geo_constant_z_min&-geo_constant_z_max)
-#geo_outer_cell = openmc.Cell(name="Unit Assembly Void Outer",fill=None,region=~(geo_assembly_cell.region))
 geo_assembly_universe = openmc.Universe(cells=[geo_assembly_cell])
 
 geometry = openmc.Geometry(geo_assembly_universe)
@@ -249,7 +304,11 @@ set_entropy_mesh.upper_right=(+geo_constant_lattice_pitch,+geo_constant_lattice_
 set_entropy_mesh.dimension=(30,30,30)
 
 settings.entropy_mesh=set_entropy_mesh
+settings.export_to_xml()
+
+## Export ##
+materials = openmc.Materials(mat_list)
+materials.export_to_xml()
 
 # Run!
-settings.export_to_xml()
 openmc.run()