OpenMC Guide
Example: Pin Cell Model
Building Your First Model
A pin cell is the basic repeating unit in a reactor core: a cylindrical fuel pellet with zircaloy cladding in a square water moderator lattice, representative of a pressurized water reactor (PWR). This example covers the full OpenMC workflow: materials, geometry, settings, tallies, and results.
┌──────────────────────────┐ │ │ pitch = 1.26 cm │ ┌──────────────────┐ │ │ │ │ │ Clad: r = 0.47 cm │ │ ┌──────────┐ │ │ │ │ │ │ │ │ Fuel: r = 0.41 cm │ │ │ UO₂ │ │ │ │ │ │ │ │ │ │ │ └──────────┘ │ │ │ │ Cladding │ │ │ └──────────────────┘ │ │ Water │ └──────────────────────────┘
Complete Model Code
Complete working model — copy and run directly:
python
import openmc
import numpy as np
# =============================================================================
# MATERIALS
# =============================================================================
# UO2 fuel with 4.5% enrichment
fuel = openmc.Material(name='UO2 Fuel')
fuel.set_density('g/cm3', 10.4)
fuel.add_nuclide('U235', 0.045) # 4.5% enriched
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 with thermal scattering
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') # Thermal scattering
# Create materials collection
materials = openmc.Materials([fuel, clad, water])
# =============================================================================
# GEOMETRY
# =============================================================================
# Define surfaces
fuel_radius = 0.4096 # cm
clad_radius = 0.4750 # cm
pitch = 1.26 # cm (square lattice)
fuel_outer = openmc.ZCylinder(r=fuel_radius)
clad_outer = openmc.ZCylinder(r=clad_radius)
# Create square boundary
left = openmc.XPlane(-pitch/2, boundary_type='reflective')
right = openmc.XPlane(pitch/2, boundary_type='reflective')
bottom = openmc.YPlane(-pitch/2, boundary_type='reflective')
top = openmc.YPlane(pitch/2, boundary_type='reflective')
# Define regions using boolean operations
fuel_region = -fuel_outer
clad_region = +fuel_outer & -clad_outer
water_region = +clad_outer & +left & -right & +bottom & -top
# Create cells
fuel_cell = openmc.Cell(fill=fuel, region=fuel_region)
clad_cell = openmc.Cell(fill=clad, region=clad_region)
water_cell = openmc.Cell(fill=water, region=water_region)
# Create geometry
universe = openmc.Universe(cells=[fuel_cell, clad_cell, water_cell])
geometry = openmc.Geometry(universe)
# =============================================================================
# SIMULATION SETTINGS
# =============================================================================
settings = openmc.Settings()
settings.particles = 10000
settings.batches = 100
settings.inactive = 20
# Set neutron source in the fuel region
source_region = openmc.stats.Box(
[-fuel_radius, -fuel_radius, -1],
[fuel_radius, fuel_radius, 1]
)
settings.source = openmc.IndependentSource(space=source_region)
# =============================================================================
# TALLIES
# =============================================================================
# Create tallies to get useful results
tallies = openmc.Tallies()
# Flux in each region
fuel_tally = openmc.Tally(name='fuel flux')
fuel_tally.filters = [openmc.CellFilter(fuel_cell)]
fuel_tally.scores = ['flux', 'nu-fission', 'absorption']
tallies.append(fuel_tally)
clad_tally = openmc.Tally(name='clad flux')
clad_tally.filters = [openmc.CellFilter(clad_cell)]
clad_tally.scores = ['flux', 'absorption']
tallies.append(clad_tally)
water_tally = openmc.Tally(name='water flux')
water_tally.filters = [openmc.CellFilter(water_cell)]
water_tally.scores = ['flux', 'absorption']
tallies.append(water_tally)
# =============================================================================
# RUN SIMULATION
# =============================================================================
# Export model files and run
model = openmc.Model(geometry, materials, settings, tallies)
# Run OpenMC (returns path to statepoint file)
statepoint_path = model.run()
# Open statepoint and extract k-effective
sp = openmc.StatePoint(statepoint_path)
keff = sp.keff
print(f"k-effective: {keff.nominal_value:.5f} ± {keff.std_dev:.5f}")
# =============================================================================
# ANALYZE RESULTS
# =============================================================================
# Load results
sp = openmc.StatePoint('statepoint.100.h5')
# Extract tally data
fuel_flux = sp.get_tally(name='fuel flux')
clad_flux = sp.get_tally(name='clad flux')
water_flux = sp.get_tally(name='water flux')
print("\nFlux Results:")
print(f"Fuel flux: {fuel_flux.mean[0, 0, 0]:.3e} ± {fuel_flux.std_dev[0, 0, 0]:.3e}")
print(f"Clad flux: {clad_flux.mean[0, 0, 0]:.3e} ± {clad_flux.std_dev[0, 0, 0]:.3e}")
print(f"Water flux: {water_flux.mean[0, 0, 0]:.3e} ± {water_flux.std_dev[0, 0, 0]:.3e}")
# Calculate reaction rates
nu_fission_rate = fuel_flux.get_slice(scores=['nu-fission'])
absorption_rate = fuel_flux.get_slice(scores=['absorption'])
print("\nReaction Rates in Fuel:")
print(f"Nu-fission rate: {nu_fission_rate.mean[0, 0, 0]:.3e}")
print(f"Absorption rate: {absorption_rate.mean[0, 0, 0]:.3e}")
print(f"Eta (η = νΣ_f/Σ_a): {nu_fission_rate.mean[0, 0, 0] / absorption_rate.mean[0, 0, 0]:.3f}")
print("\nSimulation completed successfully!")Mirrors runnable script in examples repo
Understanding the Results
The k-effective indicates whether the system is critical (k=1), subcritical (k<1), or supercritical (k>1). The reproduction factor η = νΣ_f/Σ_a gives the average number of fission neutrons produced per neutron absorbed in fuel.
Expected Results
OpenMC tallies are normalized per source particle, so the absolute values of flux and reaction rates depend on normalization. The relative ratios between regions and the derived quantities like η are the physically meaningful results.
bash
# Typical output for this model (values are per source particle):
k-effective: 1.18234 ± 0.00234
Flux Results:
Fuel flux: 2.46e+01 ± 1.23e-01
Clad flux: 1.79e+01 ± 9.88e-02
Water flux: 8.90e+00 ± 4.57e-02
Reaction Rates in Fuel:
Nu-fission rate: 5.43e-02
Absorption rate: 2.99e-02
Eta (η = νΣ_f/Σ_a): 1.818Tutorial snippet — no separate file in examples repo
Exploring Further
Try varying parameters to see their effect on k-effective:
Parameter Studies
python
# Try different enrichments
for enrichment in [0.03, 0.035, 0.04, 0.045, 0.05]:
fuel.remove_nuclide('U235')
fuel.remove_nuclide('U238')
fuel.add_nuclide('U235', enrichment)
fuel.add_nuclide('U238', 1 - enrichment)
sp = openmc.StatePoint(model.run())
keff = sp.keff
print(f"Enrichment {enrichment:.1%}: k = {keff.nominal_value:.4f}")
# Try different fuel radii
for radius in [0.35, 0.40, 0.45]:
fuel_outer.r = radius
sp = openmc.StatePoint(model.run())
keff = sp.keff
print(f"Fuel radius {radius} cm: k = {keff.nominal_value:.4f}")
# Try different lattice pitches
for p in [1.2, 1.26, 1.3, 1.4]:
left.x0 = -p/2
right.x0 = p/2
bottom.y0 = -p/2
top.y0 = p/2
sp = openmc.StatePoint(model.run())
keff = sp.keff
print(f"Pitch {p} cm: k = {keff.nominal_value:.4f}")Tutorial snippet — no separate file in examples repo
Learning Tip: As you experiment with parameters, notice how k-effective changes. This will help you understand neutron multiplication and the key factors that control reactor criticality.
Next Steps
The same workflow — materials, geometry, settings, tallies — extends to fuel assemblies, full cores, and advanced physics like depletion.