OpenMC Guide

❮ Example: Pin Cell Example: Full Reactor ❯

Example: Fuel Assembly

Building an Assembly

This example builds on the pin cell concepts to create a complete 17×17 PWR fuel assembly. We'll use universes and lattices to manage the complexity while keeping the code organized and maintainable.

The approach is hierarchical: create pin universes, arrange them in a lattice, then place the lattice in an assembly cell. This pattern scales to full core models.

Materials Definition

python

import openmc
import numpy as np

# =============================================================================
# MATERIALS
# =============================================================================

# UO2 fuel - 4.5% enriched
fuel = openmc.Material(name='UO2 Fuel')
fuel.set_density('g/cm3', 10.4)
fuel.add_nuclide('U235', 0.045)
fuel.add_nuclide('U238', 0.955)
fuel.add_element('O', 2.0)

# Zircaloy-4 cladding
clad = openmc.Material(name='Zircaloy-4')
clad.set_density('g/cm3', 6.55)
clad.add_element('Zr', 0.982)
clad.add_element('Sn', 0.015)
clad.add_element('Fe', 0.002)
clad.add_element('Cr', 0.001)

# Light water moderator
water = openmc.Material(name='Light Water')
water.set_density('g/cm3', 1.0)
water.add_nuclide('H1', 2.0)
water.add_element('O', 1.0)
water.add_s_alpha_beta('c_H_in_H2O')

# Stainless steel for guide tubes
steel = openmc.Material(name='Stainless Steel')
steel.set_density('g/cm3', 8.0)
steel.add_element('Fe', 0.68)
steel.add_element('Cr', 0.20)
steel.add_element('Ni', 0.12)

materials = openmc.Materials([fuel, clad, water, steel])

Pin Universe Definitions

python

# =============================================================================
# PIN UNIVERSES
# =============================================================================

def create_fuel_pin():
    """Standard fuel pin with cladding"""
    fuel_radius = 0.4096
    clad_inner = 0.4178
    clad_outer = 0.4750
    
    # Surfaces
    fuel_surf = openmc.ZCylinder(r=fuel_radius)
    clad_inner_surf = openmc.ZCylinder(r=clad_inner)
    clad_outer_surf = openmc.ZCylinder(r=clad_outer)
    
    # Regions
    fuel_region = -fuel_surf
    gap_region = +fuel_surf & -clad_inner_surf
    clad_region = +clad_inner_surf & -clad_outer_surf
    water_region = +clad_outer_surf
    
    # Cells
    fuel_cell = openmc.Cell(fill=fuel, region=fuel_region)
    gap_cell = openmc.Cell(fill=None, region=gap_region)
    clad_cell = openmc.Cell(fill=clad, region=clad_region)
    water_cell = openmc.Cell(fill=water, region=water_region)
    
    return openmc.Universe(cells=[fuel_cell, gap_cell, clad_cell, water_cell])

def create_guide_tube():
    """Guide tube for control rods"""
    guide_inner = 0.5610
    guide_outer = 0.6020
    
    # Surfaces
    guide_inner_surf = openmc.ZCylinder(r=guide_inner)
    guide_outer_surf = openmc.ZCylinder(r=guide_outer)
    
    # Regions
    water_inner = -guide_inner_surf
    tube_region = +guide_inner_surf & -guide_outer_surf
    water_outer = +guide_outer_surf
    
    # Cells
    water_inner_cell = openmc.Cell(fill=water, region=water_inner)
    tube_cell = openmc.Cell(fill=steel, region=tube_region)
    water_outer_cell = openmc.Cell(fill=water, region=water_outer)
    
    return openmc.Universe(cells=[water_inner_cell, tube_cell, water_outer_cell])

def create_instrumentation_tube():
    """Instrumentation tube (water-filled)"""
    tube_inner = 0.5590
    tube_outer = 0.6050
    
    # Surfaces
    tube_inner_surf = openmc.ZCylinder(r=tube_inner)
    tube_outer_surf = openmc.ZCylinder(r=tube_outer)
    
    # Regions
    water_inner = -tube_inner_surf
    tube_region = +tube_inner_surf & -tube_outer_surf
    water_outer = +tube_outer_surf
    
    # Cells
    water_inner_cell = openmc.Cell(fill=water, region=water_inner)
    tube_cell = openmc.Cell(fill=steel, region=tube_region)
    water_outer_cell = openmc.Cell(fill=water, region=water_outer)
    
    return openmc.Universe(cells=[water_inner_cell, tube_cell, water_outer_cell])

# Create the pin universes
fuel_pin = create_fuel_pin()
guide_tube = create_guide_tube()
instr_tube = create_instrumentation_tube()

Assembly Lattice

python

# =============================================================================
# ASSEMBLY LATTICE
# =============================================================================

# Create 17x17 lattice
lattice = openmc.RectLattice(name='Assembly Lattice')
lattice.lower_left = [-10.71, -10.71]  # Half of 17 * 1.26 cm pitch
lattice.pitch = [1.26, 1.26]

# Define assembly pattern
# Standard Westinghouse 17x17 assembly layout
# F = fuel pin, G = guide tube, I = instrumentation tube

assembly_pattern = [
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','F','F','G','F','F','G','F','F','G','F','F','F','F','F'],
    ['F','F','F','G','F','F','F','F','F','F','F','F','F','G','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','G','F','F','G','F','F','G','F','F','G','F','F','G','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','G','F','F','G','F','F','I','F','F','G','F','F','G','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','G','F','F','G','F','F','G','F','F','G','F','F','G','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','G','F','F','F','F','F','F','F','F','F','G','F','F','F'],
    ['F','F','F','F','F','G','F','F','G','F','F','G','F','F','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F'],
    ['F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F','F']
]

# Convert pattern to universe array
universe_map = {'F': fuel_pin, 'G': guide_tube, 'I': instr_tube}
lattice_universes = []
for row in assembly_pattern:
    universe_row = [universe_map[symbol] for symbol in row]
    lattice_universes.append(universe_row)

lattice.universes = lattice_universes

Assembly Cell and Geometry

python

# =============================================================================
# ASSEMBLY GEOMETRY
# =============================================================================

# Assembly boundary (21.42 cm square)
pitch = 1.26
assembly_size = 17 * pitch

# Create assembly boundary
assembly_left = openmc.XPlane(-assembly_size/2, boundary_type='reflective')
assembly_right = openmc.XPlane(assembly_size/2, boundary_type='reflective')
assembly_bottom = openmc.YPlane(-assembly_size/2, boundary_type='reflective')
assembly_top = openmc.YPlane(assembly_size/2, boundary_type='reflective')
assembly_lower = openmc.ZPlane(-200, boundary_type='vacuum')
assembly_upper = openmc.ZPlane(200, boundary_type='vacuum')

# Assembly region
assembly_region = (+assembly_left & -assembly_right & 
                  +assembly_bottom & -assembly_top &
                  +assembly_lower & -assembly_upper)

# Assembly cell
assembly_cell = openmc.Cell(fill=lattice, region=assembly_region)

# Root universe and geometry
root_universe = openmc.Universe(cells=[assembly_cell])
geometry = openmc.Geometry(root_universe)

Simulation Settings

python

# =============================================================================
# SIMULATION SETTINGS
# =============================================================================

settings = openmc.Settings()
settings.particles = 10000
settings.batches = 120
settings.inactive = 20

# Source distribution across the assembly
source_box = openmc.stats.Box(
    [-assembly_size/2, -assembly_size/2, -100],
    [assembly_size/2, assembly_size/2, 100]
)
settings.source = openmc.Source(space=source_box)

# Set temperature (PWR operating conditions)
fuel.temperature = 900   # K
clad.temperature = 600   # K  
water.temperature = 574  # K
steel.temperature = 574  # K

Tallies

python

# =============================================================================
# TALLIES
# =============================================================================

tallies = openmc.Tallies()

# Assembly-averaged flux and reaction rates
assembly_tally = openmc.Tally(name='assembly_avg')
assembly_tally.filters = [openmc.CellFilter(assembly_cell)]
assembly_tally.scores = ['flux', 'fission', 'absorption', 'nu-fission']
tallies.append(assembly_tally)

# Pin-by-pin power distribution
mesh = openmc.RegularMesh()
mesh.lower_left = [-assembly_size/2, -assembly_size/2]
mesh.upper_right = [assembly_size/2, assembly_size/2]
mesh.dimension = [17, 17]  # One bin per pin

mesh_tally = openmc.Tally(name='pin_powers')
mesh_tally.filters = [openmc.MeshFilter(mesh)]
mesh_tally.scores = ['kappa-fission']  # Energy released per fission
tallies.append(mesh_tally)

# Energy spectrum in fuel regions
energy_bins = np.logspace(-3, 7, 50)  # 1 meV to 10 MeV
energy_filter = openmc.EnergyFilter(energy_bins)

spectrum_tally = openmc.Tally(name='fuel_spectrum')
spectrum_tally.filters = [openmc.MaterialFilter(fuel), energy_filter]
spectrum_tally.scores = ['flux']
tallies.append(spectrum_tally)

Running the Model

python

# =============================================================================
# RUN SIMULATION
# =============================================================================

# Create complete model
model = openmc.Model(geometry, materials, settings, tallies)

# Run simulation (model.run() returns path to statepoint)
statepoint_path = model.run()
sp = openmc.StatePoint(statepoint_path)
keff = sp.keff
print(f"k-effective: {keff.nominal_value:.5f}")
print(f"Standard deviation: {keff.std_dev:.5f}")

# =============================================================================
# BASIC ANALYSIS
# =============================================================================

# Load results
sp = openmc.StatePoint('statepoint.120.h5')

# Get pin power distribution
pin_powers = sp.get_tally(name='pin_powers')
power_array = pin_powers.mean.reshape((17, 17))

# Simple power peaking analysis
max_power = np.max(power_array)
avg_power = np.mean(power_array)
peaking_factor = max_power / avg_power

print(f"Peak-to-average power ratio: {peaking_factor:.3f}")

# Find location of peak power
peak_location = np.unravel_index(np.argmax(power_array), power_array.shape)
print(f"Peak power at pin ({peak_location[0]+1}, {peak_location[1]+1})")

This complete assembly model demonstrates the hierarchical approach that makes OpenMC powerful for reactor modeling. The same pattern extends to full core models by treating assemblies as the basic building blocks.

OpenMC Guide

Example: Fuel Assembly

Building an Assembly

Materials Definition

Pin Universe Definitions

Assembly Lattice

Assembly Cell and Geometry

Simulation Settings

Tallies

Running the Model

Related Topics

Prerequisite

Example: Pin Cell

Next Steps

Example: Full Reactor

Related

OpenMC Guide Overview

Monte Carlo Basics