OpenMC Guide
Example: Pin Cell Model
What you'll learn
- Assemble a complete, runnable PWR pin cell from materials through to results.
- Tally flux and reaction rates separately in the fuel, cladding, and moderator.
- Sanity-check k-effective and the fuel reaction-rate ratio before believing either.
- Vary enrichment and pin radius to see how each moves k-effective.
Before you start
Building Your First Model
This is the payoff of Path A. A pin cell is the basic repeating unit in a PWR lattice: a cylindrical fuel pellet, cladding, and square water moderator. The materials, radii, and API patterns from the previous pages assemble here into one runnable model.
┌──────────────────────────┐ │ │ pitch = 1.26 cm │ ┌──────────────────┐ │ │ │ │ │ Clad: r = 0.4750 cm │ │ ┌──────────┐ │ │ │ │ │ │ │ │ Fuel: r = 0.4096 cm │ │ │ UO₂ │ │ │ │ │ │ │ │ │ │ │ └──────────┘ │ │ │ │ Cladding │ │ │ └──────────────────┘ │ │ Water │ └──────────────────────────┘
The clad cylinder starts right at the pellet surface — this model neglects the thin helium-filled gap that sits between fuel and cladding in a real pin.
Complete Model Code
Complete working model — copy and run directly:
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.56)
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', 0.998)
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
# =============================================================================
# Reuse the statepoint opened above. Hardcoding 'statepoint.100.h5' works
# only as long as batches stays at 100 - the path from model.run() always
# points at the file that was actually written.
# 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}")
# There is no energy filter on this tally, so this ratio is integrated over the
# whole spectrum. It is not the four-factor eta, which is a thermal-group
# quantity - see "Understanding the Results" below.
print(f"nu-fission / absorption in fuel: {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 second quantity, νΣ_f/Σ_a in the fuel, is the average number of fission neutrons produced per neutron absorbed in the fuel.
It is tempting to call that number η and compare it to the textbook value, but the tally above has a cell filter and no energy filter, so it is integrated over the whole spectrum. The four-factor η is defined on the thermal group alone. The spectrum-averaged version you get here is substantially lower, because it includes every epithermal capture in the U-238 resonances — captures that the thermal-group η excludes by construction. The two are different quantities and disagreeing with each other is not a bug. Add an EnergyFilter with a thermal cutoff if you want the four-factor quantity.
Reading the output
OpenMC tallies are normalized per source particle, so absolute flux values carry no physical meaning on their own. What you should check are the relationships between regions and the derived ratios. Your exact numbers will depend on which nuclear data library OPENMC_CROSS_SECTIONS points at, so treat the shape of the output as the thing to match, not the digits.
k-effective: <value> ± <uncertainty>
Flux Results:
Fuel flux: <value> ± <uncertainty>
Clad flux: <value> ± <uncertainty>
Water flux: <value> ± <uncertainty>
Reaction Rates in Fuel:
Nu-fission rate: <value>
Absorption rate: <value>
nu-fission / absorption in fuel: <value>Tutorial snippet — no separate file in examples repo
Sanity checks on your own run
Four checks are worth running before you trust any number from this model, and each catches a different class of mistake. The first is the easiest: k-effective should comfortably exceed 1. This is a cold, unborated, 4.5 % enriched lattice with reflective boundaries, so there is no leakage and no soluble poison holding it down. A value near or below 1 means something is wrong with the fuel or the moderator, not with the statistics.
The second you can derive rather than look up: νΣf/Σa in the fuel should come out above k-effective. Nothing escapes this model — the radial boundary is reflective and the pin is an infinite cylinder with no axial boundary at all — so every neutron produced is eventually absorbed somewhere, which makes k-effective essentially total production over total absorption. The tallied ratio has the same numerator but counts only the absorptions in the fuel, and dropping the clad and water from a denominator can only raise the quotient. Dividing k-effective by the tallied ratio recovers the fraction of all absorptions that happen in the fuel. A tallied ratio that comes out below k-effective means the tally is not measuring what you think it is.
Two caveats mark where that derivation stops being exact. OpenMC's absorption score covers every reaction that produces no secondary neutrons, plus fission, which means it excludes (n,2n) and (n,3n) — those score as scattering, so the production-equals-absorption balance is off by that small surplus. And the fuel-absorption fraction you recover is integrated over all energies, so it is not the four-factor thermal utilization f, for the same reason the ratio above is not η.
The third check is about arithmetic rather than physics: divide each flux by its region volume before comparing anything. The cladding is a thin annulus between the fuel and the water, so its flux density belongs between theirs rather than above both, and comparing the raw volume-integrated tallies instead is the single most common mistake on this page. Fourth and last, the relative errors should sit below about 1 %. Uncertainty falls as one over the square root of the total active histories, so more particles per batch and more active batches both work — particles per batch is usually the better lever, since it parallelizes cleanly and adds no correlated generations.
Exploring Further
Every object in the model above is still a live Python variable, so a parameter study is just a loop that mutates one of them and reruns — no new files, no re-parsing a deck. Before running any of the three loops below, predict the direction k-effective will move for each one; the reasoning behind each direction is worth more than the printed number.
Parameter Studies
# 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
None of these are leakage effects — every side of this model is reflective and the pin is axially infinite, so nothing escapes. What the loops move is the balance between fissile inventory, how thoroughly the water thermalises neutrons before they return to the fuel, and how much U-238 resonance capture they meet on the way. Enrichment is the straightforward one. Radius and pitch both change the water-to-fuel ratio, which has an optimum rather than a direction: a PWR lattice at this pitch sits deliberately on the under-moderated side of it, so that losing coolant costs reactivity.
Check yourself
- Assemble a complete, runnable PWR pin cell from materials through to results?
- Tally flux and reaction rates separately in the fuel, cladding, and moderator?
- Sanity-check k-effective and the fuel reaction-rate ratio before believing either?
- Vary enrichment and pin radius, and predict which way each moves k-effective?