Skip to content

Repository files navigation

MicroHTGR

A modular OpenMC framework for full-core neutronics, criticality search, spatially resolved depletion, and coupled neutronics/thermal-hydraulics analysis of prismatic high-temperature gas-cooled reactors.

License: MIT OpenMC Python Attribution: CC BY 4.0

Full-core radial cross-section of a prismatic HTGR generated by this framework: hexagonal fuel assemblies with TRISO-loaded compacts and helium coolant channels, seven B4C control rod positions, six reserve shutdown tubes, and a BeO radial reflector.

A full-core prismatic HTGR built entirely from a parameter file, plotted directly by the framework.


Contents


What this is

MicroHTGR is a research framework for analysing prismatic (block-type) high-temperature gas-cooled reactors with OpenMC. It takes a reactor described as a dictionary of parameters, builds a full Monte Carlo transport model from it, runs one of seven study types against that model, and post-processes the results into figures and CSV files.

It was written to answer the questions a fuel cycle and safety analysis actually needs to answer:

  • Where do the control rods have to sit for the core to be exactly critical, at every point in life?
  • What is the burnup distribution across the core, not just its average?
  • What are the fuel, moderator and isothermal temperature coefficients, and do they stay negative?
  • What temperatures does the fuel actually reach, once the power shape and the coolant temperature are solved consistently rather than assumed?
  • How much fast fluence does the reflector accumulate before it has to be replaced?
  • What is the shutdown margin with one rod stuck out?

It was developed for a 10 MWth microreactor design and run end to end on it (see Worked example), but nothing in the framework is specific to that reactor.

The core idea: the reactor lives in a config file

Reactor physics scripts tend to be written around one core, so changing the layout means going back into the geometry code and rewriting it.

Here the geometry is generated from a description, in which a core is written as a list of concentric rings of assemblies with each position named by a short code:

"core_rings": [
    ["rr", "f", "f"] * 6,    # Outer ring: reflector, fuel, fuel, repeated 6x
    ["f", "fc2"] * 6,        # Next ring in: fuel alternating with bank-2 control assemblies
    ["fss"] * 6,             # Fuel assemblies with reserve shutdown tubes
    ["fcp1"],                # Centre: fuel + bank-1 control rod + burnable poison
],

That layout, plus the parameters around it, is the entire core definition. The framework expands it into hexagonal lattices, axially zoned universes, TRISO particle lattices, control rod channels with continuously positioned absorber, poison rods, reflector regions and tally meshes.

The available assembly codes:

Code Assembly
f Fuelled, no rods
fp / fpa Fuelled with 6 corner burnable poison rods / 1 central poison rod
fc1 fc2 fc3 Fuelled with a central control rod on bank 1, 2 or 3
fcp1 fcp2 fcp3 Fuelled with a central banked control rod and 6 corner poison rods
fss / fssp Fuelled with a central reserve shutdown tube (optionally + poison rods)
rr Reflector block
r1 r2 r3 Reflector block with a central banked control rod
ra1 ra2 ra3 Reflector block with 3 control rods in a hexagonal ring (1/6 geometry only)
rss / rssa Reflector block with 1 / 3 reserve shutdown tubes

Everything else is a key in the same dictionary: enrichment, TRISO layer thicknesses, packing fraction, compact and coolant channel dimensions, lattice and bundle pitch, core height, axial zoning, reflector material and thickness, absorber enrichment, rod bank assignments, depletion schedule, tally configuration and Monte Carlo settings.

The practical consequence is that a different prismatic HTGR core is a config edit rather than a code change, whether it is graphite- or BeO-reflected, single- or multi-bank, poisoned or clean, and at any ring count the lattice will hold. Only CHUDR-like cores have actually been run, however, so treat larger Fort St. Vrain-style or INL HTGTR-derived layouts as supported by the code rather than as exercised. New assembly types can be added in assembly.py alongside the existing ones.

Capabilities

Transport Full-core or 1/6-symmetric hexagonal geometry, explicit TRISO particle lattices, packed-sphere absorber lattices
Double heterogeneity Reactivity-equivalent Physical Transform (RPT) with an automated calibration mode that solves for the transform radius against an explicit-TRISO reference
Criticality search Binary search on control rod insertion to a user-set k_eff tolerance, with continuous (not zone-snapped) rod positions
Depletion Spatially resolved burnup on a radial × axial zone map, reduced-chain generation, graphite depletion, restart from a partial run
Rod-following depletion Criticality search re-run at every depletion timestep, so the core is depleted at its true critical rod position throughout life
Multiphysics Converged neutronics ↔ single-channel thermal-hydraulics iteration on both k_eff and the axial heating profile
Reactivity coefficients FTC, MTC and ITC by direct perturbation, at beginning, middle or end of life
MOL/EOL analysis Re-analyse any depletion step (coefficients, heat maps, leakage spectra) without re-running the depletion
Post-processing Ten modules producing burnup, isotopics, peaking, spectra, reflector fluence, power cycle output and parametric summaries

Installation

OpenMC is distributed through conda-forge, so conda (or mamba) is the supported route:

git clone https://github.com/c-finney/MicroHTGR.git
cd MicroHTGR
conda env create -f environment.yml
conda activate microhtgr

If you already have OpenMC installed by another means, requirements.txt covers the remaining dependencies.

Cross-section data

The framework needs a continuous-energy cross-section library and a depletion chain, neither of which is bundled here. Download them from the OpenMC data library page, then point the framework at them. The reference results below were produced with ENDF/B-VIII.0:

export OPENMC_CROSS_SECTIONS="/path/to/endfb-viii.0-hdf5/cross_sections.xml"
export OPENMC_DEPLETION_CHAIN="/path/to/chain_endfb81_thermal.xml"   # optional
export MICROHTGR_OUTPUT_DIR="/path/to/output"                        # optional

Only OPENMC_CROSS_SECTIONS is required. The depletion chain defaults to a file sitting alongside the cross-section library, and output defaults to ../MicroHTGR_Output beside the repository, so nothing inside the repository needs editing to run it somewhere else.

Thermal scattering: the graphite materials use the c_Graphite S(α,β) treatment. Your library must include it, or thermal HTGR results will be badly wrong.

Quick start

Run a single steady-state eigenvalue calculation on the shipped core:

# In config.py, set:  "study_execution_mode": "SingleStudy"
python main_simulation.py

Output lands in a timestamped directory, htgr_run_MM.DD.YYYY_HH.MM.SS/, containing the OpenMC inputs, statepoints, geometry plots, a run_params.json snapshot of every parameter used, and a post-processed results folder.

To reproduce the full reference fuel cycle analysis:

cp examples/chudr_cs_depletion_config.py config.py    # back up your own config.py first
python main_simulation.py

That example is documented inline and is the configuration behind every figure below. It runs 24 depletion steps at 50,000 particles by 200 batches, each one carrying its own criticality search and thermal-hydraulic coupling loop, so expect days of wall time on a many-core machine and tens of gigabytes of output. The cost note at the top of the file explains how to scope it down for a first run.

Study modes

Set study_execution_mode in config.py:

Mode What it does
SingleStudy One steady-state eigenvalue calculation of the configured core.
ParametricStudy Sweeps any single parameter across a list of values, runs each case, and aggregates them into comparison plots. Used for pitch optimisation, reflector thickness scans, rod worth curves.
ReactivityStudy FTC, MTC and ITC by direct perturbation across a set of ΔT values.
CriticalSearch Inserts bank 1 fully, then binary-searches bank 2 until the core is critical.
DepletionStudy All-rods-out depletion over the configured timestep list.
CSDepletionStudy Depletion with a criticality search at every timestep (rods follow reactivity down as the fuel burns), optionally with thermal-hydraulic coupling at each step.
RPTCalibration Runs an explicit-TRISO reference, then Illinois regula falsi on the RPT radius until the homogenised model reproduces the explicit k_eff.

mol_eol_analysis.py runs separately against a completed depletion directory, reading depleted compositions straight out of depletion_results.h5 and injecting them into a fresh model build. That means middle- and end-of-life reactivity coefficients, heat maps and leakage spectra cost one eigenvalue solve each, not a re-run of the fuel cycle.

Building a different reactor

A worked change, converting a BeO-reflected poisoned core into a graphite-reflected core with a different rod pattern and 15.5% enrichment:

params = {
    "enrichment": 0.155,               # was 0.1975
    "triso_pf": 0.40,                  # was 0.33
    "use_BeO_reflector": False,        # graphite radial reflector instead
    "core_rings": [
        ["rr"] * 18,                   # solid reflector ring
        ["f", "f", "r1"] * 4,          # fuel with bank-1 rods in the reflector blocks
        ["f"] * 6,
        ["fc2"],                       # central bank-2 rod
    ],
    "n_ax_zones": 30,                  # was 50
    "use_1/6_geometry": False,         # model the whole core
    ...
}

No other file changes. Re-run RPTCalibration afterwards if you use the homogenised fuel treatment, since the transform radius depends on compact geometry, packing fraction and enrichment. Additionally, check that ax_zones_per_burnup_region still divides the new n_ax_zones exactly.

Things worth knowing when you do this:

  • ax_zones_per_burnup_region must divide n_ax_zones exactly.
  • ra* and rssa* assembly codes only work in 1/6 geometry.
  • With use_homogenized_fuel = True and an uncalibrated rpt_radius, results will be wrong in a way that looks plausible. Calibrate, or set it to False and model particles explicitly.
  • BeO_thickness is clamped so the reflector cannot extend past core_radius.

Coupled neutronics / thermal-hydraulics

Assuming a cosine power shape and a linear coolant temperature rise is the usual simplification, and it breaks down for a rodded core, where the axial power profile is strongly distorted and the fuel temperature feedback depends on that shape.

This framework closes the loop instead, with th_coupler() in mol_eol_analysis.py iterating five steps until the temperature field and the power shape agree:

  1. Solve the eigenvalue problem at the current temperature field.
  2. Extract the axial heating profile from the mesh tally, for the hottest and average channels.
  3. Feed that profile to the single-channel solver (ThermalHydraulics/nc_htgr.py), which solves the coolant enthalpy rise, convective film drop, graphite conduction and compact/kernel temperatures node by node, and closes a recuperated Brayton cycle.
  4. Map the resulting coolant, compact and matrix temperatures back onto the OpenMC axial zones.
  5. Repeat until both k_eff and the heating profile converge.

In CSDepletionStudy this runs inside the criticality search at every depletion timestep, so the converged temperatures are thus carried into the burn step rather than assumed.

Four temperature sources are selectable through temp_profile_source, so the expensive coupling is opt-in:

Source Use for
THCoupling Converged multiphysics. The accurate, expensive option.
IdealProfile Cosine/linear analytic profile from config.py min/max values. Good for scoping.
Isothermal One temperature everywhere. Benchmarking and code-to-code comparison.
FromCSV Reuse a converged profile from a previous coupled run. Near-free, and the right choice for parameter sweeps where the temperature field barely moves.

Post-processing

Runs automatically unless run_post_processing is disabled; each module also works standalone against a finished run directory.

Module Produces
burnup_estimation Cycle length, k_eff, leakage fraction, uranium mass and fuel volume (analytic, no stochastic volume calculation)
depletion_postprocessing k_eff and rod position histories, nuclide inventories by group, discharge burnup, conversion ratio, fissile inventory ratio, B-10 burnout, per-step net electrical output
heating_profile_extraction Hottest- and average-channel axial heating profiles, integrated heating maps, peaking factors
tally_plotter Flux, fission rate and heating distributions in XY / XZ / YZ
spectrum_thermalization Energy spectrum per unit lethargy, thermal/epithermal/fast fractions, thermalization metrics
leakage_spectrum Leakage spectra and absolute rates across every core boundary
BeO_depletion_postprocessing Reflector peak and area-averaged fast flux and cumulative fluence vs. burnup
reactivity_coefficients_postprocessing Coefficient plots and tabulated results
rod_worth_postprocessing Integral and differential control rod worth curves from an insertion sweep, with propagated uncertainties
parametric_postprocessing Cross-case comparison plots and summary tables for a parameter sweep

Every run directory also carries run_params.json, a complete snapshot of the parameters that produced it, which is what mol_eol_analysis.py reads back to reconstruct a model.

Worked example: the CHUDR microreactor

The framework was developed for CHUDR (Compact Helium-cooled Universal Defense Reactor), a 4 MWe / 10 MWth transportable prismatic HTGR designed as a senior design project at the University of Florida and entered in the ANS Student Design Competition. Everything below is output of this repository, from the configuration in examples/chudr_cs_depletion_config.py.

CHUDR is 31 fuel assemblies, 7 control rod positions, 6 reserve shutdown tubes and 6 burnable poison rods inside a 20 cm BeO radial reflector of 1.8 m active diameter and 2.379 m active height, fuelled with 19.75% HALEU UCO TRISO at a packing fraction of 0.33 across 1,218 fuel and 558 coolant channels.

Geometry

Radial material plot of the 1/6-symmetric core sector, showing hexagonal fuel assemblies with red helium coolant channels and green fuel compacts in a blue graphite matrix, bounded by the light blue BeO radial reflector. Axial material plot through the core, showing the active fuel region between upper and lower graphite reflectors.

Left: radial material plot of the 1/6 sector actually simulated. Right: axial section. Both plotted by the framework at build time.

Rod-following depletion

k-effective versus time over 3.19 years. The beginning-of-step critical search values sit on k equals 1.000 throughout the cycle, while end-of-step values fall below it, and the all-rods-out curve crosses 1.0 at discharge. Control rod bank insertion fraction versus time, showing progressive withdrawal over the fuel cycle as fuel depletes.

The orange beginning-of-step curve sits on k = 1.000 across the full 3.19 years because the criticality search re-solves the rod position at every timestep, while the blue end-of-step curve is the reactivity actually consumed by each burn interval. Cycle length is set by where the all-rods-out core falls to critical.

Fuel cycle and isotopics

Core-average, peak and minimum channel burnup in megawatt-days per metric ton of uranium versus time. Fissile actinide inventories for uranium-235, plutonium-239 and plutonium-241 versus burnup.

Spatially resolved depletion gives peak and minimum channel burnup rather than only the core average, which is the number that has to be checked against the 160,000 MWd/MtU TRISO limit.

Coupled thermal-hydraulics

Peak fuel temperature versus time over the fuel cycle, dipping at middle of life as rod withdrawal flattens the axial power profile. Net electrical output versus time, varying by under two percent across the 3.19 year cycle.

Because the temperature field is converged against the true rodded power shape at every step, peak fuel temperature comes out non-monotonic, dipping at middle of life as rod withdrawal flattens the axial profile and recovering later as radial peaking grows.

Reference results

Quantity Value
Cycle length (once-through) 3.19 years
Average discharge burnup 80,098 MWd/MtU
Peak / minimum channel burnup 91,276 / 38,639 MWd/MtU
Isothermal temperature coefficient −7.81 pcm/K (BOL) → −6.35 pcm/K (EOL)
Shutdown margin, all rods in 13,334 ± 19 pcm
Shutdown margin, one stuck rod 8,286 ± 18 pcm
Shutdown margin, reserve system 5,928 ± 17 pcm
Net electrical output 4.15 – 4.24 MWe across the cycle
BeO reflector lifetime fast fluence 5.07 × 10²⁰ n/cm²

The RPT homogenisation used for the depletion runs was validated against explicit-TRISO modelling: discharge burnup agreed to within 2.5%, with isotopic errors below 1% for every actinide and fission product of interest.

Limitations

Things worth knowing before trusting a result out of this framework.

There is no test suite and no CI. Confidence in the physics rests on the RPT homogenisation agreeing with explicit-TRISO modelling and on the reference results reproducing the values reported for CHUDR, not on a regression suite. A change to the geometry or tally code will not be caught automatically.

Linux only, in practice. Development and every result quoted here used Linux with Python 3.13.12 and OpenMC 0.15.3. Nothing in the code is deliberately platform-specific, but Windows and macOS have not been exercised.

The thermal-hydraulics model is a single channel. nc_htgr.py solves one representative channel at a time, chosen as average, hot or cold, so cross-flow between channels, inter-assembly mixing and bypass flow all sit outside it. That is adequate for the axial temperature feedback the coupling needs and for scoping, however, and it is thus not a substitute for a subchannel or CFD calculation.

Depletion runs at fixed thermal power. The power level is a single constant taken from thermal_power_MW and held across the whole cycle, so load-following and power history effects are not modelled.

The reference results come from a 1/6-symmetric model at 50,000 particles per batch over 200 batches. The quoted uncertainties belong to that configuration, and a full-core run at different statistics will not reproduce them exactly.

Only CHUDR-scale cores have been run end to end. See the note under the core idea.

Repository layout

MicroHTGR/
├── main_simulation.py          Entry point. Model build + all seven study modes.
├── config.py                   THE reactor definition. Path config, all parameters.
├── materials.py                OpenMC material definitions, built from config.
├── assembly.py                 Assembly universes and core lattice construction.
├── trisos.py                   TRISO particle packing and lattice construction.
├── b4c_spheres.py              Reserve shutdown B4C sphere packing and lattice.
├── reactivity_coefficients.py  FTC / MTC / ITC by direct perturbation.
├── mol_eol_analysis.py         MOL/EOL re-analysis + the N-TH coupler.
│
├── ThermalHydraulics/          Single-channel TH + Brayton cycle solver
│   ├── nc_htgr.py                (see NOTICE.md, originates from HTGR-SCAPC)
│   ├── nc_input.csv              TH and power cycle input deck
│   └── neutronics.csv            Example axial heating profile
│
├── PostProcessingScripts/      Ten analysis and plotting modules
│
├── examples/
│   └── chudr_cs_depletion_config.py    Reference coupled CS-depletion run
│
├── docs/figures/               Figures used in this README
├── licenses/CC-BY-4.0.txt      Licence text for the derived VTB material
├── LICENSE                     MIT
├── NOTICE.md                   Required third-party attribution
└── CITATION.cff                Citation metadata

About 16,100 lines of Python across 21 modules, with simulation output written outside the repository under the directory named by MICROHTGR_OUTPUT_DIR.

Attribution and license

This repository is released under the MIT License.

It contains material derived from work by others, and redistribution requires carrying that attribution forward. NOTICE.md is the authoritative statement, and the summary here is not a substitute for it.

  • NRIC Virtual Test Bed, prismatic HTGR assembly model. htgr/assembly, documentation. Licensed CC BY 4.0. The material definitions in materials.py, the prismatic assembly parameterisation in config.py, and the assembly construction approach in assembly.py are derived from it. This repository is a modified derivative and is not the VTB model; NOTICE.md lists the changes. Neither INL, NRIC, nor the DOE endorses or has reviewed this work.

  • HTGR-SCAPC by Bryan Huynh. github.com/bryanhhuynh/HTGR-SCAPC. ThermalHydraulics/nc_htgr.py originates from this project and is included with the modifications listed in NOTICE.md.

  • OpenMC. MIT-licensed, developed by MIT and contributors. A dependency, not a derivative.

Citing this work

If this framework is useful in your research, please cite it via CITATION.cff (GitHub's "Cite this repository" button will format it for you), and cite the VTB model and OpenMC alongside it, since both are listed as references in that file.

Acknowledgements

The CHUDR design that this framework was built to analyse was produced by a six-person senior design team in the ENU 4192 course at the University of Florida: Cade Finney, Evan Alder, Bryan Huynh, Colin Frazier, William St. Peter and Daniel Fernandez, under the guidance of Dr. DuWayne Schubring and the faculty of the Department of Materials Science & Engineering, Nuclear Engineering Program.

Thanks to the OpenMC development team, and to Idaho National Laboratory and NRIC for publishing the Virtual Test Bed openly, since having a credible open prismatic HTGR reference model to start from is why this project began where it did.

About

Modular OpenMC neutronics, depletion, and N-TH coupling framework for prismatic HTGRs

Topics

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages