Modeling a CANDU BundleΒΆ

candu

In this example, we will create a typical CANDU bundle with rings of fuel pins. At present, OpenMC does not have a specialized lattice for this type of fuel arrangement, so we must resort to manual creation of the array of fuel pins.

In [1]:
%matplotlib inline
from math import pi, sin, cos
import numpy as np
import openmc

Let's begin by creating the materials that will be used in our model.

In [2]:
fuel = openmc.Material(name='fuel')
fuel.add_element('U', 1.0)
fuel.add_element('O', 2.0)
fuel.set_density('g/cm3', 10.0)

clad = openmc.Material(name='zircaloy')
clad.add_element('Zr', 1.0)
clad.set_density('g/cm3', 6.0)

heavy_water = openmc.Material(name='heavy water')
heavy_water.add_nuclide('H2', 2.0)
heavy_water.add_nuclide('O16', 1.0)
heavy_water.add_s_alpha_beta('c_D_in_D2O')
heavy_water.set_density('g/cm3', 1.1)

With out materials created, we'll now define key dimensions in our model. These dimensions are taken from the example in section 11.1.3 of the Serpent manual.

In [3]:
# Outer radius of fuel and clad
r_fuel = 0.6122
r_clad = 0.6540

# Pressure tube and calendria radii
pressure_tube_ir = 5.16890
pressure_tube_or = 5.60320
calendria_ir = 6.44780
calendria_or = 6.58750

# Radius to center of each ring of fuel pins
ring_radii = np.array([0.0, 1.4885, 2.8755, 4.3305])

To begin creating the bundle, we'll first create annular regions completely filled with heavy water and add in the fuel pins later. The radii that we've specified above correspond to the center of each ring. We actually need to create cylindrical surfaces at radii that are half-way between the centers.

In [4]:
# These are the surfaces that will divide each of the rings
radial_surf = [openmc.ZCylinder(R=r) for r in
               (ring_radii[:-1] + ring_radii[1:])/2]

water_cells = []
for i in range(ring_radii.size):
    # Create annular region
    if i == 0:
        water_region = -radial_surf[i]
    elif i == ring_radii.size - 1:
        water_region = +radial_surf[i-1]
    else:
        water_region = +radial_surf[i-1] & -radial_surf[i]
        
    water_cells.append(openmc.Cell(fill=heavy_water, region=water_region))

Let's see what our geometry looks like so far. In order to plot the geometry, we create a universe that contains the annular water cells and then use the Universe.plot() method. While we're at it, we'll set some keyword arguments that can be reused for later plots.

In [5]:
plot_args = {'width': (2*calendria_or, 2*calendria_or)}
bundle_universe = openmc.Universe(cells=water_cells)
bundle_universe.plot(**plot_args)

Now we need to create a universe that contains a fuel pin. Note that we don't actually need to put water outside of the cladding in this universe because it will be truncated by a higher universe.

In [6]:
surf_fuel = openmc.ZCylinder(R=r_fuel)

fuel_cell = openmc.Cell(fill=fuel, region=-surf_fuel)
clad_cell = openmc.Cell(fill=clad, region=+surf_fuel)

pin_universe = openmc.Universe(cells=(fuel_cell, clad_cell))
In [7]:
pin_universe.plot(**plot_args)

The code below works through each ring to create a cell containing the fuel pin universe. As each fuel pin is created, we modify the region of the water cell to include everything outside the fuel pin.

In [8]:
num_pins = [1, 6, 12, 18]
angles = [0, 0, 15, 0]

for i, (r, n, a) in enumerate(zip(ring_radii, num_pins, angles)):
    for j in range(n):
        # Determine location of center of pin
        theta = (a + j/n*360.) * pi/180.
        x = r*cos(theta)
        y = r*sin(theta)
        
        pin_boundary = openmc.ZCylinder(x0=x, y0=y, R=r_clad)
        water_cells[i].region &= +pin_boundary
        
        # Create each fuel pin -- note that we explicitly assign an ID so 
        # that we can identify the pin later when looking at tallies
        pin = openmc.Cell(fill=pin_universe, region=-pin_boundary)
        pin.translation = (x, y, 0)
        pin.id = (i + 1)*100 + j
        bundle_universe.add_cell(pin)
In [9]:
bundle_universe.plot(**plot_args)

Looking pretty good! Finally, we create cells for the pressure tube and calendria and then put our bundle in the middle of the pressure tube.

In [10]:
pt_inner = openmc.ZCylinder(R=pressure_tube_ir)
pt_outer = openmc.ZCylinder(R=pressure_tube_or)
calendria_inner = openmc.ZCylinder(R=calendria_ir)
calendria_outer = openmc.ZCylinder(R=calendria_or, boundary_type='vacuum')

bundle = openmc.Cell(fill=bundle_universe, region=-pt_inner)
pressure_tube = openmc.Cell(fill=clad, region=+pt_inner & -pt_outer)
v1 = openmc.Cell(region=+pt_outer & -calendria_inner)
calendria = openmc.Cell(fill=clad, region=+calendria_inner & -calendria_outer)

root_universe = openmc.Universe(cells=[bundle, pressure_tube, v1, calendria])

Let's look at the final product. We'll export our geometry and materials and then use plot_inline() to get a nice-looking plot.

In [11]:
geom = openmc.Geometry(root_universe)
geom.export_to_xml()

mats = openmc.Materials(geom.get_all_materials().values())
mats.export_to_xml()
In [12]:
p = openmc.Plot.from_geometry(geom)
p.color_by = 'material'
p.colors = {
    fuel: 'black',
    clad: 'silver',
    heavy_water: 'blue'
}
openmc.plot_inline(p)

Interpreting Results

One of the difficulties of a geometry like this is identifying tally results when there was no lattice involved. To address this, we specifically gave an ID to each fuel pin of the form 100*ring + azimuthal position. Consequently, we can use a distribcell tally and then look at our DataFrame which will show these cell IDs.

In [13]:
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 10
settings.source = openmc.Source(space=openmc.stats.Point())
settings.export_to_xml()
In [14]:
fuel_tally = openmc.Tally()
fuel_tally.filters = [openmc.DistribcellFilter(fuel_cell)]
fuel_tally.scores = ['flux']

tallies = openmc.Tallies([fuel_tally])
tallies.export_to_xml()
In [15]:
openmc.run(output=False)
Out[15]:
0

The return code of 0 indicates that OpenMC ran successfully. Now let's load the statepoint into a openmc.StatePoint object and use the Tally.get_pandas_dataframe(...) method to see our results.

In [16]:
sp = openmc.StatePoint('statepoint.{}.h5'.format(settings.batches))
In [17]:
t = sp.get_tally()
t.get_pandas_dataframe()
Out[17]:
level 1 level 2 level 3 distribcell nuclide score mean std. dev.
univ cell univ cell univ cell
id id id id id id
0 10002 10043 10000 100 10001 10004 0 total flux 0.207805 0.007037
1 10002 10043 10000 200 10001 10004 1 total flux 0.197197 0.005272
2 10002 10043 10000 201 10001 10004 2 total flux 0.190310 0.007816
3 10002 10043 10000 202 10001 10004 3 total flux 0.194736 0.006469
4 10002 10043 10000 203 10001 10004 4 total flux 0.191097 0.006431
5 10002 10043 10000 204 10001 10004 5 total flux 0.189910 0.004891
6 10002 10043 10000 205 10001 10004 6 total flux 0.182203 0.003851
7 10002 10043 10000 300 10001 10004 7 total flux 0.165922 0.005815
8 10002 10043 10000 301 10001 10004 8 total flux 0.168933 0.008300
9 10002 10043 10000 302 10001 10004 9 total flux 0.159587 0.003085
10 10002 10043 10000 303 10001 10004 10 total flux 0.159158 0.005910
11 10002 10043 10000 304 10001 10004 11 total flux 0.148537 0.005308
12 10002 10043 10000 305 10001 10004 12 total flux 0.150945 0.006654
13 10002 10043 10000 306 10001 10004 13 total flux 0.154237 0.003665
14 10002 10043 10000 307 10001 10004 14 total flux 0.165888 0.004733
15 10002 10043 10000 308 10001 10004 15 total flux 0.156777 0.006540
16 10002 10043 10000 309 10001 10004 16 total flux 0.165277 0.005935
17 10002 10043 10000 310 10001 10004 17 total flux 0.156528 0.005732
18 10002 10043 10000 311 10001 10004 18 total flux 0.159610 0.004584
19 10002 10043 10000 400 10001 10004 19 total flux 0.096597 0.004466
20 10002 10043 10000 401 10001 10004 20 total flux 0.118214 0.005451
21 10002 10043 10000 402 10001 10004 21 total flux 0.106167 0.004722
22 10002 10043 10000 403 10001 10004 22 total flux 0.110814 0.004208
23 10002 10043 10000 404 10001 10004 23 total flux 0.112319 0.005079
24 10002 10043 10000 405 10001 10004 24 total flux 0.110232 0.004153
25 10002 10043 10000 406 10001 10004 25 total flux 0.099967 0.005085
26 10002 10043 10000 407 10001 10004 26 total flux 0.095444 0.003615
27 10002 10043 10000 408 10001 10004 27 total flux 0.092620 0.003997
28 10002 10043 10000 409 10001 10004 28 total flux 0.095517 0.004022
29 10002 10043 10000 410 10001 10004 29 total flux 0.113737 0.009530
30 10002 10043 10000 411 10001 10004 30 total flux 0.108368 0.007241
31 10002 10043 10000 412 10001 10004 31 total flux 0.106990 0.005716
32 10002 10043 10000 413 10001 10004 32 total flux 0.112050 0.005002
33 10002 10043 10000 414 10001 10004 33 total flux 0.115054 0.006239
34 10002 10043 10000 415 10001 10004 34 total flux 0.114394 0.004919
35 10002 10043 10000 416 10001 10004 35 total flux 0.114352 0.005322
36 10002 10043 10000 417 10001 10004 36 total flux 0.110890 0.005051

We can see that in the 'level 2' column, the 'cell id' tells us how each row corresponds to a ring and azimuthal position.