🔥 Firefly Reading Circle — Build a Nuclear Reactor in OpenMC¶
Goal (1 hour): build, run, and visualize a simplified nuclear reactor model, understand k-effective, and experiment with what makes a reactor more or less reactive.
Run each cell top-to-bottom with Shift+Enter. Every explanation lives in its own markdown (text) cell, right above the code it describes — code cells are kept short and mostly comment-free so the two stay visually distinct as you scroll.
Companion doc: READING_CIRCLE_GUIDE.md in this folder has the full
facilitator run-of-show, timing, and discussion prompts.
0. Setup check¶
Run this first. It loads the libraries and points OpenMC at the nuclear
data that ships with this folder (xs_data/), so there's nothing else to
configure. It should print the OpenMC version and "ready".
In [1]:
import os
import glob
from pathlib import Path
XS = Path("xs_data/cross_sections.xml").resolve()
os.environ["OPENMC_CROSS_SECTIONS"] = str(XS)
import openmc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
openmc.config["cross_sections"] = XS
In [2]:
print("OpenMC version:", openmc.__version__)
print("Nuclear data :", XS, "✅ ready" if XS.exists() else "❌ MISSING — open this notebook from the project folder")
OpenMC version: 0.16.0 Nuclear data : /home/samarth/work/openmc-reactor-circle/xs_data/cross_sections.xml ✅ ready
1. The big idea¶
A nuclear reactor is just three ingredients in careful balance:
- Fuel — releases neutrons when it fissions. We use UO2 enriched to 4% U-235 (typical PWR fuel).
- Moderator — slows fast neutrons down so they're more likely to cause another fission. We use ordinary (lightly borated) water.
- Control rods — deliberately absorb neutrons to control the chain reaction. We use boron carbide (B4C). Push them in to slow the reaction down, pull them out to speed it up.
The chain reaction is summarized by one number, k-effective (keff): the average number of neutrons from one fission that go on to cause another fission.
| keff | meaning |
|---|---|
| = 1.000 | critical — steady, self-sustaining chain reaction (what a running plant targets) |
| > 1.000 | supercritical — neutron population (and power) growing |
| < 1.000 | subcritical — neutron population (and power) shrinking |
OpenMC is an open-source Monte Carlo particle transport code. It literally simulates many thousands of individual neutrons bouncing through the geometry we define, tracking every scatter/absorption/fission event, and from that statistics estimates keff and where the power is produced.
No code in this section — just read, then continue to Section 2.
In [3]:
fuel = openmc.Material(name="UO2 fuel (4.0% enriched)")
fuel.add_element('U', 1.0, enrichment=4.0)
fuel.add_element('O', 2.0)
fuel.set_density('g/cm3', 10.4)
2b. Cladding — Zircaloy¶
Simplified to pure zirconium for clarity in a teaching model.
In [4]:
clad = openmc.Material(name="Zircaloy cladding")
clad.add_element('Zr', 1.0)
clad.set_density('g/cm3', 6.55)
2c. Moderator/coolant — borated water¶
~500 ppm natural boron and a hot PWR operating density (0.74 g/cm³) are both realistic plant values. Boron here is a deliberate neutron poison — "experiment c" removes it.
In [5]:
water = openmc.Material(name="Borated water moderator")
water.add_element('H', 2.0)
water.add_element('O', 1.0)
water.add_element('B', 500e-6)
water.set_density('g/cm3', 0.74)
water.add_s_alpha_beta('c_H_in_H2O')
2d. Control rod absorber — boron carbide (B4C)¶
Deliberately a strong neutron absorber.
In [6]:
b4c = openmc.Material(name="B4C control rod absorber")
b4c.add_element('B', 4.0)
b4c.add_element('C', 1.0)
b4c.set_density('g/cm3', 2.52)
2e. Air — withdrawn control-rod gas gap¶
Fills the thin gap inside the pin cladding.
In [7]:
air = openmc.Material(name="Air (void/instrument channel)")
air.add_element('N', 0.78)
air.add_element('O', 0.21)
air.add_element('Ar', 0.01)
air.set_density('g/cm3', 0.001225)
2f. Collect and export¶
OpenMC reads materials from an XML file it writes from this Python object.
In [8]:
materials = openmc.Materials([fuel, clad, water, b4c, air])
materials.export_to_xml()
print("Materials defined:", [m.name for m in materials])
Materials defined: ['UO2 fuel (4.0% enriched)', 'Zircaloy cladding', 'Borated water moderator', 'B4C control rod absorber', 'Air (void/instrument channel)']
3. Geometry: pin cell → lattice → finite core¶
We build a single fuel pin cell (concentric cylinders: fuel pellet → gas gap → zircaloy clad → water), then repeat it in a 7×7 lattice to make a mini assembly. A cross of control-rod pins runs through the middle, with a couple of water-filled guide tubes. The whole thing sits inside a water reflector with a vacuum boundary — meaning neutrons that escape the edge are gone for good. This is a real finite reactor model, not an idealized infinite one.
3a. Pin-cell dimensions and cylinders¶
Standard PWR pin pitch and radii (cm).
In [9]:
pin_pitch = 1.26
fuel_or = 0.39
clad_ir = 0.40
clad_or = 0.46
fuel_cyl = openmc.ZCylinder(r=fuel_or)
clad_in_cyl = openmc.ZCylinder(r=clad_ir)
clad_out_cyl = openmc.ZCylinder(r=clad_or)
3b. Three pin-cell "universes"¶
A fuel pin, a control-rod pin (same clad/water, B4C instead of fuel), and a water guide tube (no fuel or absorber at all).
In [10]:
def make_fuel_pin():
"""Standard UO2 fuel pin cell."""
c_fuel = openmc.Cell(fill=fuel, region=-fuel_cyl)
c_gap = openmc.Cell(fill=air, region=+fuel_cyl & -clad_in_cyl)
c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)
c_water = openmc.Cell(fill=water, region=+clad_out_cyl)
return openmc.Universe(cells=[c_fuel, c_gap, c_clad, c_water])
def make_control_pin():
"""Control rod pin: B4C absorber instead of fuel, same clad/water."""
c_b4c = openmc.Cell(fill=b4c, region=-fuel_cyl)
c_gap = openmc.Cell(fill=air, region=+fuel_cyl & -clad_in_cyl)
c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)
c_water = openmc.Cell(fill=water, region=+clad_out_cyl)
return openmc.Universe(cells=[c_b4c, c_gap, c_clad, c_water])
def make_water_pin():
"""Pure water pin (guide tube / instrument location)."""
c_water = openmc.Cell(fill=water, region=+clad_out_cyl)
c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)
c_inner = openmc.Cell(fill=water, region=-clad_in_cyl)
return openmc.Universe(cells=[c_inner, c_clad, c_water])
In [11]:
fuel_u = make_fuel_pin()
ctrl_u = make_control_pin()
wtr_u = make_water_pin()
print("Pin-cell universes built: fuel, control rod, water guide tube")
Pin-cell universes built: fuel, control rod, water guide tube
3c. The core pattern¶
F = fuel pin, C = control rod pin, W = water guide tube.
This is the main thing you'll edit in Section 8 to run experiments —
try changing a few Cs to Fs later and see what happens to keff!
In [12]:
N = 7
F, C, W = fuel_u, ctrl_u, wtr_u
pattern = [
[F, F, F, C, F, F, F],
[F, F, F, C, F, F, F],
[F, F, W, C, W, F, F],
[C, C, C, C, C, C, C],
[F, F, W, C, W, F, F],
[F, F, F, C, F, F, F],
[F, F, F, C, F, F, F],
]
lattice = openmc.RectLattice(name="7x7 simplified assembly")
lattice.lower_left = (-N * pin_pitch / 2, -N * pin_pitch / 2)
lattice.pitch = (pin_pitch, pin_pitch)
lattice.universes = pattern
3d. Core boundary, reflector, and axial extent¶
A 5 cm water reflector surrounds the lattice; outside that is vacuum (real leakage). Top/bottom are reflective, approximating an infinitely tall core so this stays a fast 2D-ish problem.
In [13]:
core_half = N * pin_pitch / 2
reflector_half = core_half + 5.0
core_box = openmc.model.RectangularPrism(
width=2 * core_half, height=2 * core_half, boundary_type='transmission'
)
reflector_box = openmc.model.RectangularPrism(
width=2 * reflector_half, height=2 * reflector_half, boundary_type='vacuum'
)
z_min = openmc.ZPlane(z0=-10, boundary_type='reflective')
z_max = openmc.ZPlane(z0=10, boundary_type='reflective')
3e. Assemble and export the geometry¶
In [14]:
core_cell = openmc.Cell(fill=lattice, region=-core_box & +z_min & -z_max)
reflector_cell = openmc.Cell(
fill=water, region=+core_box & -reflector_box & +z_min & -z_max
)
root_universe = openmc.Universe(cells=[core_cell, reflector_cell])
geometry = openmc.Geometry(root_universe)
geometry.export_to_xml()
print(f"Core half-width: {core_half:.2f} cm")
print(f"Reflector half-width: {reflector_half:.2f} cm")
Core half-width: 4.41 cm Reflector half-width: 9.41 cm
4. Settings: a short criticality (k-eigenvalue) run¶
We use a small particle count on purpose (2000 per batch) so this runs in a couple of seconds — perfect for a live demo where you'll rerun it many times. A production calculation would use millions of particles for a tighter uncertainty on keff.
In [15]:
settings = openmc.Settings()
settings.batches = 60
settings.inactive = 15
settings.particles = 2000
settings.run_mode = 'eigenvalue'
bounds = [-core_half, -core_half, -9, core_half, core_half, 9]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.IndependentSource(space=uniform_dist)
settings.export_to_xml()
print("Settings exported: 60 batches (15 inactive), 2000 particles/batch")
Settings exported: 60 batches (15 inactive), 2000 particles/batch
/home/samarth/miniconda3/envs/openmc-env/lib/python3.13/site-packages/openmc/stats/multivariate.py:1115: FutureWarning: The 'only_fissionable' has been deprecated. Use the 'constraints' argument when defining a source instead.
warn("The 'only_fissionable' has been deprecated. Use the "
5. Tallies: where do neutrons go, and where does power come from?¶
We lay a 100×100 mesh over the whole model and ask OpenMC to record flux (how many neutrons pass through each cell) and fission (how much power is produced in each cell).
In [16]:
mesh = openmc.RegularMesh()
mesh.dimension = [100, 100]
mesh.lower_left = [-reflector_half, -reflector_half]
mesh.upper_right = [reflector_half, reflector_half]
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally(name="flux_and_fission")
tally.filters = [mesh_filter]
tally.scores = ['flux', 'fission']
tallies = openmc.Tallies([tally])
tallies.export_to_xml()
print("Tallies exported")
Tallies exported
6. Define the geometry plots¶
Two 2D slices: one colored by material, one by cell — useful for checking the geometry looks the way you expect before spending time on a full transport run.
In [17]:
plot_xy = openmc.Plot()
plot_xy.filename = 'plot_xy_materials'
plot_xy.width = (2 * reflector_half + 2, 2 * reflector_half + 2)
plot_xy.pixels = (1000, 1000)
plot_xy.color_by = 'material'
plot_xy.colors = {
fuel: (255, 60, 60),
clad: (150, 150, 150),
water: (120, 190, 255),
b4c: (40, 40, 40),
air: (255, 255, 255),
}
plot_xy.basis = 'xy'
plot_xy_cell = openmc.Plot()
plot_xy_cell.filename = 'plot_xy_cells'
plot_xy_cell.width = plot_xy.width
plot_xy_cell.pixels = (1000, 1000)
plot_xy_cell.color_by = 'cell'
plot_xy_cell.basis = 'xy'
plots = openmc.Plots([plot_xy, plot_xy_cell])
plots.export_to_xml()
openmc.plot_geometry()
print("Geometry plots rendered to plot_xy_materials.png / plot_xy_cells.png")
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%%%%%%%%%
################## %%%%%%%%%%%%%%%%%%%%%%%
################### %%%%%%%%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%%%%%%
##################### %%%%%%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%
################# %%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%
############ %%%%%%%%%%%%%%%
######## %%%%%%%%%%%%%%
%%%%%%%%%%%
| The OpenMC Monte Carlo Code
Copyright | 2011-2026 MIT, UChicago Argonne LLC, and contributors
License | https://docs.openmc.org/en/latest/license.html
Version | 0.16.0
Commit Hash | 617d35a5063c57796b43428bc401e627d2011046
Date/Time | 2026-10-10 10:20:21
OpenMP Threads | 6
Reading settings XML file...
Reading materials XML file...
Reading geometry XML file...
Reading tallies XML file...
Preparing distributed cell instances...
Reading plot XML file...
=======================> PLOTTING SUMMARY <========================
Plot ID: 1
Plot file: plot_xy_materials.png
Universe depth: -1
Plot Type: Slice
Origin: 0 0 0
Width: 20.82 20.82
Coloring: Materials
Basis: XY
Pixels: 1000 1000
Plot ID: 2
Plot file: plot_xy_cells.png
Universe depth: -1
Plot Type: Slice
Origin: 0 0 0
Width: 20.82 20.82
Coloring: Cells
Basis: XY
Pixels: 1000 1000
Processing plot 1: plot_xy_materials.png...
Processing plot 2: plot_xy_cells.png...
Geometry plots rendered to plot_xy_materials.png / plot_xy_cells.png
/home/samarth/miniconda3/envs/openmc-env/lib/python3.13/site-packages/openmc/plots.py:1389: FutureWarning: The Plot class is deprecated. Use SlicePlot for 2D slice plots or VoxelPlot for 3D voxel plots. warnings.warn(
7. Run the simulation¶
This launches the actual Monte Carlo transport calculation. Watch the live tally output scroll by — each line is one "batch" of 2000 neutron histories. The last several ("active") batches are averaged to estimate keff.
In [18]:
openmc.run()
sp_path = sorted(glob.glob('statepoint.*.h5'))[-1]
print("\nStatepoint written to:", sp_path)
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%%%%%%%%%
################## %%%%%%%%%%%%%%%%%%%%%%%
################### %%%%%%%%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%%%%%%
##################### %%%%%%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%%
####################### %%%%%%%%%%%%%%%%%
###################### %%%%%%%%%%%%%%%%%
#################### %%%%%%%%%%%%%%%%%
################# %%%%%%%%%%%%%%%%%
############### %%%%%%%%%%%%%%%%
############ %%%%%%%%%%%%%%%
######## %%%%%%%%%%%%%%
%%%%%%%%%%%
| The OpenMC Monte Carlo Code
Copyright | 2011-2026 MIT, UChicago Argonne LLC, and contributors
License | https://docs.openmc.org/en/latest/license.html
Version | 0.16.0
Commit Hash | 617d35a5063c57796b43428bc401e627d2011046
Date/Time | 2026-10-10 10:20:22
OpenMP Threads | 6
Reading settings XML file...
Reading cross sections XML file...
Reading materials XML file...
Reading geometry XML file...
Reading U234 from
/home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_U234.h5
Reading U235 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_U235.h5 Reading U238 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_U238.h5 Reading U236 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_U236.h5
Reading O16 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_O16.h5
Reading O17 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_O17.h5 Reading Zr90 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Zr90.h5 Reading Zr91 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Zr91.h5
Reading Zr92 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Zr92.h5
Reading Zr94 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Zr94.h5
Reading Zr96 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Zr96.h5 Reading H1 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_H1.h5 Reading H2 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_H2.h5 Reading B10 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_B10.h5
Reading B11 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_B11.h5 Reading C0 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_C0.h5 Reading N14 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_N14.h5 Reading N15 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_N15.h5 Reading Ar36 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Ar36.h5 WARNING: Negative value(s) found on probability table for nuclide Ar36 at 294K Reading Ar38 from /home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Ar38.h5
Reading Ar40 from
/home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_Ar40.h5
Reading c_H_in_H2O from
/home/samarth/work/openmc-reactor-circle/xs_data/ENDFB-7.1-NNDC_c_H_in_H2O.h5
Minimum neutron data temperature: 294 K
Maximum neutron data temperature: 294 K
Reading tallies XML file...
Preparing distributed cell instances...
Reading plot XML file...
Writing summary.h5 file...
Maximum neutron transport energy: 20000000 eV for U235
Initializing source particles...
====================> K EIGENVALUE SIMULATION <====================
Bat./Gen. k Average k
========= ======== ====================
1/1 0.12954
2/1 0.14234
3/1 0.12837
4/1 0.12916
5/1 0.16880
6/1 0.13001
7/1 0.14507
8/1 0.15345
9/1 0.14145
10/1 0.14912
11/1 0.14615
12/1 0.14094
13/1 0.12562
14/1 0.15082
15/1 0.14965
16/1 0.15604
17/1 0.14723 0.15164 +/- 0.00441
18/1 0.14983 0.15104 +/- 0.00261
19/1 0.14182 0.14873 +/- 0.00295
20/1 0.14450 0.14788 +/- 0.00244
21/1 0.15098 0.14840 +/- 0.00206
22/1 0.13234 0.14611 +/- 0.00288
23/1 0.14854 0.14641 +/- 0.00251
24/1 0.14428 0.14617 +/- 0.00223
25/1 0.13288 0.14484 +/- 0.00240
26/1 0.13807 0.14423 +/- 0.00225
27/1 0.14296 0.14412 +/- 0.00206
28/1 0.14254 0.14400 +/- 0.00190
29/1 0.13572 0.14341 +/- 0.00185
30/1 0.13480 0.14284 +/- 0.00182
31/1 0.15393 0.14353 +/- 0.00184
32/1 0.15000 0.14391 +/- 0.00177
33/1 0.12838 0.14305 +/- 0.00188
34/1 0.12413 0.14205 +/- 0.00204
35/1 0.12830 0.14136 +/- 0.00205
36/1 0.12545 0.14061 +/- 0.00209
37/1 0.14287 0.14071 +/- 0.00200
38/1 0.16425 0.14173 +/- 0.00217
39/1 0.13390 0.14141 +/- 0.00210
40/1 0.14210 0.14143 +/- 0.00201
41/1 0.12812 0.14092 +/- 0.00200
42/1 0.12823 0.14045 +/- 0.00198
43/1 0.14014 0.14044 +/- 0.00191
44/1 0.13782 0.14035 +/- 0.00184
45/1 0.13153 0.14006 +/- 0.00181
46/1 0.16606 0.14090 +/- 0.00194
47/1 0.14315 0.14097 +/- 0.00188
48/1 0.14016 0.14094 +/- 0.00182
49/1 0.13804 0.14086 +/- 0.00177
50/1 0.13814 0.14078 +/- 0.00172
51/1 0.11028 0.13993 +/- 0.00187
52/1 0.13132 0.13970 +/- 0.00184
53/1 0.16183 0.14028 +/- 0.00188
54/1 0.15313 0.14061 +/- 0.00186
55/1 0.13521 0.14048 +/- 0.00182
56/1 0.14060 0.14048 +/- 0.00177
57/1 0.12346 0.14007 +/- 0.00178
58/1 0.11959 0.13960 +/- 0.00180
59/1 0.14099 0.13963 +/- 0.00176
60/1 0.13616 0.13955 +/- 0.00172
Creating state point statepoint.60.h5...
=======================> TIMING STATISTICS <=======================
Total time for initialization = 9.0494e-01 seconds
Reading cross sections = 8.9567e-01 seconds
Total time in simulation = 9.8339e-01 seconds
Time in transport only = 9.6919e-01 seconds
Time in inactive batches = 2.3441e-01 seconds
Time in active batches = 7.4898e-01 seconds
Time synchronizing fission bank = 5.1255e-03 seconds
Sampling source sites = 4.3489e-03 seconds
SEND/RECV source sites = 7.7066e-04 seconds
Time accumulating tallies = 8.2718e-04 seconds
Time writing statepoints = 3.2119e-03 seconds
Total time for finalization = 1.3319e-02 seconds
Total time elapsed = 1.9086e+00 seconds
Calculation Rate (inactive) = 127981 particles/second
Calculation Rate (active) = 120163 particles/second
============================> RESULTS <============================
k-effective (Collision) = 0.13797 +/- 0.00176
k-effective (Track-length) = 0.13955 +/- 0.00172
k-effective (Absorption) = 0.13856 +/- 0.00192
Combined k-effective = 0.13900 +/- 0.00163
Leakage Fraction = 0.70222 +/- 0.00150
Statepoint written to: statepoint.60.h5
In [19]:
fig, axes = plt.subplots(1, 2, figsize=(13, 6.5))
for ax, fname, title in [
(axes[0], 'plot_xy_materials.png', 'Materials\n(red=fuel, grey=clad, blue=water, black=B4C)'),
(axes[1], 'plot_xy_cells.png', 'Cells (pin-cell structure)'),
]:
img = mpimg.imread(fname)
ax.imshow(img)
ax.set_title(title, fontsize=10)
ax.axis('off')
plt.tight_layout()
plt.show()
8b. Flux and fission (power) maps, plus keff¶
Flux shows where neutrons travel; fission shows where power is actually produced — notice how the control-rod cross suppresses both locally.
In [20]:
with openmc.StatePoint(sp_path) as sp:
tally_result = sp.get_tally(name='flux_and_fission')
flux = tally_result.get_slice(scores=['flux']).mean.reshape(100, 100)
fission = tally_result.get_slice(scores=['fission']).mean.reshape(100, 100)
keff = sp.keff
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
im0 = axes[0].imshow(flux, origin='lower', cmap='viridis')
axes[0].set_title('Neutron flux distribution')
axes[0].axis('off')
plt.colorbar(im0, ax=axes[0], fraction=0.046)
im1 = axes[1].imshow(fission, origin='lower', cmap='hot')
axes[1].set_title('Fission (power) distribution')
axes[1].axis('off')
plt.colorbar(im1, ax=axes[1], fraction=0.046)
plt.tight_layout()
plt.show()
In [21]:
print(f"k-effective = {keff.nominal_value:.5f} +/- {keff.std_dev:.5f}")
if keff.nominal_value > 1.0:
print(" -> SUPERCRITICAL (neutron population growing)")
elif keff.nominal_value < 1.0:
print(" -> SUBCRITICAL (neutron population shrinking)")
else:
print(" -> CRITICAL (steady chain reaction)")
k-effective = 0.13900 +/- 0.00163 -> SUBCRITICAL (neutron population shrinking)
9. 🧪 Now experiment!¶
Pick ONE change below, re-run from the relevant section downward (noted per experiment) — not from the very top, that just redefines things that don't need to change — and see what happens to keff and the power map.
Discuss with your group before changing: do you think keff will go up or down, and why?
| # | Change | Where | Re-run from |
|---|---|---|---|
| a | Pull a control rod: change one C to F in the pattern grid |
§3c | Section 3 |
| b | Change enrichment: enrichment=4.0 → 2.0 or 5.0 |
§2a | Section 2 |
| c | Remove boron: delete water.add_element('B', 500e-6) |
§2c | Section 2 |
| d | Bigger core: N = 7 → 9 or 11 (extend the pattern grid too!) |
§3c | Section 3 |
| e | Water density: set_density('g/cm3', 1.0) instead of 0.74 |
§2c | Section 2 |
After editing the target cell, just re-run every cell from that section through Section 8 again (Cell menu → "Run All Below" works well here).
10. Wrap-up discussion¶
- Why do real reactors need active, continuous control — not just "reach critical and done"?
- What's the difference between what we modeled (a static k-eigenvalue snapshot) and a real reactor transient — time-dependent behavior, fuel temperature feedback, xenon poisoning building up over hours? OpenMC can do depletion/burnup calculations too — that's a "part 2" topic.
- This model is a teaching toy: the pin dimensions and general pattern are realistic PWR-ish numbers, but it is not a licensed, validated reactor design. Real reactor physics calculations require validated nuclear data, full 3D geometry, depletion, thermal-hydraulic feedback, and rigorous QA — well beyond what fits in an hour.
Troubleshooting¶
ModuleNotFoundError: No module named 'openmc'→ Jupyter wasn't started bysetup.sh. Close it and run./setup.shagain.- "Nuclear data MISSING" in Section 0 → the notebook was opened from
outside the project folder. Open it from the folder that contains
xs_data/. - Geometry plot is all one color → check the
colorsdict in Section 6 matches the materials you actually defined.