core_1568_pebbles.py 18.7 KB
Newer Older
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
55
# To scale up tp 3 cm pebble radius ,  I should scale Nek mesh and CAD file by a factor
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
61
tolerance             = 0.00                      # Tolerance for pebbles universe filling
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")
65
pebble_centers        = np.loadtxt('list_pebbles_1568.csv', delimiter=',', skiprows=1)[:, 1:]
66
print("File pebble_centers.txt reading completed")
67
rescaled_file = 'pebble_centers.txt'
68 69 70
np.savetxt(rescaled_file, pebble_centers.reshape((-1,3)))
print("Saved rescaled centers to " + rescaled_file)
# Vessel
71 72
extra_thickness       =  12.125  # change this value to set keff=1
x_vessel              =  0.0  # x position vessel
73
y_vessel              =  0.0  # y position vessel
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
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))]
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))]
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)
249
#  pitch                   = (upper_right - lower_left)/shape
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
255
u_zero                  = openmc.Geometry([c_vessel] + list(c_pebble))
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
265
settings.inactive  = 50
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()