Building a Mesh of the PBFHR

The PBFHR design document will define the dimensions for mesh generation of the core of the PBFHR. This mesh and the materials described will be used to inform Serpent, MCNP, and MOOSE analyses of the PBFHR neutronics.


In [1]:
import numpy as np
from pyne import mesh, material 
from pyne.xs.channels import sigma_t
from yt.config import ytcfg; ytcfg["yt","suppressStreamLogging"] = "True"
from yt.mods import *
from reactor_element import *


/Users/khuff/anaconda/lib/python2.7/site-packages/pyne/mesh.py:12: VnVWarning: pyne.mesh is not yet V&V compliant.
  warn(__name__ + " is not yet V&V compliant.", VnVWarning)
/Users/khuff/anaconda/lib/python2.7/site-packages/pyne/mesh.py:18: VnVWarning: the PyTAPS optional dependency could not be imported. Some aspects of the mesh module may be incomplete.
  "Some aspects of the mesh module may be incomplete.", VnVWarning)
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-8b1ce18de60c> in <module>()
      1 import numpy as np
----> 2 from pyne import mesh, material
      3 from pyne.xs.channels import sigma_t
      4 from yt.config import ytcfg; ytcfg["yt","suppressStreamLogging"] = "True"
      5 from yt.mods import *

/Users/khuff/anaconda/lib/python2.7/site-packages/pyne/mesh.py in <module>()
     18          "Some aspects of the mesh module may be incomplete.", VnVWarning)
     19 
---> 20 from pyne.material import Material, MaterialLibrary, MultiMaterial
     21 
     22 if sys.version_info[0] > 2:

ImportError: dlopen(/Users/khuff/anaconda/lib/python2.7/site-packages/pyne/material.so, 2): Library not loaded: @loader_path/../../../libhdf5.7.dylib
  Referenced from: /Users/khuff/anaconda/lib/python2.7/site-packages/pyne/material.so
  Reason: image not found

In [ ]:
def in_sphere(coord, radius, position):
    x0=position[0]
    y0=position[1]
    z0=position[2]
    rsq = (coord[0]-x0)**2 + (coord[1]-y0)**2 + (coord[2]-z0)**2
    if (rsq <= radius**2):
        return True
    else:
        return False
    
def add_sphere(m, radius, position, material):
    for i, mat, ve in m:
        coord = m.mesh.getVtxCoords(ve)
        if in_sphere(coord,radius,position):
            m.mats[i] = material

Core Geometry

The vessel has :

  • an outer graphite reflector,
  • an annular core full of pebbles,
  • salt, and control elements,
  • and a central cylindrical reflector.

For simplicity, the mesh will represent the reactor vessel as a cube. The outer reflector blocks will occupy the space outside the cylindrical core. The core has five regions. From the bottom up, there is a fuelling chute, an expasion region, the core, the compression region, and a defuelling chute. All are full of pebbles and salt.


In [ ]:
def slant_r(pt1, pt2, z):
    r1, z1 = pt1
    r2, z2 = pt2
    r = (float((r2-r1))/float((z2-z1)))*z + r1
    
def inner_radius(z):
    # fueling chute
    if z <= 1.118 : 
        return 0.45
    # expansion region
    elif z <=1.12912:
        return slant_r([1.118,0.45],[1.12912,0.35],z)
    # core
    elif z <= 4.148:
        return 0.35
    # compression
    elif z <= 4.77154:
        return slant_r([4.148,0.35],[4.77154,0.71],z)
    # defueling
    else :
        return 0.71
    
def outer_radius(z):
       # fueling chute
    if z <= 1.118 : 
        return 0.8574
    # expansion region
    elif z <=1.648:
        return slant_r([1.118,0.8574],[1.648,1.25],z)
    # core
    elif z <= 4.148:
        return 1.25
    # compression
    elif z <= 4.77154:
        return slant_r([4.148,1.25],[4.77154,0.89],z)
    # defueling
    else :
        return 0.89

height = 5.7715
total_radius = 1.65

Material 1 : Flibe

"A key issue for the use of flibe is the requirement to use lithium that has been enriched from the natural concentration of 92.4% to approximately 99.995% 7Li. The U.S. ceased production of enriched lithium in the 1960s after producing a large inventory of 6Li and a smaller inventory of 7Li."

Isotopic Composition

The composition of the flibe coolant is Li2 Be F4 with 99.995% enriched 7Li.

  • That means 7Li= 99.995% and 6Li=0.005%
  • F natural abundance = 100% 19F
  • Be natural abundance = 100% 9Be

From Yousry Gohar ( http://aries.ucsd.edu/LIB/MEETINGS/0103-TRANSMUT/gohar/Gohar-present.pdf ):

Physical Properties

  • Melting point, °C 459.1
  • Thermal conductivity, W/cm. °C 0.010
  • Viscosity, centipoises 0.116 exp [3755/T(°k)]
  • Electrical conductivity at 500 , ohm-1cm-1 9.2 10-3
  • Heat capacity, cal/g. °C 0.57
  • Density, g/cm3 2.214 - 4.2 x 10-4 T(°C)

In [ ]:
flibe = material.from_atom_frac({'F19': 2./3., 
                                 'Li6': 0.00005/3., 
                                 'Li7': 0.99995/3., 
                                 'Be9':1./6.}, 
                                mass=1.0, 
                                density=2.214)

Fuel

The fuel is a composite material made of triso pebbles. The composite material composition depends heavily on the triso packing density and the enrichment of the uranium.

The canonical makeup of the Triso kernels is $$UC_{0.5}O_{1.5}$$. The kernels are embedded in a graphite matrix. The fuel and graphite matrix exists in an annular shell around a pure graphite core.

From PB-FHR-DESIGN-CALC-002-00-Reactivity-Effect-Flibe-Impurities.pdf

The FHR pebble unit cell transport model is composed of 3cm OD TRISO pebbles in a BCC lattice with mirror reflection and 40% flibe by volume. The 500μm OD fuel kernels are composed of UO1.5C0.5 enriched to 19.9% 235U. Kernels are spaced to bring about a 300 carbon-to-heavy-metal atomic ratio within the pebble. All material cross-sections are broadened to 900K except the graphite thermal scattering kernels, which are broadened to 1000K. 6Li is depleted to 50wppm on a lithium mass basis within the flibe. The MCNP input file for the clean salt is included in Appendix A.

Pebble Geometry

The pebbles are 3.0 cm in diameter.


In [18]:
d_pebble = 0.03 # [m]
r_pebble = d_pebble/2. # [m]
d_pebble_core = 0.025 # [m]
r_pebble_core = d_pebble_core/2. # [m]

CHM

The carbon to heavy metal atomic ratio is taken to be about 300. The base case for many previous neutronics calculations (Cisneros) is a carbon to heavy metal ratio of 300.8.


In [19]:
chm = 300.8 # [carbon atoms / heavy metal atoms]

Kernel Diameter

$$d_{kernel} = [400,500] \mu m$$

In [20]:
d_kernel_um = 500 # [micrometers]
d_kernel = d_kernel_um*1e-6 # [m]
r_kernel = d_kernel/2. # [m]

Pebble lattice

The pebbles are in a BCC lattice with 40% flibe volume fraction.


In [21]:
fuel = material.from_atom_frac({'U235': 0.045, 
                                'U238': 0.955, 
                                'O16': 2.0}, 
                               mass=1.0, 
                               density=10.7)

In [22]:
graphite = material.from_atom_frac({'C12': 0.989, 
                                    'C13': 0.011}, 
                                   mass=1.0, 
                                   density=2.267)

salt = material.from_atom_frac({'C12': 0.989, 
                                    'C13': 0.011}, 
                                   mass=1.0, 
                                   density=1.267)

In [16]:
centers = []
centers.append([1,2,3])
centers.append([4,5,6])
centers.append([7,8,9])
centers


Out[16]:
[[1, 2, 3], [4, 5, 6], [7, 8, 9]]

In [17]:
m = mesh.Mesh(structured_coords=[np.linspace(0.0, total_radius, 100*total_radius+1), np.linspace(0.0, height, 100*height+1), [0.0, 1.0]], structured=True)
for i, mat, ve in m:
    m.mats[i] = salt
    
x_vals = np.linspace(0.0,1.5,4)
y_vals = np.linspace(0.0,6.0,13)
z_vals = [0.5]

auto_centroids = []
for x in x_vals:
    for y in y_vals:
        for z in z_vals:
            r=[x,y,z]
            auto_centroids.append(r)
            
centroids = [[0,0,0.5],
             [0,0.5,0.5],
             [0,1.0,0.5],
             [0,1.5,0.5],
             [0,2.0,0.5],
             [0,2.5,0.5],
             [0,3.0,0.5],
             [0,3.5,0.5],
             [0,4.0,0.5],
             [0,4.5,0.5],
             [0,5.0,0.5],
             [0,5.5,0.5],
             [0,6.0,0.5],
             [0.5,0.5,0.5],
             [0.5,1.0,0.5],
             [0.5,1.5,0.5],
             [0.5,2.0,0.5],
             [0.5,2.5,0.5],
             [0.5,3.0,0.5],
             [0.5,3.5,0.5],
             [0.5,4.0,0.5],
             [0.5,4.5,0.5],
             [0.5,5.0,0.5],
             [0.5,5.5,0.5],
             [0.5,6.0,0.5],
             [1.0,0.5,0.5],
             [1.0,1.0,0.5],
             [1.0,1.5,0.5],
             [1.0,2.0,0.5],
             [1.0,2.5,0.5],
             [1.0,3.0,0.5],
             [1.0,3.5,0.5],
             [1.0,4.0,0.5],
             [1.0,4.5,0.5],
             [1.0,5.0,0.5],
             [1.0,5.5,0.5],
             [1.0,6.0,0.5],
             [1.5,0.5,0.5],
             [1.5,1.0,0.5],
             [1.5,1.5,0.5],
             [1.5,2.0,0.5],
             [1.5,2.5,0.5],
             [1.5,3.0,0.5],
             [1.5,3.5,0.5],
             [1.5,4.0,0.5],
             [1.5,4.5,0.5],
             [1.5,5.0,0.5],
             [1.5,5.5,0.5],
             [1.5,6.0,0.5]]

for centroid in centroids :
    pebble = Sphere(0.2, centroid)
    pebble.add_element(m,fuel)
    pebble_core = Sphere(0.1, centroid)
    pebble_core.add_element(m,graphite)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-93cd85eb29c8> in <module>()
----> 1 m = mesh.Mesh(structured_coords=[np.linspace(0.0, total_radius, 100*total_radius+1), np.linspace(0.0, height, 100*height+1), [0.0, 1.0]], structured=True)
      2 for i, mat, ve in m:
      3     m.mats[i] = salt
      4 
      5 x_vals = np.linspace(0.0,1.5,4)

NameError: name 'total_radius' is not defined

In [55]:
#m.sigma_t = mesh.ComputedTag(lambda mesh, i: sigma_t(mesh.mats[i], group_struct=[1.0, 1e-6], phi_g=[1.0])[0])
#m.u_mass = mesh.ComputedTag(lambda mesh, i: max(mesh.mats[i]['Th':'Pu'].mass, 0.0))
pf = PyneMoabHex8StaticOutput(m)
SlicePlot(pf, 'z', 'density', origin='native').display()
#print "Pu Content of Mesh:", sum(m.u_mass[:])




In [25]:


In [44]: