OpenMC Guide
Example: Pin Cell Model
Building Your First Model
Let's build a complete fuel pin cell model from scratch. This example will demonstrate the essential OpenMC workflow: define materials, create geometry, set up the simulation, and analyze results.
A pin cell represents the basic repeating unit in a reactor core. Our model includes a cylindrical fuel pellet, zircaloy cladding, and water moderator in a square lattice. This is representative of what you'd find in a pressurized water reactor (PWR).
┌──────────────────────────┐ │ │ pitch = 1.26 cm │ ┌──────────────────┐ │ │ │ │ │ Clad: r = 0.47 cm │ │ ┌──────────┐ │ │ │ │ │ │ │ │ Fuel: r = 0.41 cm │ │ │ UO₂ │ │ │ │ │ │ │ │ │ │ │ └──────────┘ │ │ │ │ Cladding │ │ │ └──────────────────┘ │ │ Water │ └──────────────────────────┘
Complete Model Code
Here's the complete working model. You can copy this code and run it directly - it includes everything needed for a basic pin cell calculation.
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.Source(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', '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
fission_rate = fuel_flux.get_slice(scores=['fission'])
absorption_rate = fuel_flux.get_slice(scores=['absorption'])
print("\nReaction Rates in Fuel:")
print(f"Fission rate: {fission_rate.mean[0, 0, 0]:.3e}")
print(f"Absorption rate: {absorption_rate.mean[0, 0, 0]:.3e}")
print(f"Eta (η): {fission_rate.mean[0, 0, 0] / absorption_rate.mean[0, 0, 0]:.3f}")
print("\nSimulation completed successfully!")Understanding the Results
When you run this model, you'll get several key results that tell you about the neutron physics. The k-effective tells you if your system is critical (k=1), subcritical (k<1), or supercritical (k>1).
The flux results show how neutrons are distributed in different regions. Typically, you'll see the highest flux in the fuel region where fission reactions occur. The eta (η) value represents the reproduction factor - the average number of fission neutrons produced per thermal neutron absorbed in fuel.
Expected Results
# Typical output for this model:
k-effective: 1.18234 ± 0.00234
Flux Results:
Fuel flux: 2.456e+13 ± 1.234e+11
Clad flux: 1.789e+13 ± 9.876e+10
Water flux: 8.901e+12 ± 4.567e+10
Reaction Rates in Fuel:
Fission rate: 5.432e+11
Absorption rate: 2.987e+11
Eta (η): 1.818Exploring Further
Now that you have a working model, try modifying it to understand how different parameters affect the results. Here are some exercises to build your intuition:
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}")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
This pin cell model demonstrates the fundamental OpenMC workflow. From here, you can explore more complex geometries like fuel assemblies, add more sophisticated tallies, or investigate advanced physics like depletion calculations.
The key concepts you've learned - materials, geometry, settings, and tallies - form the foundation for all OpenMC models. As you tackle more complex problems, you'll use these same building blocks in more sophisticated ways.