core_1568_pebbles.py 18.7 KB
Newer Older
 Patrick Shriwise committed May 31, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 ``````# -------------- Legend --------------------- # b_ = box # c_ = cell # cz_ = cylinder z # l_ = lattice # m_ = material # px_ = plane x # py_ = plane y # pz_ = plane z # r_ = region # s_ = sphere # u_ = universe # -------------- Modules -------------------- import os import math import numpy as np import openmc import openmc.model import scipy.stats import pandas as pd from math import pi from IPython.display import Image import matplotlib.pyplot as plt # -------------- Functions ------------------ def density_flibe(p, T): # assumes the default value for drho_dp of 1.7324e-7 rho = -0.4884*T + 1.7324e-7*(p - 101325.0) + 2413.0 # kg/m3 return ('kg/m3', rho) # -------------- Arbitrary Parameters ------- debug = 0 # Set 1 to print debug info voxel = 0 # Set 1 to generate 3d voxel plots and 0 for 2d ppm plots random_distribution = 0 # Set 1 for TRISO particles random distribution and 0 for regular lattice # -------------- Units Parameters ----------- mm = 0.1 # default cm μm = 1e-4 # default cm # -------------- Geometry Parameters -------- # TRISO particle radius_fuel = 400.0*μm/2 # UCBTH-14-002, Table 2-1 thickness_c_buffer = 100.0*μm # UCBTH-14-002, Table 2-1 thickness_pyc_inner = 35.0*μm # UCBTH-14-002, Table 2-1 thickness_sic = 35.0*μm # UCBTH-14-002, Table 2-1 thickness_pyc_outer = 35.0*μm # UCBTH-14-002, Table 2-1 radius_c_buffer = radius_fuel + thickness_c_buffer radius_pyc_inner = radius_c_buffer + thickness_pyc_inner radius_sic = radius_pyc_inner + thickness_sic radius_pyc_outer = radius_sic + thickness_pyc_outer packing_fraction = 0.40 # UCBTH-14-002, Table 2-1 if random_distribution==0: pitch_triso_lattice = (radius_fuel+thickness_c_buffer+thickness_pyc_inner+thickness_sic+thickness_pyc_outer)*(4*pi/3/packing_fraction)**0.33333333333333333333333 # Pebble - inner zone graphite; central active zone triso; outer zone graphite # hyuan@anl.gov email April 2, 2019 8:40 PM # In the actual TAMU experiment, the pebble diameter is 22.225mm. # But in CAD, we made it 22.225mm*0.95 = 21.11375 mm To avoid any pebble contacts. # So the pebble radius is 21.11375 mm/2 = 10.556875 mm = 1.0556875 cm `````` Patrick Shriwise committed May 31, 2020 55 ``````# To scale up tp 3 cm pebble radius , I should scale Nek mesh and CAD file by a factor `````` Patrick Shriwise committed May 31, 2020 56 57 58 59 60 ``````# 3/1.0556875 = 2.8417500444 tamu_exp_factor = 1/1.0556875 # Factor to scale geometry to TAMU experiment with pebble radius 1.0556875 cm radius_pebble_inner = 2.5/2 # UCBTH-14-002, Table 2-1 (differs slightly from Cisneros, Table 5-2) ; scaled by TAMU experiment factor radius_pebble_outer = 3.0/2 # UCBTH-14-002, Table 2-1; Cisneros, Table 5-2 ; scaled by TAMU experiment factor radius_pebble_central = radius_pebble_outer - 0.1 # UCBTH-14-002, Table 2-1; Cisneros, Table 5-2 ; scaled by TAMU experiment factor `````` Patrick Shriwise committed May 31, 2020 61 ``````tolerance = 0.00 # Tolerance for pebbles universe filling `````` Patrick Shriwise committed May 31, 2020 62 63 64 ``````radius_pebble_flibe = radius_pebble_outer+tolerance # Pebble centers coordinates (x,y,z) print("Reading pseudo-random pebble centers from file pebble_centers.txt") `````` Patrick Shriwise committed May 31, 2020 65 ``````pebble_centers = np.loadtxt('list_pebbles_1568.csv', delimiter=',', skiprows=1)[:, 1:] `````` Patrick Shriwise committed May 31, 2020 66 67 68 69 70 ``````print("File pebble_centers.txt reading completed") rescaled_file = 'pebble_centers_rescaled.txt' np.savetxt(rescaled_file, pebble_centers.reshape((-1,3))) print("Saved rescaled centers to " + rescaled_file) # Vessel `````` Patrick Shriwise committed May 31, 2020 71 72 ``````extra_thickness = 12.125 # change this value to set keff=1 x_vessel = 0.0 # x position vessel `````` Patrick Shriwise committed May 31, 2020 73 ``````y_vessel = 0.0 # y position vessel `````` Patrick Shriwise committed May 31, 2020 74 75 76 77 78 79 ``````z1_vessel = np.min(pebble_centers[:, 2]) - radius_pebble_outer z2_vessel = np.max(pebble_centers[:, 2]) + radius_pebble_outer radius_vessel = np.max(np.linalg.norm(pebble_centers[:, :-1], axis=1)) + radius_pebble_outer z1_vessel -= extra_thickness # bottom plane vessel ; scaled by TAMU experiment factor z2_vessel += extra_thickness # top plane vessel ; scaled by TAMU experiment factor radius_vessel += extra_thickness # scaled by TAMU experiment factor `````` Patrick Shriwise committed May 31, 2020 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 ``````height_vessel = z2_vessel - z1_vessel # -------------- Printing Parameters --------- print ("GEOMETRY PARAMETERS") print ("TRISO particles radius_fuel [cm] = ", radius_fuel) print ("TRISO particles radius_c_buffer [cm] = ", radius_c_buffer) print ("TRISO particles radius_pyc_inner [cm] = ", radius_pyc_inner) print ("TRISO particles radius_sic [cm] = ", radius_sic) print ("TRISO particles radius_pyc_outer [cm] = ", radius_pyc_outer) print ("TRISO particles packing fraction = ", packing_fraction) if random_distribution==0: print ("TRISO particles regular lattice pitch_triso_lattice [cm] = ", pitch_triso_lattice) print ("Pebble radius_pebble_inner [cm] = ", radius_pebble_inner) print ("Pebble radius_pebble_outer [cm] = ", radius_pebble_outer) print ("Pebble radius_pebble_central [cm] = ", radius_pebble_central) print ("Pebble radius_pebble_flibe [cm] = ", radius_pebble_flibe) print ("Vessel radius_vessel [cm] = ", radius_vessel) print ("Vessel z1_vessel (min z) [cm] = ", z1_vessel) print ("Vessel height_vessel [cm] = ", height_vessel) print ("Pebbles centers [cm] :") print (pebble_centers) # -------------- Materials Parameters -------- # TRISO particle enrichment_uranium = 0.199 fraction_u234 = 2e-3 # U234/(U234+U235) mass from Giacint facility fuel density_fuel = ('g/cm3', 10.5) # UCBTH-14-002, Table 2-1 density_c_buffer = ('g/cm3', 1.0) # Cisneros , Table 3-1 density_pyc = ('g/cm3', 1.87) # Cisneros , Table 3-1 density_sic = ('g/cm3', 3.2) # Cisneros , Table 3-1 # Pebble density_graphite_inner = ('g/cm3', 1.59368) # Cisneros , Table 3-1 density_graphite_outer = ('g/cm3', 1.74) # Cisneros , Table 3-1 # FLiBe coolant enrichment_li7 = 0.99995 temperature_inlet = 273.15 + 600.0 # UCBTH-14-002, Table 1-1 temperature_outlet = 273.15 + 700.0 # UCBTH-14-002, Table 1-1 temperature_flibe = (temperature_inlet+temperature_outlet)/2 # -------------------------------------------------- # NO PARAMETERS DEFINITION BEYOND THIS LINE # -------------- Materials Definition -------------- # TRISO particle # Fuel from Nagley et al. Fabrication of Uranium Oxycarbide Kernels for HTR Fuel https://inldigitallibrary.inl.gov/sites/sti/sti/4886646.pdf Table 2 m_fuel = openmc.Material(name='m_fuel - uranium oxycarbide - triso partciles') m_fuel.add_nuclide('U234', 89.58* enrichment_uranium* fraction_u234 , percent_type='wo') m_fuel.add_nuclide('U235', 89.58* enrichment_uranium*(1-fraction_u234), percent_type='wo') m_fuel.add_nuclide('U238', 89.58*(1-enrichment_uranium) , percent_type='wo') m_fuel.add_nuclide('C0' , 1.80 , percent_type='wo') m_fuel.add_element('O' , 8.62 , percent_type='wo') m_fuel.set_density(*density_fuel) # m_graphite_c_buffer = openmc.Material(name='m_graphite_c_buffer - triso partciles') m_graphite_c_buffer.set_density(*density_c_buffer) m_graphite_c_buffer.add_nuclide('C0', 1.0) m_graphite_c_buffer.add_s_alpha_beta('c_Graphite') # m_graphite_pyc = openmc.Material(name='m_graphite_pyc - triso particles') m_graphite_pyc.set_density(*density_pyc) m_graphite_pyc.add_nuclide('C0', 1.0) m_graphite_pyc.add_s_alpha_beta('c_Graphite') # m_sic = openmc.Material(name='sic - triso partciles') m_sic.add_nuclide('C0' , 1.0) m_sic.add_element('Si', 1.0) m_sic.set_density(*density_sic) # TODO: Add SiC S(alpha, beta) ; ENDF/B-VIIIb4 has data for both Si and C in SiC # m_graphite_matrix = openmc.Material(name='m_graphite_matrix - triso particles') m_graphite_matrix.set_density('g/cm3', 1.6) m_graphite_matrix.add_nuclide('C0', 1.0) m_graphite_matrix.add_s_alpha_beta('c_Graphite') # TODO: Need real specifications for the carbon filler # Pebble m_graphite_inner = openmc.Material(name='m_graphite_inner - pebble inner zone') m_graphite_inner.set_density(*density_graphite_inner) m_graphite_inner.add_nuclide('C0', 1.0) m_graphite_inner.add_s_alpha_beta('c_Graphite') # m_graphite_outer = openmc.Material(name='m_graphite_outer - pebble outer zone') m_graphite_outer.set_density(*density_graphite_outer) m_graphite_outer.add_nuclide('C0', 1.0) m_graphite_outer.add_s_alpha_beta('c_Graphite') # FLiBe coolant - From Cisneros, appendix B, material 24 m_flibe = openmc.Material(name='m_flibe - 2LiF-BeF2') m_flibe.set_density(*density_flibe(101.325e3, temperature_flibe)) m_flibe.add_nuclide('Li7', 2.0* enrichment_li7) m_flibe.add_nuclide('Li6', 2.0*(1 - enrichment_li7)) m_flibe.add_element('Be' , 1.0) m_flibe.add_element('F' , 4.0) # TODO: FLiBe coolant - no S(alpha, beta) data available up to ENDF/B-VIIIb4 # -------------- Geometry Definition -------------- # TRISO particle s_fuel = openmc.Sphere(R=radius_fuel) s_c_buffer = openmc.Sphere(R=radius_c_buffer) s_pyc_inner = openmc.Sphere(R=radius_pyc_inner) s_sic = openmc.Sphere(R=radius_sic) s_pyc_outer = openmc.Sphere(R=radius_pyc_outer) c_triso_fuel = openmc.Cell(name='c_triso_fuel' , fill=m_fuel, region=-s_fuel) c_triso_c_buffer = openmc.Cell(name='c_triso_c_buffer' , fill=m_graphite_c_buffer, region=+s_fuel & -s_c_buffer) c_triso_pyc_inner = openmc.Cell(name='c_triso_pyc_inner', fill=m_graphite_pyc, region=+s_c_buffer & -s_pyc_inner) c_triso_sic = openmc.Cell(name='c_triso_sic' , fill=m_sic, region=+s_pyc_inner & -s_sic) c_triso_pyc_outer = openmc.Cell(name='c_triso_pyc_outer', fill=m_graphite_pyc, region=+s_sic & -s_pyc_outer) c_triso_matrix = openmc.Cell(name='c_triso_matrix' , fill=m_graphite_matrix, region=+s_pyc_outer) u_triso = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer, c_triso_matrix]) # Pebble s_pebble_inner = openmc.Sphere(R=radius_pebble_inner) s_pebble_central = openmc.Sphere(R=radius_pebble_central) s_pebble_outer = openmc.Sphere(R=radius_pebble_outer) c_pebble_inner = openmc.Cell(name='c_pebble_inner' , fill=m_graphite_inner, region=-s_pebble_inner) c_pebble_central = openmc.Cell(name='c_pebble_central', region=+s_pebble_inner & -s_pebble_central) c_pebble_outer = openmc.Cell(name='c_pebble_inner' , fill=m_graphite_outer, region=+s_pebble_central & -s_pebble_outer) c_pebble_flibe = openmc.Cell(name='c_flibe' , fill=m_flibe , region=+s_pebble_outer) # Fill c_pebble_central with TRISO particles if random_distribution==0: # TRISO particles regular lattice using 'pitch_triso_lattice' l_triso = openmc.RectLattice(name='l_triso') l_triso.lower_left = (-pitch_triso_lattice/2, -pitch_triso_lattice/2, -pitch_triso_lattice/2) l_triso.pitch = ( pitch_triso_lattice , pitch_triso_lattice , pitch_triso_lattice) l_triso.outer = u_triso l_triso.universes = np.tile(u_triso, (1,1,1)) # l_triso.universes[0, 0, 0] = universe_triso else: # TRISO particles random distribution using 'packing_fraction' c_triso_pyc_outer_random = openmc.Cell(name='c_triso_pyc_outer_random', fill=m_graphite_pyc, region=+s_sic) u_triso_random = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer_random]) if not os.path.exists('triso_centers.npy'): print("Writing random TRISO centers to file triso_centers.npy") r_triso_random = +s_pebble_inner & -s_pebble_central spheres_random = openmc.model.pack_spheres(radius=radius_pyc_outer, region=r_triso_random, pf=packing_fraction, initial_pf=0.15) triso_random = [openmc.model.TRISO(radius_pyc_outer, u_triso_random, i) for i in spheres_random] centers = np.vstack([i.center for i in triso_random]) if debug==1: print(triso_random[0]) print(triso_random[1]) print(centers.min(axis=0)) print(centers.max(axis=0)) print("Calculated packing fraction of the TRISO particles random distribution = ",len(triso_random)*radius_pyc_outer**3/(radius_pebble_central**3-radius_pebble_inner**3)) np.save('triso_centers.npy', triso_random) print("Reading random TRISO centers from file triso_centers.npy") triso_random = np.load('triso_centers.npy') lower_left, upper_right = c_pebble_central.region.bounding_box shape = (5, 5, 5) pitch = (upper_right - lower_left)/shape l_triso = openmc.model.create_triso_lattice(triso_random, lower_left, pitch, shape, m_graphite_matrix) print("File triso_centers.npy reading completed") c_pebble_central.fill = l_triso u_pebble = openmc.Universe(cells=[c_pebble_inner, c_pebble_central, c_pebble_outer, c_pebble_flibe]) # Vessel cz_vessel = openmc.ZCylinder(x0=x_vessel, y0=y_vessel, R=radius_vessel, boundary_type='reflective') # transmission | reflective | vacuum pz_vessel_bottom = openmc.ZPlane(z0=z1_vessel, boundary_type='reflective') pz_vessel_top = openmc.ZPlane(z0=z2_vessel, boundary_type='reflective') cell_name = ['cell_pebble_' + str(i) for i in range (len(pebble_centers))] `````` Patrick Shriwise committed May 31, 2020 230 ``````s_pebble = [openmc.Sphere(x0=pebble_centers[i,0], y0=pebble_centers[i,1], z0=pebble_centers[i,2], R=radius_pebble_flibe, boundary_type='transmission') for i in range (len(pebble_centers))] `````` Patrick Shriwise committed May 31, 2020 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 ``````c_pebble = [openmc.Cell(name=cell_name[i], fill=u_pebble, region=-s_pebble[i]) for i in range (len(pebble_centers))] r_vessel = -cz_vessel & +pz_vessel_bottom & -pz_vessel_top for i in range(len(pebble_centers)): c_pebble[i].translation = pebble_centers[i] r_vessel = r_vessel & +s_pebble[i] c_vessel = openmc.Cell(name='c_vessel', fill=m_flibe, region = r_vessel) # # Pebbles pseudo-random distribution using openmc.model.TRISO function - IT DOES NOT WORK WITH TALLY CARDS # c_pebble_outer_random = openmc.Cell(name='c_pebble_outer_random', fill=m_graphite_outer, region=+s_pebble_outer) # u_pebble_random = openmc.Universe(cells=[c_pebble_inner, c_pebble_central, c_pebble_outer]) # print("Reading random pebble centers from file pebble_centers.txt") # spheres_random = np.loadtxt('pebble_centers.txt') # spheres_random = spheres_random*radius_pebble_outer # pebbles centers have cm unit and must be scaled by TAMU experiment factor # print("File pebble_centers.txt reading completed") # print(spheres_random) # pebble_random = [openmc.model.TRISO(radius_pebble_outer, u_pebble_random, i) for i in spheres_random] # lower_left, upper_right = c_vessel.region.bounding_box # shape = (5, 5, 5) `````` Patrick Shriwise committed May 31, 2020 249 ``````# pitch = (upper_right - lower_left)/shape `````` Patrick Shriwise committed May 31, 2020 250 251 252 253 254 ``````# l_pebble = openmc.model.create_triso_lattice(pebble_random, lower_left, pitch, shape, m_flibe) # c_vessel.fill = l_pebble # # Global geometry # c_all = [[c_vessel, c_pebble[i]] for i in range (len(pebble_centers))] # The for loop does not work `````` Patrick Shriwise committed May 31, 2020 255 ``````u_zero = openmc.Geometry([c_vessel] + list(c_pebble)) `````` Patrick Shriwise committed May 31, 2020 256 257 258 259 260 261 262 263 264 ``````# -------------- Printing Cells -------------- print ("Cell of the vessel:") print (c_vessel) print ("Cells of the pebbles:") print (c_pebble) # -------------- Settings -------------------- settings = openmc.Settings() settings.source = openmc.Source(space=openmc.stats.Box(*c_vessel.bounding_box, only_fissionable=True)) settings.particles = 10000 `````` Patrick Shriwise committed May 31, 2020 265 ``````settings.inactive = 50 `````` Patrick Shriwise committed May 31, 2020 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 ``````settings.batches = 150 settings.temperature = dict(default=573, method='interpolation', multipole=True, range=(300.0, 1500.0), tolerance=1000.0) # Fuel volume calculation volume_fuel = openmc.VolumeCalculation([m_fuel], 10000000, *c_pebble_central.bounding_box) settings.volume_calculations = [volume_fuel] model = openmc.model.Model(geometry=u_zero, settings=settings) model.export_to_xml() # -------------- Tally ----------------------- # print(l_pebble) # print(c_pebble_central) # print(u_pebble_random) # print(pebble_random) # # filter1 = openmc.DistribcellFilter([13]) # filter1 = openmc.UniverseFilter([u_pebble_random]) # filter2 = openmc.CellFilter([13]) # tally = openmc.Tally(name='tally kappa-fission') # tally.filters = [filter1, filter2] # tally.scores = ['kappa-fission'] # # tally = openmc.Tally(name='tally kappa-fission') # tally.filters = [openmc.CellFilter([c_pebble_1, c_pebble_2])] # tally.scores = ['kappa-fission'] # tallies = openmc.Tallies([tally]) # tallies.export_to_xml() # -------------- Plots -------------- # Plot parameters plot_width = max(2*radius_vessel,height_vessel) plot_zcenter = z1_vessel+height_vessel/2 m_colors = {m_fuel: 'brown', m_graphite_c_buffer: 'LightSteelBlue', m_graphite_pyc: 'blue', m_sic: 'orange', m_graphite_matrix: 'cyan', m_graphite_inner: 'DeepSkyBlue', m_graphite_outer: 'Navy', m_flibe: 'yellow'} # plot1 = openmc.Plot() plot1.filename = 'plot1' plot1.width = (2*radius_vessel, 2*radius_vessel) plot1.basis = 'xy' plot1.origin = (0, 0, 0) plot1.pixels = (1000, 1000) plot1.color_by = 'material' plot1.colors = m_colors # plot2 = openmc.Plot() plot2.filename = 'plot2' plot2.width = (2*radius_vessel, 2*radius_vessel) plot2.basis = 'xy' plot2.origin = (0, 0, 1) plot2.pixels = (1000, 1000) plot2.color_by = 'material' plot2.colors = m_colors # plot3 = openmc.Plot() plot3.filename = 'plot3' plot3.width = (2*radius_vessel, 2*radius_vessel) plot3.basis = 'xy' plot3.origin = (0, 0, 1) plot3.pixels = (1000, 1000) plot3.color_by = 'material' plot3.colors = m_colors # plot4 = openmc.Plot() plot4.filename = 'plot4' plot4.width = (plot_width, plot_width) plot4.basis = 'xz' plot4.origin = (0, 0 ,plot_zcenter) plot4.pixels = (1000, 1000) plot4.color_by = 'material' plot4.colors = m_colors # plots = openmc.Plots([plot1, plot2, plot3, plot4]) if voxel==1: plotv1 = openmc.Plot() plotv1.filename = 'plotv1' plotv1.width = (2*radius_pebble_outer, 2*radius_pebble_outer, 2*radius_pebble_outer) plotv1.origin = (0, 0 ,0) plotv1.pixels = (300, 300, 300) plotv1.color_by = 'material' plotv1.colors = m_colors plotv1.type = 'voxel' plots = openmc.Plots([plotv1]) plots.export_to_xml()``````