Multiphysics Coupling

Every OpenMC model up to this point has been isothermal, and every one of them has been a fiction. A reactor at power has fuel near 900 K at the pellet surface and well above 1200 K on the centerline, cladding a few hundred degrees cooler, and coolant that heats and expands as it rises through the core. Each of those affects the transport calculation, and the transport calculation is what set the power that caused them.

Fuel temperature acts through Doppler broadening. The absorption resonances of U-238 are narrow features in energy, and thermal motion of the target nuclei smears each one out: the peak falls and the wings rise, conserving the area under the resonance. Because a tall narrow resonance is self-shielded — the flux is depressed at exactly the energies where the cross section is largest, so the interior of the fuel barely sees the peak — spreading the resonance out increases the number of absorptions even though the integral of the cross section did not change. Hotter fuel therefore captures more neutrons, and k falls. This is the fuel temperature coefficient, it is negative, and it is the fastest-acting mechanism that makes a power reactor stable.

Coolant temperature acts through density. Water expands as it heats, so a hotter channel holds fewer hydrogen atoms per cubic centimeter, moderates less effectively, and shifts the spectrum harder. In a pressurized water reactor, which is designed to be under-moderated for exactly this reason, less moderation means fewer thermal fissions and k falls again. The under-moderated design is what makes the sign negative — an over-moderated lattice would gain reactivity as the coolant heated, which is why the ratio of moderator to fuel is a safety parameter rather than a convenience.

The two together close a loop. Power sets temperatures, temperatures set cross sections, cross sections set the power distribution. Neither half can be solved without the other, and coupling is the business of finding the state where both are simultaneously satisfied.

OpenMC accepts a temperature in three places, and the widely repeated claim that temperature belongs on cells rather than materials is not right as stated. All three are valid, and they form a precedence chain: a temperature on a cell wins over one on the material filling it, which wins over the global default in the settings.

import openmc

fuel = openmc.Material(name='UO2', temperature=900.0)   # valid: material-level
fuel.set_density('g/cm3', 10.4)
fuel.add_nuclide('U235', 0.045)
fuel.add_nuclide('U238', 0.955)
fuel.add_element('O', 2.0)

fuel_or = openmc.ZCylinder(r=0.4096)
fuel_cell = openmc.Cell(fill=fuel, region=-fuel_or)
fuel_cell.temperature = 1200.0    # wins over the material's 900 K

settings = openmc.Settings()
settings.temperature = {'default': 293.6}   # used only where neither is given

What is true, and what the folk rule is reaching for, is narrower: a temperature that varies between cells sharing one material has to live on those cells, because the material object can only hold one value. A pin cell whose fuel is at a single temperature can carry it on the material perfectly well — OpenMC's own openmc.model.borated_water() helper passes temperature= straight to the material constructor.

In a coupled calculation, though, the practical answer is cells, since that is the resolution at which the thermal solver will hand you a field. A radial temperature profile through a pellet needs the pellet divided into annular cells, and the number of rings you choose is the resolution of your feedback.

Continuous-energy libraries are distributed at a handful of temperatures — commonly 294, 600, 900, 1200, and 2500 K. A cell at 1147 K falls between two of them, and the method entry decides what happens.

settings.temperature = {
    'default': 293.6,
    'method': 'interpolation',   # 'nearest' is the default
    'tolerance': 10.0,           # how far 'nearest' may reach, in K
    'range': (294.0, 1500.0),    # preload this span, for runtime temperature changes
    'multipole': True,           # use windowed multipole data where available
}

With nearest, OpenMC snaps to the closest library temperature and raises an error if the nearest one is further away than tolerance. That error is a feature — it stops a model from quietly running 1147 K fuel on 900 K data. With interpolation, OpenMC loads the bracketing temperatures and interpolates stochastically between them, at the cost of holding both in memory.

Windowed multipole data is the more elegant option where it exists. Rather than storing pre-broadened tables, it stores the resonance parameters in a form that can be Doppler broadened analytically at whatever temperature is asked for, which makes temperature a continuous variable and removes the memory cost of holding several tables. It requires a WMP library alongside the usual HDF5 data, and it covers the resolved resonance range, so it supplements the standard tables rather than replacing them. The range entry matters when temperatures will change during a run: it tells OpenMC to load data spanning that interval up front.

One thing temperature does not do automatically is change density. Doppler broadening and thermal expansion are separate physics, and OpenMC will not thin your coolant because you told it the coolant is hot. If the density should change, you have to change it — materials covers set_density.

Almost all production coupling is operator splitting — a Picard iteration in which each physics is solved to convergence in turn with the other held fixed, repeating until neither moves. Transport gives a power distribution, the thermal solver turns that into a temperature and density field, those go back into the transport model, and you stop when the temperatures stop changing.

The step people omit is the one where the change actually reaches the solver. Setting cell.temperature on a Python object does nothing to a run unless the model is re-exported, because OpenMC's transport solver reads XML. A loop that mutates the geometry and calls run() again will iterate cheerfully and converge to the answer it started with.

import numpy as np
import openmc

model = build_model()      # geometry, materials, settings, tallies

for iteration in range(max_iterations):
    # 1. Transport. model.run() returns the path to the statepoint, not a
    #    StatePoint object.
    sp_path = model.run()

    # 2. Recover the volumetric heat source. kappa-fission scores the
    #    recoverable energy released per fission, in eV per source particle.
    with openmc.StatePoint(sp_path) as sp:
        tally = sp.get_tally(name='power')
        heat = tally.get_reshaped_data(value='mean').squeeze()

    # 3. Thermal solve, in whatever code owns it
    temps = solve_heat_transfer(heat)          # K, one per fuel cell
    densities = coolant_density(temps)         # g/cm3

    # 4. Push the new state back into the model
    for cell_id, T in temps.items():
        model.geometry.get_all_cells()[cell_id].temperature = T
    coolant.set_density('g/cm3', densities.mean())

    # 5. The step that is easy to forget: without this, run() re-reads the
    #    XML from the previous iteration and nothing has changed.
    model.export_to_model_xml()

    # 6. Convergence on the field, not on k — k is noisy at this precision
    if np.allclose(list(temps.values()), previous, rtol=1e-3):
        break
    previous = list(temps.values())

For tight coupling, where re-exporting XML and restarting the solver every iteration is too expensive, OpenMC exposes its internals through openmc.lib. That interface keeps the simulation resident in memory and lets you change a cell temperature between batches without touching a file, which is how the coupled frameworks below drive it.

import openmc.lib

with openmc.lib.run_in_memory():
    for iteration in range(max_iterations):
        openmc.lib.reset()
        openmc.lib.run()

        heat = openmc.lib.tallies[1].mean
        temps = solve_heat_transfer(heat)

        # In memory: no XML, no restart, and this takes effect on the next batch
        for cell_id, T in temps.items():
            openmc.lib.cells[cell_id].set_temperature(T)

The coupling needs a heat source, and which score you want depends on how much of the energy you intend to account for. kappa-fission gives the recoverable energy released per fission, which is the standard choice and is what most thermal solvers expect. It attributes that energy to the fission site, which is a good approximation for the roughly 85 % carried by fission fragments and a poor one for the gamma rays.

mesh = openmc.RegularMesh()
mesh.dimension = [1, 1, 20]                 # 20 axial levels through the pin
mesh.lower_left = (-0.63, -0.63, 0.0)
mesh.upper_right = (0.63, 0.63, 366.0)

power = openmc.Tally(name='power')
power.filters = [openmc.MeshFilter(mesh)]
power.scores = ['kappa-fission']            # eV per source particle

If the gamma transport matters — and it does for a fine axial power shape, since a fraction of the fission energy is deposited centimeters away from where it was released — run coupled neutron–photon transport and score heating instead, which tallies actual energy deposition wherever it happens. It costs more, both in run time and in requiring photon data, and it moves the peak.

A tally in either case is normalized per source particle, so converting to watts needs the reactor power: multiply by the total power and divide by the total energy scored across the whole model. Getting this normalization wrong is the most common reason a first coupled calculation produces temperatures that are physically absurd rather than merely inaccurate.

The loop above is worth writing once, for understanding. Production work uses code that already handles the mesh mapping, the parallel data movement, and the restart logic, since none of that is where the physics is.

Cardinal, from Argonne, wraps OpenMC and the NekRS computational fluid dynamics code as applications inside the MOOSE framework, which supplies the finite-element heat conduction and the transfer machinery between meshes. Because it drives OpenMC through openmc.lib rather than through files, it can iterate without restarting the transport solve. ENRICO takes a narrower approach — a coupling driver built specifically for neutronics with either OpenMC or Shift, against Nek5000 or NekRS or a simpler heat surrogate — and is correspondingly easier to stand up.

For a first coupled calculation, neither is the right starting point. A single pin with twenty axial cells, a one-dimensional conduction model in NumPy, and the Picard loop above will teach you where the convergence problems are and cost an afternoon. Reach for a framework when the mesh mapping between physics becomes the hard part, which is roughly when you move from a pin to an assembly.