Neutronics–Thermal-Hydraulics Coupling: An SMR Example

Neutronics–thermal-hydraulics coupling connects the neutron calculation to the heat-and-flow calculation. Neutronics supplies a heat source; thermal-hydraulics returns temperatures and coolant density that change the neutron calculation.

Imagine a small water-cooled reactor operating steadily. A small withdrawal of a control rod adds positive reactivity: it shifts the neutron balance toward a growing population. Fission power begins to rise. The fuel heats up. That temperature change affects the neutron balance again.

Now give the neutron calculation to one program and the heat-and-flow calculation to another. Which temperature should the first program use for the power it is calculating: the temperature before the power rise, or the temperature that rise is about to produce?

That is the coupling problem. It sits inside simulations of large reactors and small modular reactors alike. I worked on methods for this exchange with Enzo Dari; our 2017 ENIEF contribution on neutronic–thermal-hydraulic coupling is one starting point for this explanation. Here, we will build the intuition with a much smaller model you can manipulate.

What is neutronics–thermal-hydraulics coupling?

The neutronics calculation describes the neutron population and, in a spatial model, where fissions occur. It supplies a power distribution: the heat source for the other calculation.

The thermal-hydraulics calculation describes how that heat moves through the fuel and into the coolant, and how the coolant carries it away. It returns temperatures and density. These change the cross sections used in the neutron calculation, closing the loop.

Two return paths are especially useful to picture in a water-cooled reactor. As fuel temperature rises, Doppler broadening changes the resonance absorption of neutrons. For the operating state assumed here, its net reactivity contribution is negative. As water heats and its density changes, moderation and absorption change too. The sign and strength of that second contribution depend on the design and operating state; we choose a negative temperature coefficient for this example.

Negative feedback means the temperature response opposes the initial power increase. It does not mean that every numerical way of exchanging those temperatures will behave well.

Which variables do the two solvers exchange?

Direction	Field	Meaning
Neutronics → TH	q‴, W/m³	Volumetric heat source
TH → neutronics	Tf and Tc, K	Fuel and coolant temperatures
TH → neutronics	ρc, kg/m³	Coolant mass density

Neutron flux φ is not coolant mass flow. Coolant density ρc is not reactivity ρ; the symbols refer to different quantities. The simplified experiment below represents density feedback through a coolant-temperature coefficient.

Put it in an SMR

Think of an integral pressurised water SMR: the core and major primary-system components share a reactor vessel. CAREM is one real example of the integral, natural-circulation approach. Other water-cooled SMRs use forced circulation. The IAEA SMR catalogue describes these different designs.

In a detailed natural-circulation calculation, heating also changes buoyancy, flow and heat removal. A full model must let those quantities participate in the exchange. That makes an SMR a useful setting for discussing coupling, but small size alone does not tell us whether a numerical method will converge.

For the experiment below, keep just the core's thermal feedback: one average fuel temperature, one average coolant temperature, and a neutron population with one effective delayed-neutron group. Hold coolant flow and inlet temperature fixed. Choose a nominal 100 MW thermal operating point and a small, sustained +50 pcm reactivity step; one pcm is $10^{-5}$ in dimensionless reactivity.

These are illustrative parameters. This is a teaching model of a feedback loop relevant to water-cooled SMRs, rather than a calculation of CAREM or of natural circulation. It has no spatial power shape, pressure dynamics, boiling or plant control system. Those details belong in a plant model; leaving them out here lets us watch the exchange itself.

Why can a single exchange miss the feedback?

A transient simulation divides time into intervals. Suppose the current state is known at $t_n$, and we want the state at $t_{n+1}=t_n+\Delta t$.

The simplest partitioned approach is a single sweep:

1. Give the neutron solver the old temperatures and calculate the new power.
2. Give the thermal solver that new power and calculate the new temperatures.
3. Accept the result and advance time.

It is easy to connect two existing codes this way. Each accepts inputs and returns outputs. With a sufficiently small time step, the exchanged fields may change little enough that the approximation meets the accuracy required.

But the two calculations have not agreed on the new state. The heat calculation may say the fuel has warmed by several kelvin, while the neutron calculation used the fuel temperature from before that warming. For our negative-feedback example, that lag can overpredict the first power increase.

Both individual solvers in the experiment use backward Euler, an implicit time integrator. Their exchange is still lagged. An implicit solver inside each program does not automatically make the coupling between them implicit.

02 / Same interval, different agreement

ONE PASSUse the old temperatures

State at tₙ → N → TH → accept tₙ₊₁

The heat calculation gets new power. Neutronics still used temperatures from tₙ.

CONVERGED COUPLINGAgree on the new temperatures

State at tₙ → [ N ⇄ TH ] → accept tₙ₊₁

Repeat or solve the coupled equations. Every trial belongs to the same tₙ₊₁.

Can Picard iteration solve an implicit coupling?

Yes. Picard can converge an implicit coupling by repeatedly exchanging fields at the same new time level. Whether it converges depends on the coupling map and the numerical settings.

An implicit coupling uses temperatures and power that belong to the same new time level. Neither field is known initially, so we solve for them together.

One accessible way is Picard iteration. Guess a power, calculate temperatures, use those temperatures to calculate another power, and repeat. Keep the old physical state fixed throughout these trials. Only advance time after the exchange meets a specified tolerance.

This distinction matters: the iteration counter is not a clock. Ten coupling iterations may all be trying to determine the state at the same instant, say four seconds after the perturbation.

Picard can solve a fully implicit coupling. It is the default fixed-point method in MOOSE's multiphysics coupling framework, whose documentation explicitly separates iterations within a time step from lagging data across steps.

Sometimes the guesses settle quickly. Sometimes a trial overcorrects the previous one: higher guessed power produces hotter fuel, which returns lower power, which produces cooler fuel, which returns higher power. If successive corrections grow, the iteration fails even though the physical feedback is negative.

Try the exchange

Press Watch the comparison. Picard and Newton start with the same guess: watch which one brings the disagreement down to the dashed line. Try Small · 0.25 s, then Large · 4 s. The physics stays the same; only the numerical time step changes.

Switch to See the power to compare a single exchange with a converged implicit calculation over 60 seconds. A finer-step solution of the same model provides a reference. You can return to Compare the methods at any time.

Try it yourselfCan the two solvers agree?

Same reactor. Same time step. Two ways to find the answer.

01 · Compare the methods02 · See the power

Time step

4.00 s

Small · 0.25 sLarge · 4 s

▶ Watch the comparisonOne step →Time stays at 4.00 s

● Picard repeat the exchange● Newton use the sensitivity

Picard23 updatesIteration failed

The numerical guesses do not settle.

Newton3 updatesAgreed ✓

Power: 105.02% of initial.

Lower is better. Each update is another guess at the same instant—not another second.

What the graph measures

Both methods solve the same implicit first step. The graph shows the actual residual |p − N(TH(p))| on a logarithmic scale; agreement means ≤ 10⁻¹⁰, a teaching tolerance. Values below 10⁻¹² sit at the plot floor. The horizontal range adapts to the longer run. The Picard slope at the solution is -1.17. A failed iteration is not a reactor instability. Newton uses an analytic derivative here, so update counts are not production-code cost comparisons.

The two views answer different questions. Compare the methods asks whether trial fields agree at one particular time. See the power asks what trajectory the time-stepping procedure produces. A failed iteration in the comparison is a numerical convergence problem, not evidence of a reactor oscillation.

Also notice that the charcoal curve need not lie on the reference. Converging the coupling removes the disagreement between the two blocks for the chosen time discretization. A coarse time step can still miss the early power peak or distort the transient. Time-step refinement remains necessary.

When does Newton help?

Picard asks the two calculations to take turns. Newton also uses information about how their answers change when a trial changes.

In this small model, the thermal calculation turns a guessed power $p$ into temperatures, and the neutron calculation turns those temperatures back into a power. Call that full round trip $\mathcal G(p)$. Agreement means

$$R(p)=p-\mathcal G(p)=0.$$

Picard takes the returned answer as its next guess. Newton uses a local linearization of a coupled residual to choose a correction. In our implementation the temperature and precursor variables can be eliminated, leaving a quadratic equation in power, so its derivative is inexpensive to calculate.

At the larger step sizes in this example, Newton reaches agreement while undamped Picard does not. If both converge to the same root, they produce the same implicit step. Newton's advantage here is reaching that answer reliably and in fewer updates, rather than defining a different physical solution.

For a large coupled model, that calculation is more demanding. Newton needs a suitable initial guess and a sufficiently accurate linear solve, and may need damping or other safeguards. Jacobian-free Newton–Krylov methods can obtain derivative information through perturbed residual evaluations, but those evaluations and preconditioning have a cost. The scalar iteration counts above are not production-code speedups. Published comparisons in TINTE discuss both fully implicit Picard and JFNK, including those implementation choices.

Under-relaxation or acceleration can also improve Picard. Choosing between them is an engineering decision based on robustness, accuracy and total computational work.

When is the extra conversation worth it?

It is worth investigating converged coupling when a time-step refinement study shows that lagged feedback materially changes a quantity you care about: peak power, fuel temperature, or the timing of a response. Rapid transients, strong feedback and additional interacting fields can make that issue more visible. A small reactor with flow participating in the feedback is one relevant application, but the same question applies to a large reactor.

If Picard converges quickly, it may already be enough. If it stalls or oscillates, compare relaxation, acceleration and Newton-type methods. Judge the comparison at matched tolerances, and count calls to the physics solvers or elapsed time as well as outer iterations.

Three checks keep the interpretation clear:

• Coupling convergence: do the fields satisfy the coupled discrete equations to the requested tolerance?
• Time and mesh refinement: do the quantities of interest stop changing materially as the discretization is refined?
• Physical validation: do the models represent the relevant phenomena and agree with suitable evidence?

These checks complement each other. A small residual answers the first question. It cannot answer the other two.

The model, if you want to check it

The formulas and parameters are included here so the experiment can be reproduced. They are deliberately simple enough that the implicit step also has an algebraic solution, which provides an independent check on the iterative solver.

Equations, units and numerical methods

Let $p=P/P_0$ be normalized fission power and $c$ the delayed-neutron precursor inventory normalized to its initial equilibrium value. A one-group point-kinetics approximation gives

$$\dot p=\frac{(\rho-\beta)p+\beta c}{\Lambda}, \qquad \dot c=\lambda(p-c).$$

This normalized form follows from the usual point-kinetics equations. The effective one-group parameters here are illustrative, not a fitted delayed-neutron dataset for a specific core.

Write $\theta_f$ and $\theta_c$ for departures from the initial average fuel and coolant temperatures, in kelvin. Subtracting the initial steady heat balance gives

$$C_f\dot\theta_f=P_0(p-1)-H(\theta_f-\theta_c),$$$$C_c\dot\theta_c=H(\theta_f-\theta_c)-W\theta_c,$$$$\rho=\rho_{\rm ext}+\alpha_f\theta_f+\alpha_c\theta_c.$$

$H$ is fuel-to-coolant conductance, and $W=\dot m c_p$ is an effective heat-removal coefficient for a mixed coolant node with fixed inlet temperature. The fuel-to-coolant transfer cancels when the energy balances are added: it moves heat internally. This model treats fission power as deposited heat and omits decay-heat dynamics. The coolant feedback coefficient represents a local temperature effect, including density changes at assumed fixed pressure; density is not an independent state here.

| Quantity | Illustrative value | | --- | --- | | Nominal thermal power $P_0$ | 100 MW | | Fuel heat capacity $C_f$ | $4/3$ MJ/K | | Coolant heat capacity $C_c$ | 20 MJ/K | | Conductance $H$ | $1/3$ MW/K | | Heat-removal coefficient $W$ | 2.5 MW/K | | Delayed fraction $\beta$ | 0.0065 | | Prompt generation time $\Lambda$ | $10^{-4}$ s | | Effective precursor decay constant $\lambda$ | 0.08 s$^{-1}$ | | Fuel / coolant feedback | −3 / −10 pcm/K | | Applied reactivity step | +50 pcm at $t=0$ |

Initially $p=c=1$ and $\theta_f=\theta_c=0$. The initial absolute temperatures do not enter these constant-coefficient perturbation equations. They would have to be specified, with material properties and pressure, in a physical plant model.

All blocks use backward Euler. The one-pass scheme evaluates reactivity with old temperatures, solves kinetics, then solves the two thermal balances with the new power. The converged scheme evaluates feedback at the new time level. Picard and Newton hold the same old state fixed while solving that step. In a reference-refinement check from 0.001 s to 0.0005 s, the maximum power difference at shared times over 60 s was about 0.045 percentage points of initial power. This checks numerical resolution, not physical validation.

After eliminating the thermal variables and $c$, the coupled equation has the form $Dp^2+Bp-A=0$. The charcoal curve uses Newton on that equation. Its derivative is analytic, with a safeguard against leaving the positive, evaluable branch. The reference uses the positive algebraic root at a 0.001 s step. Both iteration methods are checked against the same residual $|p-\mathcal G(p)|\leq10^{-10}$; that is a dimensionless numerical tolerance chosen for this demonstration.

For the first step, the signed Picard slope at the solution is about −0.091 at 0.25 s and −1.170 at 4 s. A magnitude below one gives local convergence under the usual smoothness assumptions; a magnitude above one makes the fixed point locally unstable for undamped Picard. This is a property of this iteration map at this state and step size. It is not a classification of SMR designs.

Why this is useful in practice

In coupled-code work, the useful skill is knowing what the exchanged variables mean, which time level they belong to, and how to tell an unconverged exchange from an inaccurate model. That is the perspective behind my work on implicit neutronic–thermal-hydraulic coupling with Enzo Dari, presented at ENIEF 2017 and published as a summary in Mecánica Computacional, XXXV(32), p. 1827.

The same questions appear in the physical and numerical models inside Criticality, a reactor simulator built for exploring behaviour. For an SMR study, the next step beyond this teaching example is to bring in the actual geometry, cross sections, flow model and validation cases—and keep asking whether the two calculations are describing the same state.

Try the browser-based nuclear reactor simulator to explore a broader plant model.

Working on a coupled reactor model? Let’s discuss the coupling assumptions and verification strategy.