openmc_model.py 19.4 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 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 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 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 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
import os
import math
import numpy as np
import openmc
import openmc.model
import pandas as pd
import scipy.stats
from   math   import pi
from   IPython.display import Image
import matplotlib.pyplot as plt

# -------------- Unit Conversions -----------
mm                    = 0.1   # default cm
μm                    = 1e-4  # default cm

# -------------- Output Parameters -------
voxel = False                  # Set 1 to generate 3d voxel plots and 0 for 2d ppm plots
random_distribution = False    # Set 1 for TRISO particles random distribution and 0 for regular (simple cubic) lattice
verbose = True

# -------------- Functions ------------------
def density_flibe(p, t):
    """
    Returns FLiBe density at a pressure, p,
    and temperature, t in units of kg/m3
    """
    # 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 rho

# -------------- 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_pyc_outer)*(4*pi/3/packing_fraction)**(1/3)

# Pebble
radius_pebble_inner   = 2.5/2                     # UCBTH-14-002, Table 2-1
radius_pebble_outer   = 3.0/2                     # UCBTH-14-002, Table 2-1
radius_pebble_central = radius_pebble_outer - 0.1 # UCBTH-14-002, Table 2-1
tolerance             = 0.00                      # Tolerance for pebbles universe filling
radius_pebble_flibe   = radius_pebble_outer+tolerance

# Pebble centers coordinates (x,y,z)
pebble_centers_file = 'pebble_centers_rescaled.txt'
print("************************************************************")
print("This script expects pebble centers scaled to radius = 1.5 cm")
print("************************************************************")
print("Reading pebble centers from file {}...".format(pebble_centers_file), end="")
pebble_centers = np.loadtxt(pebble_centers_file, delimiter=' ')
print("Read {} pebble positions\n".format(len(pebble_centers)))

# Scale the pebble dimensions (radius of 1.0) to actual dimensions
for i in range(len(pebble_centers)):
  pebble_centers[i] *= radius_pebble_outer

print("  Min pebble (center) z coordinate = {}".format(np.min(pebble_centers[:,2])))
print("  Max pebble (center) z coordinate = {}".format(np.max(pebble_centers[:,2])))
print("  Max pebble (center) r coordinate = {}".format(np.max(np.linalg.norm(pebble_centers[:, :-1], axis=1))))
print("  Min pebble (center) r coordinate = {}\n".format(np.min(np.linalg.norm(pebble_centers[:, :-1], axis=1))))

# Vessel Dimensions - we add a small tolerance to make sure no pebbles are cut off. The vesel
# is the region that actually holds the pebbles
tol = 0.00001
vessel_x, vessel_y = (0.0, 0.0)
vessel_z_min = np.min(pebble_centers[:, 2]) - radius_pebble_outer - tol
vessel_z_max = np.max(pebble_centers[:, 2]) + radius_pebble_outer + tol
vessel_outer_radius = np.max(np.linalg.norm(pebble_centers[:, :-1], axis=1)) + radius_pebble_outer + tol
vessel_inner_radius = np.min(np.linalg.norm(pebble_centers[:, :-1], axis=1)) - radius_pebble_outer - tol
vessel_height = vessel_z_max - vessel_z_min

# Add outer reflector of thickness 40 cm, UCBTH-14-002. We assume that there is then also a reflector
# of the same thickness on the top and bottom (an arbitrary selection). This reflector surrounds the pebble bed
# around its outer radius as well as forms the center reflector column.
reflector_thickness = 40.0
reflector_z_min = vessel_z_min - reflector_thickness
reflector_z_max = vessel_z_max + reflector_thickness
reflector_outer_radius = vessel_outer_radius + reflector_thickness

# -------------- Printing Parameters ---------
if verbose:
    print ("GEOMETRY PARAMETERS")
    print ("  TRISO particle fuel radius       [cm] = {}".format(radius_fuel))
    print ("  TRISO particle buffer radius     [cm] = {}".format(radius_c_buffer))
    print ("  TRISO particle iPyC radius       [cm] = {}".format(radius_pyc_inner))
    print ("  TRISO particle SiC radius        [cm] = {}".format(radius_sic))
    print ("  TRISO particle oPyC radius       [cm] = {}".format(radius_pyc_outer))
    print ("  TRISO particles packing fraction [--] = {}".format(packing_fraction))
    print ("  TRISO lattice pitch              [cm] = {}".format(pitch_triso_lattice))
    if random_distribution==0:
      print ("  TRISO particles regular lattice  [cm] = {}".format(pitch_triso_lattice))
    print ("  Pebble core radius               [cm] = {}".format(radius_pebble_inner))
    print ("  Pebble shell outer radius        [cm] = {}".format(radius_pebble_outer))
    print ("  Pebble shell inner radius        [cm] = {}".format(radius_pebble_central))
    print ("  Vessel outer radius              [cm] = {}".format(vessel_outer_radius))
    print ("  Vessel inner radius              [cm] = {}".format(vessel_inner_radius))
    print ("  Vessel min z                     [cm] = {}".format(vessel_z_min))
    print ("  Vessel height                    [cm] = {}".format(vessel_height))

# -------------- 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 Graphite
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
rho_flibe            = density_flibe(101.325e3, temperature_flibe)
density_flibe        = ('kg/m3', rho_flibe)

# Reflector, a mixture of graphite and flibe.
reflector_porosity   = 0.112                # Novak thesis, Table 6.4, assuming 5 mm gaps b/w blocks
rho_graphite         = 1632.0               # Novak thesis, page 234
density_reflector    = ('kg/m3', reflector_porosity * rho_flibe + (1 - reflector_porosity) * rho_graphite)

# --------------------------------------------------
# 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 particles')
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)

# Buffer
m_graphite_c_buffer = openmc.Material(name='m_graphite_c_buffer - triso particles')
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')

# Pyrolitic carbon
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')

# Silicon carbide
m_sic = openmc.Material(name='sic - triso particles')
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

# Graphite matrix
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 inner graphite core
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')

# Pebble outer graphite shell
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)
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

# Graphite-flibe homogenized reflector
flibe_molar_mass = m_flibe.average_molar_mass * 7
graphite_molar_mass = m_graphite_outer.average_molar_mass

# moles of flibe in 1 m^3:
#            0.112 m^3          * g/m^3                * mol / g
mols_flibe = reflector_porosity * (rho_flibe * 1000.0) / flibe_molar_mass

# moles of graphite in 1 m^3:
#            0.888 m^3                  * g/m^3                   * mol / g
mols_graphite = (1.0 - reflector_porosity) * (rho_graphite * 1000.0) / graphite_molar_mass

# Then, the absolute number of moles is irrelevant since we're adding nuclides based on
# the atomic representation
m_reflector = openmc.Material(name='m_reflector')
m_reflector.add_nuclide('Li7', 2.0 * enrichment_li7 * mols_flibe)
m_reflector.add_nuclide('Li6', 2.0 * (1 - enrichment_li7) * mols_flibe)
m_reflector.add_element('Be' , 1.0 * mols_flibe)
m_reflector.add_element('F'  , 4.0 * mols_flibe)
m_reflector.add_nuclide('C0', 1.0 * mols_graphite)
m_reflector.add_s_alpha_beta('c_Graphite')
m_reflector.set_density(*density_reflector)

# -------------- Geometry Definition --------------
# TRISO particle universe
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='TRICO Fuel', fill=m_fuel, region=-s_fuel)
c_triso_c_buffer   = openmc.Cell(name='TRISO Graphite Buffer', fill=m_graphite_c_buffer, region=+s_fuel & -s_c_buffer)
c_triso_pyc_inner  = openmc.Cell(name='TRISO Pyrolitic Graphite Inner', fill=m_graphite_pyc, region=+s_c_buffer  & -s_pyc_inner)
c_triso_sic        = openmc.Cell(name='TIRSO Silicone Carbide', fill=m_sic, region=+s_pyc_inner & -s_sic)
c_triso_pyc_outer  = openmc.Cell(name='TRISO Pyrolitic Graphite Outer', fill=m_graphite_pyc, region=+s_sic & -s_pyc_outer)
c_triso_matrix     = openmc.Cell(name='TRISO Graphite Matrix', fill=m_graphite_matrix, region=+s_pyc_outer)

triso_cells = [c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer, c_triso_matrix]
u_triso = openmc.Universe(cells=triso_cells)

# Pebble Geometry
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='Pebble graphite inner region', fill=m_graphite_inner, region=-s_pebble_inner)
c_pebble_central = openmc.Cell(name='Pebble central region (TRISOs)', region=+s_pebble_inner & -s_pebble_central)
c_pebble_outer = openmc.Cell(name='Pebble graphite outer region', fill=m_graphite_outer, region=+s_pebble_central & -s_pebble_outer)
c_pebble_flibe = openmc.Cell(name='Pebble exterior (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='Pebble TRISO Lattice')
    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, (20, 20, 20))
# TRISO particles random distribution using 'packing_fraction'
else:
    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 verbose:
            packing_fraction = len(triso_random)*radius_pyc_outer**3/(radius_pebble_central**3-radius_pebble_inner**3)
            print("Calculated packing fraction of the TRISO particles random distribution = {}".format(packing_fraction))
            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                      = (20, 20, 20)
        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")

# Fil central pebble cell with the TRISO lattice
c_pebble_central.fill   = l_triso

pebble_cells = [c_pebble_inner, c_pebble_central, c_pebble_outer, c_pebble_flibe]
u_pebble = openmc.Universe(cells=pebble_cells)

# Vessel cell - this is the pebble-containing region
vessel_outer = openmc.ZCylinder(x0=vessel_x, y0=vessel_y, r=vessel_outer_radius)
vessel_inner = openmc.ZCylinder(x0=vessel_x, y0=vessel_y, r=vessel_inner_radius)
vessel_bottom = openmc.ZPlane(z0=vessel_z_min)
vessel_top = openmc.ZPlane(z0=vessel_z_max)
vessel_region = -vessel_outer & +vessel_bottom & -vessel_top & +vessel_inner
vessel_cell = openmc.Cell(name='Pebble Vessel', region = vessel_region)

# Outer reflector cell
reflector_outer = openmc.ZCylinder(x0=vessel_x, y0=vessel_y, r=reflector_outer_radius, boundary_type = 'reflective')
reflector_bottom = openmc.ZPlane(z0=reflector_z_min, boundary_type = 'reflective')
reflector_top = openmc.ZPlane(z0=reflector_z_max, boundary_type = 'reflective')
outer_reflector_region = -reflector_outer & +vessel_bottom & -vessel_top & +vessel_inner
inner_reflector_region = -vessel_inner & +vessel_bottom & -vessel_top
top_reflector_region = +vessel_top & -reflector_top & -reflector_outer
bottom_reflector_region = -vessel_bottom & +reflector_bottom & -reflector_outer
outer_reflector_cell = openmc.Cell(name = 'Outer Reflector', region = outer_reflector_region, fill = m_reflector)
inner_reflector_cell = openmc.Cell(name = 'Inner Reflector', region = inner_reflector_region, fill = m_reflector)
top_reflector_cell = openmc.Cell(name = 'Top Reflector', region = top_reflector_region, fill = m_reflector)
bottom_reflector_cell = openmc.Cell(name = 'Bottom Reflector', region = bottom_reflector_region, fill = m_reflector)

# Creating TRISOs for the pebbles to pack them into a lattice for efficiency
cell_name = ['cell_pebble_' + str(i) for i in range (len(pebble_centers))]
pebble_trisos = []
for center, name in zip(pebble_centers, cell_name):
    pebble_trisos.append(openmc.model.TRISO(radius_pebble_outer, u_pebble, center))
    pebble_trisos[-1].name = name

# Place pebbles into a lattice
llc_vessel, urc_vessel = vessel_region.bounding_box
313
l_pebble_shape = np.asarray((40, 40, 40))
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 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391
l_pebble_pitch = (urc_vessel - llc_vessel) / l_pebble_shape
l_pebble = openmc.model.create_triso_lattice(pebble_trisos, llc_vessel, l_pebble_pitch, l_pebble_shape, m_flibe)

# fill vessel with pebble lattice
vessel_cell.fill = l_pebble

geom = openmc.Geometry([vessel_cell, outer_reflector_cell, inner_reflector_cell, top_reflector_cell, bottom_reflector_cell])

if verbose:
    print ("Cell of the vessel:")
    print (vessel_cell)

# -------------- Settings --------------------
settings           = openmc.Settings()
settings.source    = openmc.Source(space=openmc.stats.Box(*vessel_cell.bounding_box))
settings.particles = 50000
settings.inactive  = 50
settings.batches   = 250
settings.temperature = dict(default=573, method='interpolation', multipole=True, range=(300.0, 1500.0), tolerance=1000.0)
settings.material_cell_offsets = False
settings.output = { 'summary' : False,
                    'tallies' : False }
# 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=geom, settings=settings)
model.export_to_xml()

# -------------- Plots --------------
# Plot parameters
vessel_diameter = 2.0 * reflector_outer_radius
pebble_diameter = 2.0 * radius_pebble_outer
extra = 1.1

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',
            m_reflector:         'slategray'}
#
plot1          = openmc.Plot()
plot1.filename = 'plot1'
plot1.width    = (extra * vessel_diameter, extra * vessel_diameter)
plot1.basis    = 'xy'
plot1.origin   = (0, 0, 0)
plot1.pixels   = (2000, 2000)
plot1.color_by = 'material'
plot1.colors   = m_colors
#
plot2          = openmc.Plot()
plot2.filename = 'plot2'
plot_width   = extra * max(vessel_diameter, vessel_height)
plot2.width    = (plot_width, plot_width)
plot2.basis    = 'xz'
plot_zcenter = vessel_z_min + vessel_height/2.0
plot2.origin   = (0, 0, plot_zcenter)
plot2.pixels   = (2000, 2000)
plot2.color_by = 'material'
plot2.colors   = m_colors
#
plots          = openmc.Plots([plot1, plot2])
if voxel==1:
   plotv1 = openmc.Plot()
   plotv1.filename = 'plotv1'
   plotv1.width    = (pebble_diameter, pebble_diameter, pebble_diameter)
   plotv1.origin   = (pebble_centers[0,0], pebble_centers[0,1] ,pebble_centers[0,2])
   plotv1.pixels   = (300, 300, 300)
   plotv1.color_by = 'material'
   plotv1.colors   = m_colors
   plotv1.type     = 'voxel'
   plots.append(plotv1)

plots.export_to_xml()